- Email: [email protected]
A melt clusterization within the interfacial boundary layer and its hydrodynamics modelling at the microgravity semiconductor single crystal growth
ARTICLE IN PRESS
Journal of Crystal Growth 270 (2004) 329–339 www.elsevier.com/locate/jcrysgro
A melt clusterization within the interfacial boundary layer and its hydrodynamics modelling at the microgravity semiconductor single crystal growth Vladimir Ginkina, Andrey Kartavykhb,, Mikhail Zabudkoa a Institute for Physics and Power Engineering, Bondarenko sq.1, 249033 Obninsk, Kaluga region, Russia Institute of Chemical Problems for Microelectronics, Department of Space Materials Science, B. Tolmachevsky per.5, 119017 Moscow, Russia
b
Received 7 April 2004; accepted 28 June 2004 Communicated by M. Uwaha Available online 25 August 2004
Abstract The paper gives a brief analysis of current concepts on the process of melt ordering and structural self-organization at the temperature close to melting point, including that within the interface area when growing semiconductor single crystals. The experimental data presented confirm the existence of ordered multi-atomic growth units (clusters) in the melt under the conditions close to undisturbed equilibrium crystallization. The conclusion is justified on the prospects of studying the clustering processes within a melt in microgravity. A mathematical model of convective mass transfer is proposed as an independent investigation tool. This model includes the equations of hydrodynamics, impurity transfer and convective flow within the interface transient area. Materials properties and the temperature are presented as a single-value function of enthalpy. For the first time the melt structural model near the crystallization front considers the availability of clusters formation (CF) that cause resistance to the melt flow. This resistance of medium is described by means of specific double-phase coefficient introduction depending on the enthalpy value in each assigned mesh microvolume under calculation. Basing onto our previous space growth experiments aboard the Photon series satellites (J. Crystal Growth 205 (1999) 497; 223 (2001) 29), the necessary empirical model parameters are evaluated. The results of CF-model testing are demonstrated, including the description of real orbital and ground-based GaSb:In crystallization experiments. r 2004 Elsevier B.V. All rights reserved. PACS: 81.10.Fq; 81.10.Mx; 47.11.+j; 47.27.Te; 64.70.Dv Keywords: A1. Computer simulation; A1. Growth models; A1. Heat transfer; A1. Mass transfer; A1. Nucleation; A2. Microgravity conditions
Corresponding author. Tel.: +7-095-239-9063; fax: +7-095-953-8869.
E-mail address: [email protected] (A. Kartavykh). 0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.06.047
ARTICLE IN PRESS 330
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
1. Review of the problem Investigation of a substance structural transformations close to the melting point temperature ðT Þ at phase transitions is of the highest priority for understanding the fundamental processes and mechanisms of crystallization. In particular, there are no radical problems for the modern semiconductor material science either in grown single crystal structure/properties characterization or in the description and modelling of hydrodynamics and heat-mass transfer (HMT) in a melt bulk during crystallization. But the substance state in the narrow interfacial phase transition region, the HMT processes within the pre-crystallization layer still remain ‘‘blank points’’ to a significant extent. The difficulties in this study are aggravated by the fact that the transient area in the process of onground growth has small geometrical dimension and is permanently destructing (‘‘smoothing’’) from the melt side with the natural and artificially created convection flows. For instance, in Czochralski technique the typical interface boundary layer thickness is only just 0.1–0.5 mm. Gubenko [1] experimentally found that close to the liquidus temperature the doped metal and semiconductor melts do not follow the theory of strongly diluted (ideal) solutions and their properties cannot be described with the linear laws of Henry, Raul, etc. He has formulated the theory of melt structural self-organization that contains the following key statements: (1) along with free atoms in a melt there exist the ordered atom groups (OAG or clusters), (2) clusters have statistically localizedin-time, directed bonds between neighbouring atoms and the definite dynamically-equilibrium near order, (3) melt overheating above the liquidus line results in continuous cluster destruction and in weakening of interatomic interaction forces. That is accompanied in macroscale with melt viscosity and density decreasing as well as changing of other structure-dependent properties, (4) in the crystallization course the local heterogeneous equilibrium is established in the conjugate phase interface, with chemical potentials of components in conjugate phase layers being equal. So any changes in the melt transient layer structure result in local variations of distribution coefficients of impurities
presenting here, and in structural changes of adjacent solid phase layer. Description of the microscopic nature of the liquid state is one of the most difficult problems in condensed matter physics in general. Owing to the dynamical nature of the liquid state, it is not possible to discuss a particular microscopic structure; only ensemble averages can be specified. Such averages can be performed via well-crafted molecular dynamics simulations: the length of the simulation, the size of the ensemble and the nature of the interatomic forces must all be carefully analysed. The molecular dynamics approach is now one of the most resulting numerical modelling method of ‘‘cluster’’ mechanisms of pre-crystallization state evolution in microscale. Significant methodological progress has been recently achieved in this direction. In paper [2] when modelling GaP single crystal growth from stoichiometric melt it was shown that the total cluster concentration in transient layer in A3B5 compounds can achieve (2–3)1020 cm3 , including 1019 cm3 fraction of the biggest clusters, whose number of associated atoms being equal to 50–75. This estimation agrees well with typical values of dopants concentration in near-the-front diffusion layers and testifies that dopant atoms can really serve as initial cluster-originating centers. However it should be noted that calculation results obtained with the molecular dynamics methods depend strongly on the choice of the specific type of atom interaction potentials in the melt and crystal. The Lennard-Johns and Born-Mayer potentials are usually used a priori. Mathematical expressions for potential description in an explicit form have a number of uncertain empirical parameters, thus resulting in unjustified simplification during crystallization process modelling. In view of that it is necessary to thoroughly check the consistency and agreement of calculation and experimental data. The most of disadvantages above can be avoided using ab initio pseudopotential-density functional method that is highly accurate and well tested for semiconductors in the solid state, but has only recently been applied to liquids. Within the pseudopotential method only the valence electrons are explicitly treated, and within the density functional theory exchange and
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
correlation terms are mapped onto an effective one-electron potential. Good review of recent progress in that direction is done in paper of Chelikowsky et al. [3]. It is shown numerically that A3B5 and A2B6 compound melts, being in adiabatic conditions at the temperature near T , are able to form both anion and cation clusters having chain- and even crown-like pseudo-molecular dynamical microstructure. Strong and complex melt ordering in Ge, Si, GaAs, CdTe and GeTe, especially at the supercooled metastable state, is confirmed also in [3] with direct X-ray and neutron diffraction measurements. Some experimental facts confirm the existence and stability of structural self-organization effects in semiconductor melt under conditions of extremely weak convection. Ostrogorsky and Dutta [4] used a specially developed fine-controlled submerged heater method and pre-synthesized GaSb and InAs (x=2–5 at.%) as a load, and managed to get a set of bulk single crystals of quasi-binary composition ðGaSbÞ1x ðInAsÞx , which drastically differed from usual quadruple solid solutions Gað1xÞ Inx Asx Sbð1xÞ in their fundamental properties. For example, quasi-binary crystals have a smaller band gap width with the equal atomic composition. Thus experimentally it has been proved that double-atomic molecular dipole In–As clusters exist in GaSb melt, and they may built-in into the crystal lattice as an independent structural growth unit. During the study of doped Ge and Si crystallization by precision differential thermal analysis (DTA) technique it is possible to observe weak thermal effects related to clusterization at the temperature very close to that of phase transition. In papers [5,6] the melt restructuring temperature (and thus, the cluster formation specific energy) was shown to be different with the presence of B, P and Sb dopants in Si. The similar extreme temperature points reflecting in macroscale the process of melt self-ordering are accompanied with the significant changes in thermophysical properties of silicon melt. In this respect the best known are, probably, results of Kimura Meta Melt Project [7–9]. An anomalous density jump with the magnitude depending on initial doping nature was found near the melting point of Si. Surface
331
tension and viscosity also showed an anomalous behaviour. However, the recent detail investigation of Si physical characteristics, using the modern technique of electromagnetic levitation of molten drop, does not observe above-mentioned effects even under deep supercooling [10,11]. The reasons of that discordance are not definitely clear. One can note only that melt stirring inside the levitating drop is extremely intensive due to the strong external actions of supporting electromagnetic field, IR or RF-microwave surface heating etc. That may disturb completely the fine melt ordering process. As the fundamental effects of structure ordering in the transient layer are extremely sensitive to any energy-inserting influence—thermal, gravitational or mechanic (convection, stirring, vibration), electromagnetic etc., the microgravity state realizing aboard a space vehicle must be the ideal experimental medium for their manifestation and study. During crystallization in conditions maximum close to weightlessness and under growth near to diffusion-controlled mode, the transient region thickness in semiconductor melt can achieve 3–6 mm [12–14]. The purposeful study of the melts state in space has not been performed yet. However there are some experimental facts, which indirectly confirm the existence of strong structuring effects in a melt in conditions of orbital space flight. For instance, the liquid-phase diffusion coefficients (D) of metallic alloys components measured from space experimental data are known to be always 2–5 times lower relatively to those in normal gravity ðg0 Þ. When studying the diffusion processes in narrow capillaries aboard the Space Shuttle and Mir orbiters, Smith et al. [15,16] discovered the abnormal linear temperature dependence D T for all molten systems under study, including Sb–In, Sb–Ga, Bi–Sb etc. (11 systems overall were tested including more than 200 single-point experiments). In contrast, in ground-based reference tests the usual Arrhenius relation D expð1=TÞ was observed. Let us note that D is extremely sensitive parameter to the structure of medium where diffusion occurs. Evidently, diffusion in space occurred in undisturbed ordered melts whose structure significantly changed under the temperature. It is significant
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
332
that superimposing of controlled forced vibration (0.1 Hz, 4 103 g0 ) onto diffusion cell aboard orbiter led to another surprising relationship revealing D T 2 [15]. That result is still waiting for its detail explanation. However, it may indirectly indicative of the mechanically inserted reconstruction of a melt structure. The data like above-mentioned serve as the objective prerequisite for theory development of ordered transient layer with the attraction of space-based crystallization experiments. Macroscopic changes of melt structure-dependent characteristics in the front-adjoining layer, which reflect the dynamics of cluster formation at the micro-level, give the grounds for developing hydrodynamics models that can become an additional and independent exploration tool. The idea of developing and testing such a model as applied to the process of ground and orbital semiconductor crystallization is partially implemented by the authors in this work. The model development makes it possible to consider the fundamental peculiarities of crystal melt growth processes, in particular, to correctly calculate the HMT both in the core liquid region and in interface vicinity.
2. Mathematical model of melt crystallization The laminar flow in the melt is described by the Navier–Stokes equations in their Boussinesq approximation form: @h ~ r þ V rh ¼ rl rT; ð1Þ @t
r0
! ~ @V ~ rÞ V ~ ¼ rp þ ðV @t ~r0 ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ ~ þ rm rV g;
ð2Þ
In the solidified sections of the calculated area ~ ¼ 0, and Eqs. (1)–(4) will take the form: V r
@h ¼ rl rT; @t
@C ¼ rD rC: @t
ð5Þ
ð6Þ
In Eqs. (1)–(6) r; m; l; D are density, viscosity, thermal conductivity and diffusion coefficients respectively; bT ; bC are thermal expansion and concentration expansion coefficients; h; T; C are deflections of enthalpy, temperature and concentration from certain stationary values h0 ; T 0 ; C 0 ; r0 is the melt density at T ¼ T 0 ; C ¼ C 0 ; and ~ gthe gravity acceleration. For simplicity’s sake h; T; C will subsequently be referred to as enthalpy, temperature and concentration. Eqs. (1)–(6) are supplemented by initial and boundary conditions. Let us break down the region over which the solution is to be computed into meshes and average the enthalpy values in each of the mesh cells. Let us express the average mesh temperature T in terms of enthalpy h: 8 ðh h1 Þ > > Tn þ for hph1 ; > > > C ps ðhÞ < T ¼ Tn for h1 ohoh1 þ x; > > > h h x 1 > > for hXh1 þ x; : Tn þ C pl ðhÞ where C ps ðhÞ is heat capacity of the solid phase of a substance; C pl ðhÞ is the heat capacity for the liquid phase; h1 is the maximum value of enthalpy for the solid phase, x is phase transition latent heat, T n is crystallization temperature. The state of a substance into mesh cells with enthalpy values h1 ohoh1 þ x is intermediate (transient) from thermodynamics viewpoint between solid and liquid. For pure substances T ¼ T 0 ¼ const. For substances containing an impurity T n is dependent on C: T n ¼ T n0 þ mC;
~ ¼ 0; rV
ð3Þ
@C ~ þ V rC ¼ rD rC: @t
ð4Þ
where m is the slope of the solidus line in the phase diagram. Let us use L and S to denote the fraction of liquid and solid phase respectively in the inter-
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
mediate-state region h1 ohoh1 þ x: 8 0 for hph1 ; > > x > : 1 for hXh1 þ x;
333
and porosity coefficient K, and is opposed to the ~ in direction: vector V ~ @V ~rÞ V ~ ¼ 1 rp þ n DV ~ þ ðV @t r0 ~ nV K; ð7Þ S where n ¼ rm ; S is square of mesh section normally 0 ~, a is an empirical ‘‘clusterizato velocity vector V tion constant’’ which together with K describes the double-phase medium nature. In physical sense a will depend firstly on the shape and averaged size of clusters, i.e. on the surface/volume ratio of solid particles putting up resistance to the flow motion. When K ¼ 1 (in solid state) Eq. (7) changes into ~ ¼ 0, and when K ¼ 0 it goes over to usual V hydrodynamics Eq. (2). Impurity transfer in meshes where the substance is in transient state will be described by the following Equation [18]: ga ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ~
Let us assume that the quantities r; l; C p ; D are constants for the solid and liquid phases of the substance. Meshes with transient enthalpy values will also have intermediate values of r; l; C p ; D which we will assume to be linear with h. If we use X s to denote a parameter value for the solid phase and X L for the liquid phase, the value of X for the intermediate state can be obtained from the formula: X ¼ X S L þ X L S ; T¼
ðh L xÞ ; Cp
where the parameter X is r; l; C p ; D. Let us assume that the intermediate melt region is filled with clusters, i.e. population of atoms (molecules) associated in solid structures. The size and number of clusters can be arbitrary. The clusters’ heat movement velocity is significantly lower than that for single melt atoms because they have much bigger sizes. The results of moleculardynamic calculations confirm this assumption. For example, the estimation of cluster diffusion coefficient with the number of associated atoms of about 70 in GaP melt given in [2] shows the value D ¼ 8:6 107 cm2 =s, whereas for single phosphorus atoms D ¼ 8:1 105 cm2 =s. The difference is seen consists 2 orders of magnitude. The given value allows the multiatomic cluster movement velocity to be neglected within the model. Then the melt flow in the medium that is filled not tightly with motionless clusters can be presented as a liquid flow in a porous double-phase body [17]. The melt flow in meshes where the substance is in intermediate state will be described using the porous solid approximation with the double-phase (porosity) coefficient K ¼ LS . To be able to make such a description let us make Eq. (2) include a resisting force due to conditions of porosity. This force varies as the melt viscosity, as flow velocity
@C ~ @S þ V rC ¼ rD rC þ ð1 K p ÞC; ð8Þ @t @t where K p o1 is impurity segregation coefficient. Finally, the mathematical model of heat-mass transfer for the whole area, where the crystallization process takes place, includes Eqs. (1), (3), (7), (8) in which the temperature T and values of r; l; C p ; D; L ; S are single-valued and continuous one-variable functions of the enthalpy h. The hallmark of the proposed model is that there is no localized phase transition interface. This boundary is fluctuated (‘‘smudged’’) forming a region of meshes whose enthalpy values refer to an intermediate state of the substance. Nevertheless the crystallization front velocity can be defined as the rate at which the proportion of the solid/liquid phases of the substance changes. This approach allows one to describe impurity segregation in the process of crystallization together with the impurity diffusion–convection transfer in the melt. As far as the hydrodynamics equations are concerned the adhesion condition for the flow velocity is laid down at the boundary between the solid phase and the intermediate state of the substance being crystallized. And the additional resistance to the flow in the intermediate state region is described in terms of porous body
L
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
334
approximation in which the porosity coefficient is determined by the solid-to-liquid phase ratio.
3. Computation algorithm Let us introduce a time grid with variable pitch t and use an implicit discretization scheme. This yields the following equation system: ~ rTðhÞ þ h ¼ h ; raðhÞ rTðhÞ þ V Cpt Cpt ~ þ ðV ~ rÞ V ~þ n DV ¼
ð9Þ
~ ~ V nV þa K t S
~ V 1 rp ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ~ gþ ; t r0 ð10Þ
~ ¼ 0; rV
ð11Þ
S ð1 K p ÞC þ L C t S ð1 K p ÞC þ L C ; ð12Þ ¼ t ~ where the symbols V ; S ; C refer to the values of those quantities at the preceding point of time. The system of equations (9)–(12) is linearized and solved by the method proposed in [19] and described in [20]. This method consists from the Patankar method [21], exponential transformation method and the conjugate gradient method with preconditioning by the incomplete factorisation method [22]. ~rC þ rDrC þ V
melting point). Meanwhile, the given model does not consider true (molecular) viscosity of the melt (n), but some effective (integral) viscosity (n ) of substance having been in transient state between fluid and solid. Let us call this model of integral viscosity increase one (IVI-model). The cause of integral viscosity increase near the interface, as it was mentioned above, could be the formation of high-drag low-mobile multiatomic clusters resisting the melt flow within the transient area (just is cluster formation or CF-model). In the present work the last term in the Eq. (7) describes the resistive force like as Darcy law for the fluids current in a porous medium. In CF-model the melt viscosity is assumed to be equal to its true value, and is either constant or slightly depending function of temperature. For establishing the quantitative relationship between IVI- and CF-models it is necessary to evaluate the empirical parameter a in the last term of the Eq. (7). For this purpose we have used for a long-time-knowing fundamental effect of compositional melt stratification in hermetic ampoule under the action of weak diametrical gravitational force (see theory e.g. in [24]). The action of stationary axis-perpendicular microacceleration leads to formation of monotonously increasing diametrical concentration profile of dopant, having K p o1, in the growing crystal along the gravity vector direction. Such the effect was observed repeatedly by various authors in the course of space-grown crystal study (see, e.g. [25,26]). As the test basis we have chose the paper of Duffar et al. [26], where GaSb:In crystallization experiments were performed by Bridgman technique both aboard the EURECA orbiter and on-ground. For sketch depicted in Fig. 1 the parametric
4. Empirical parameter evaluation In paper [23] the mathematical model is proposed allowing describe adequately the experimentally revealed anomalous distribution of impurity in low-doped semiconductor crystals grown under conditions of Photon series satellites orbital flight. In the foundation of this model the insight lays on the sharp increasing of melt viscosity in growth interface vicinity within the transient layer limiting by [T þ 10K; T ] isotherms (where T is
Fig. 1. The calculation scheme (indium dopant concentrations are given for initial GaSb charge parts before the melting). Growth rate—0:1 mm=s, initial axial thermal gradient—20 K/ cm, final axial thermal gradient-30 K/cm.
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339 3.1 3.0 2.9
∆C,× 1017at/cc
2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0 10.0
(1/v)⋅(dv*/dT), K-1
Fig. 2. Relationship between diametrical inhomogeneity ratio of In-dopant distribution in GaSb crystal, and integral viscosity parameter in transient boundary layer. Calculation is done for axis-perpendicular stationary microgravity level of 104 g0 , after 106 s since solidification beginning.
335
according both IVI- and CF-models. In the former the increase of integral viscosity was varied up to 100 times within transient area, and in the latter one the empirical clusterization constant a was modified. In both cases the resulting concentration difference (DC) of In impurity was calculated between the edges of crystal diameter. Fig. 2 shows the DC dependence on relative thermal gradient of integral viscosity ð1=nÞ ðdn =dTÞ, obtained by IVI-model. Fig. 3 shows the similar dependence of a, for which both models (IVI- and CF-) give identical diametrical non-uniformity DC in impurity distribution. The good agreement between calculated and experimental results was achieved in the work [23] when ð1=nÞ ðdn =dTÞ ¼ 6–7. These values correspond in Fig. 3 to a ð6:2–7:0Þ 1017 estimation. Assuming as the first approximation that quite conservative CF-model is valid for the most of semiconductors, in the current work we used for orbital crystallization process simulation the a ¼ 7 1017 value.
8.0
5. Model testing
α × 1017
6.0
In order to test the computer code, the numerical simulation was performed of the real crystallization experiments [26] within the frame of proposed CF-model both in space and on-ground. The calculation scheme is given in Fig. 1. Calculation was carried out in 2-D (X,Y) geometry. The values of physical parameters used in calculations are given in Table 1.
4.0
2.0
Table 1 Physical parameters of GaSb used in calculations
0.0 2.0
4.0
6.0 (1/v)⋅(dv*/dT),
8.0
10.0
K-1
Fig. 3. Relationship between the clusterization constant a and integral viscosity parameter in transient boundary layer. Calculation is done for the same conditions as that of Fig. 2.
calculations were made during 106 s at residual gravity vector direction of 104 g0 perpendicularly to growth axis. Calculations were performed
Melting point, K Density of the melt, kg/m3 Density of the solid, kg/m3 Thermal conductivity, Wt/(m K) Thermal capacity, J/(kg K) Kinematic viscosity, m2/s Diffusion coefficient of In in GaSb melt, m2/s Segregation coefficient of In in GaSb Temperature-expansion coefficient, K 1 Prandtl number
T rm rs lm Cp n D
985 6010 5600 8 175 0:38 106 1:2 108
Kp bT Pr
0.15 2:0 104 0.05
ARTICLE IN PRESS 336
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
Fig. 4 shows the calculated and experimental distributions of In isovalent dopant concentration along the reference crystal axis for the conditions of earth gravity. The calculation is made with a ¼ 0 in the Eq. (7). The agreement between the calculated results and experimental data is rather satisfactory. It obviously points to that the melt layer containing clusters has very small width and its influence on the convective heat and masstransfer in earth conditions can be neglected. Quite another matter is crystal grown in microgravity. The good agreement of calculated
and experimental axial distributions of In was achieved in Fig. 5 under zero gravity conditions only with the a ¼ 7 1017 value application. That indicates a really low level of convection and the existence of highly pronounced transient boundary layer during the experiments aboard the EURECA orbiter. In paper [26] the impurity radial distributions had different slopes for different cross-sections.
Fig. 4. Longitudinal In concentration distribution in GaSb for the conditions of earth gravity: 1—measured [26], 2—calculated using the clusterization parameter ratio a ¼ 0.
Fig. 5. Longitudinal In concentration distribution in GaSb for the conditions of zero gravity: 1—measured [26], 2—calculated using the clusterization parameter ratio a ¼ 7 1017 .
Fig. 6. Distribution of In impurity concentration in C 1019 at=cm3 (a) and equal In concentration lines in the longitudinal crystal section (b).
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
337
Fig. 7. Spatial distribution of longitudinal and transverse components of melt flow rate vector and melt flow rate fields. (a)—at the moment of 5 104 s, (b)—at the moment of 105 s.
The authors of [26] suppose that it indicates the change in the direction of residual gravity vector. In order to check this assumption the calculation was made of GaSb crystallization at g ¼ 104 g0 when the residual gravity vector assumed been perpendicular to the growth axis during the first 5 104 s of crystallization, and during the next 5104 s it had the opposite direction. Fig. 6 shows the spatial distributions of In concentration in the solid GaSb being crystallized, as well as concentration isolines obtained in result of this calcula-
tion according the CF-model. The presented data show that isolines firstly are curving towards one side, then are rectified, and finally are curving towards the other side. Thus confirming the model performance at the given assumption on the existence of residual gravitation and precession of its vector during solidification. Fig. 7 shows spatial distributions of flow rate vector components for axes X and Y at the time moments of 5 104 and 105 s in order to illustrate the flow nature in the considered case. These data
ARTICLE IN PRESS 338
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
testify that the vortex reverses its direction when the direction of residual gravitation vector changes. The maximum value of flow rate module in the first case was equal to 0:34 105 m=s, and in the second case 0:39 105 m=s.
6. Conclusion The model presented in the paper allows the description of melt crystallization process in view of moving interface, convection heat and mass transfer and impurity segregation. For the first time the hydrodynamics model considers the structural pre-crystallization state of the boundary layer in the interface area. The transient region is firstly described as a double-phase medium containing three-dimensional solid phase clusters with a melt together. The presented model does not simulate the quantity and the size of the clusters near the growth interface, and takes into account only the influence of cluster-containing transient area on hydrodynamics of melt and distribution of impurity. The objective of the current work was the estimation of this influence. The calculations show that this influence in earth conditions is insignificant and it can be only detected when the external forces are small, e.g. in condition of microgravity. The specific feature of the model consists in the description of clusterization processes and those of substance crystallization in enthalpy (thermodynamic) variables, thus giving further the possibility to cross-over the results with the data obtained with other independent approaches, i.e. physical chemistry, molecular dynamics, etc. So there is a potential for further development and improvement of model representations. The next step proposes a quantitative description of clusters distribution by their sizes using molecular-dynamic approach and establishing of these data connection with local values of enthalpy. In the suggested version of model the semiempirical consideration was only given for the atmost case of melts structure ordering, i.e. for the formation of solid-state associated particles. Nevertheless the model gives quite adequate description of crystallization system response to
the change of key external conditions—temperature, crystallization rate, magnitude and direction of microacceleration vector, etc. The results of Indoped GaSb crystallization test calculations show a good agreement with the space experimental data. The most pronounced effect from the use of model can be expected for description of crystallization processes performed in the conditions of reduced gravity. That is the case when the boundary segregation phenomena near the crystallization front the most prominently affects the quality of the crystal growing.
Acknowledgements We would like to acknowledge support from the Russian Foundation for Basic Research, grant No. 03-02-16282.
References [1] A.Ya. Gubenko, High purity substances 4 (1994) 110 (in Russian). [2] S.V. Kotov, A.R. Lyutikov, Yu.P. Khukhryansky, et al., Tech. Phys. Lett. 28 (2002) 584. [3] J.R. Chelikowsky, J. Derby, V. Godlewsky, et al., J. Phys.: Condens. Matter 13 (2001) R817. [4] A.G. Ostrogorsky, P.S. Dutta, J. Crystal Growth 197 (1999) 1. [5] A.Ya. Gubenko, Microelectronic Eng. 45 (1999) 161. [6] A.Ya. Gubenko, Crystallogr. Rep. 46 (2001) 81. [7] H. Sasaki, E. Tokizaki, K. Terashima, S. Kimura, Japan J. Appl. Phys. 33 (1994) 3803. [8] S. Kawanishi, H. Sasaki, K. Terashima, S. Kimura, Japan J. Appl. Phys. 34 (1995) L1509. [9] S. Kimura, K. Terashima, J. Crystal Growth 180 (1997) 323. [10] M. Langen, T. Hibiya, M. Eguchi, I. Egry, J. Crystal Growth 186 (1998) 550. [11] W.-K. Rhim, K. Ohsaka, J. Crystal Growth 208 (2000) 313. [12] A.V. Kartavykh, E.S. Kopeliovich, M.G. Milvidskii, V.V. Rakov, Crystallogr. Rep. 43 (1998) 1075. [13] A.V. Kartavykh, E.S. Kopeliovich, M.G. Mil’vidskii, V.V. Rakov, Microgravity Sci. Technol. XII (1999) 16. [14] A.V. Kartavykh, E.S. Kopeliovich, M.G. Mil’vidskii, V.V. Rakov, J. Crystal Growth 205 (1999) 497. [15] R.W. Smith, X. Zhu, M.C. Tunnicliffe, et al., Proceedings of the International Conference ‘‘Spacebound-2000’’.
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
[16]
[17]
[18] [19]
Banff, Canada, 2000. p. 96 (also www.space.gc.ca/asc/pdf/ 2000 SMITH1.PDF). R.W. Smith, X. Zhu, M. C. Tunnicliffe, et al., Proceedings of the First International Symposium on Microgravity Research and Application in Physical Sciences and Biotechnology, Sorrento, Italy, Vol. II. 2000. ESA publ. SP-454, p. 887. V. Ginkin, O. Naumenko, M. Zabudko, et al., in: Proceedings of the International Symposium on Advances in Computational Heat Transfer, (CHT-04). Midnatsol Cruise, Norway, 2004, p. 226. V. Timchenko, P.Y.P. Chen, D.G. de Vahl, et al., Int. J. Heat and Mass Transfer 43 (2000) 963. V.P. Ginkin, in: Proceedings of the Second International Symposium on Advances in Computational Heat Transfer, Palm Cove, Australia, 2001, p. 1161.
339
[20] V.P. Ginkin, S.M. Ganina, Math. Model. 15 (2003) 53 (in Russian). [21] S.V. Patankar, D.B. Spalding, J. Heat Mass Transfer 15 (1972) 1787. [22] V.P. Ginkin, A.V. Kulik, O.M. Naumenko, in: Proceedings of the 4th International Conference on Supercomputing in Nuclear Applications (SNA 2000), Tokyo, Japan, 2000, p. 78. [23] V.K. Artemyev, V.I. Folomeev, V.P. Ginkin, et al., J. Crystal Growth 223 (2001) 29. [24] S.A. Nikitin, V.I. Polezhaev, A.I. Fedyushkin, J. Crystal Growth 52 (1981) 471. [25] V.S. Zemskov, V.N. Kubasov, I.N. Belokurova, et al., Apollo-Soyuz Test Project. Summary Science Report. V.1. Washington, NASA edit., 1977, p. 539–553. [26] T. Duffar, M.D. Serrano, C.D. Moore, et al., J. Crystal Growth 192 (1998) 63.
Journal of Crystal Growth 270 (2004) 329–339 www.elsevier.com/locate/jcrysgro
A melt clusterization within the interfacial boundary layer and its hydrodynamics modelling at the microgravity semiconductor single crystal growth Vladimir Ginkina, Andrey Kartavykhb,, Mikhail Zabudkoa a Institute for Physics and Power Engineering, Bondarenko sq.1, 249033 Obninsk, Kaluga region, Russia Institute of Chemical Problems for Microelectronics, Department of Space Materials Science, B. Tolmachevsky per.5, 119017 Moscow, Russia
b
Received 7 April 2004; accepted 28 June 2004 Communicated by M. Uwaha Available online 25 August 2004
Abstract The paper gives a brief analysis of current concepts on the process of melt ordering and structural self-organization at the temperature close to melting point, including that within the interface area when growing semiconductor single crystals. The experimental data presented confirm the existence of ordered multi-atomic growth units (clusters) in the melt under the conditions close to undisturbed equilibrium crystallization. The conclusion is justified on the prospects of studying the clustering processes within a melt in microgravity. A mathematical model of convective mass transfer is proposed as an independent investigation tool. This model includes the equations of hydrodynamics, impurity transfer and convective flow within the interface transient area. Materials properties and the temperature are presented as a single-value function of enthalpy. For the first time the melt structural model near the crystallization front considers the availability of clusters formation (CF) that cause resistance to the melt flow. This resistance of medium is described by means of specific double-phase coefficient introduction depending on the enthalpy value in each assigned mesh microvolume under calculation. Basing onto our previous space growth experiments aboard the Photon series satellites (J. Crystal Growth 205 (1999) 497; 223 (2001) 29), the necessary empirical model parameters are evaluated. The results of CF-model testing are demonstrated, including the description of real orbital and ground-based GaSb:In crystallization experiments. r 2004 Elsevier B.V. All rights reserved. PACS: 81.10.Fq; 81.10.Mx; 47.11.+j; 47.27.Te; 64.70.Dv Keywords: A1. Computer simulation; A1. Growth models; A1. Heat transfer; A1. Mass transfer; A1. Nucleation; A2. Microgravity conditions
Corresponding author. Tel.: +7-095-239-9063; fax: +7-095-953-8869.
E-mail address: [email protected] (A. Kartavykh). 0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.06.047
ARTICLE IN PRESS 330
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
1. Review of the problem Investigation of a substance structural transformations close to the melting point temperature ðT Þ at phase transitions is of the highest priority for understanding the fundamental processes and mechanisms of crystallization. In particular, there are no radical problems for the modern semiconductor material science either in grown single crystal structure/properties characterization or in the description and modelling of hydrodynamics and heat-mass transfer (HMT) in a melt bulk during crystallization. But the substance state in the narrow interfacial phase transition region, the HMT processes within the pre-crystallization layer still remain ‘‘blank points’’ to a significant extent. The difficulties in this study are aggravated by the fact that the transient area in the process of onground growth has small geometrical dimension and is permanently destructing (‘‘smoothing’’) from the melt side with the natural and artificially created convection flows. For instance, in Czochralski technique the typical interface boundary layer thickness is only just 0.1–0.5 mm. Gubenko [1] experimentally found that close to the liquidus temperature the doped metal and semiconductor melts do not follow the theory of strongly diluted (ideal) solutions and their properties cannot be described with the linear laws of Henry, Raul, etc. He has formulated the theory of melt structural self-organization that contains the following key statements: (1) along with free atoms in a melt there exist the ordered atom groups (OAG or clusters), (2) clusters have statistically localizedin-time, directed bonds between neighbouring atoms and the definite dynamically-equilibrium near order, (3) melt overheating above the liquidus line results in continuous cluster destruction and in weakening of interatomic interaction forces. That is accompanied in macroscale with melt viscosity and density decreasing as well as changing of other structure-dependent properties, (4) in the crystallization course the local heterogeneous equilibrium is established in the conjugate phase interface, with chemical potentials of components in conjugate phase layers being equal. So any changes in the melt transient layer structure result in local variations of distribution coefficients of impurities
presenting here, and in structural changes of adjacent solid phase layer. Description of the microscopic nature of the liquid state is one of the most difficult problems in condensed matter physics in general. Owing to the dynamical nature of the liquid state, it is not possible to discuss a particular microscopic structure; only ensemble averages can be specified. Such averages can be performed via well-crafted molecular dynamics simulations: the length of the simulation, the size of the ensemble and the nature of the interatomic forces must all be carefully analysed. The molecular dynamics approach is now one of the most resulting numerical modelling method of ‘‘cluster’’ mechanisms of pre-crystallization state evolution in microscale. Significant methodological progress has been recently achieved in this direction. In paper [2] when modelling GaP single crystal growth from stoichiometric melt it was shown that the total cluster concentration in transient layer in A3B5 compounds can achieve (2–3)1020 cm3 , including 1019 cm3 fraction of the biggest clusters, whose number of associated atoms being equal to 50–75. This estimation agrees well with typical values of dopants concentration in near-the-front diffusion layers and testifies that dopant atoms can really serve as initial cluster-originating centers. However it should be noted that calculation results obtained with the molecular dynamics methods depend strongly on the choice of the specific type of atom interaction potentials in the melt and crystal. The Lennard-Johns and Born-Mayer potentials are usually used a priori. Mathematical expressions for potential description in an explicit form have a number of uncertain empirical parameters, thus resulting in unjustified simplification during crystallization process modelling. In view of that it is necessary to thoroughly check the consistency and agreement of calculation and experimental data. The most of disadvantages above can be avoided using ab initio pseudopotential-density functional method that is highly accurate and well tested for semiconductors in the solid state, but has only recently been applied to liquids. Within the pseudopotential method only the valence electrons are explicitly treated, and within the density functional theory exchange and
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
correlation terms are mapped onto an effective one-electron potential. Good review of recent progress in that direction is done in paper of Chelikowsky et al. [3]. It is shown numerically that A3B5 and A2B6 compound melts, being in adiabatic conditions at the temperature near T , are able to form both anion and cation clusters having chain- and even crown-like pseudo-molecular dynamical microstructure. Strong and complex melt ordering in Ge, Si, GaAs, CdTe and GeTe, especially at the supercooled metastable state, is confirmed also in [3] with direct X-ray and neutron diffraction measurements. Some experimental facts confirm the existence and stability of structural self-organization effects in semiconductor melt under conditions of extremely weak convection. Ostrogorsky and Dutta [4] used a specially developed fine-controlled submerged heater method and pre-synthesized GaSb and InAs (x=2–5 at.%) as a load, and managed to get a set of bulk single crystals of quasi-binary composition ðGaSbÞ1x ðInAsÞx , which drastically differed from usual quadruple solid solutions Gað1xÞ Inx Asx Sbð1xÞ in their fundamental properties. For example, quasi-binary crystals have a smaller band gap width with the equal atomic composition. Thus experimentally it has been proved that double-atomic molecular dipole In–As clusters exist in GaSb melt, and they may built-in into the crystal lattice as an independent structural growth unit. During the study of doped Ge and Si crystallization by precision differential thermal analysis (DTA) technique it is possible to observe weak thermal effects related to clusterization at the temperature very close to that of phase transition. In papers [5,6] the melt restructuring temperature (and thus, the cluster formation specific energy) was shown to be different with the presence of B, P and Sb dopants in Si. The similar extreme temperature points reflecting in macroscale the process of melt self-ordering are accompanied with the significant changes in thermophysical properties of silicon melt. In this respect the best known are, probably, results of Kimura Meta Melt Project [7–9]. An anomalous density jump with the magnitude depending on initial doping nature was found near the melting point of Si. Surface
331
tension and viscosity also showed an anomalous behaviour. However, the recent detail investigation of Si physical characteristics, using the modern technique of electromagnetic levitation of molten drop, does not observe above-mentioned effects even under deep supercooling [10,11]. The reasons of that discordance are not definitely clear. One can note only that melt stirring inside the levitating drop is extremely intensive due to the strong external actions of supporting electromagnetic field, IR or RF-microwave surface heating etc. That may disturb completely the fine melt ordering process. As the fundamental effects of structure ordering in the transient layer are extremely sensitive to any energy-inserting influence—thermal, gravitational or mechanic (convection, stirring, vibration), electromagnetic etc., the microgravity state realizing aboard a space vehicle must be the ideal experimental medium for their manifestation and study. During crystallization in conditions maximum close to weightlessness and under growth near to diffusion-controlled mode, the transient region thickness in semiconductor melt can achieve 3–6 mm [12–14]. The purposeful study of the melts state in space has not been performed yet. However there are some experimental facts, which indirectly confirm the existence of strong structuring effects in a melt in conditions of orbital space flight. For instance, the liquid-phase diffusion coefficients (D) of metallic alloys components measured from space experimental data are known to be always 2–5 times lower relatively to those in normal gravity ðg0 Þ. When studying the diffusion processes in narrow capillaries aboard the Space Shuttle and Mir orbiters, Smith et al. [15,16] discovered the abnormal linear temperature dependence D T for all molten systems under study, including Sb–In, Sb–Ga, Bi–Sb etc. (11 systems overall were tested including more than 200 single-point experiments). In contrast, in ground-based reference tests the usual Arrhenius relation D expð1=TÞ was observed. Let us note that D is extremely sensitive parameter to the structure of medium where diffusion occurs. Evidently, diffusion in space occurred in undisturbed ordered melts whose structure significantly changed under the temperature. It is significant
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
332
that superimposing of controlled forced vibration (0.1 Hz, 4 103 g0 ) onto diffusion cell aboard orbiter led to another surprising relationship revealing D T 2 [15]. That result is still waiting for its detail explanation. However, it may indirectly indicative of the mechanically inserted reconstruction of a melt structure. The data like above-mentioned serve as the objective prerequisite for theory development of ordered transient layer with the attraction of space-based crystallization experiments. Macroscopic changes of melt structure-dependent characteristics in the front-adjoining layer, which reflect the dynamics of cluster formation at the micro-level, give the grounds for developing hydrodynamics models that can become an additional and independent exploration tool. The idea of developing and testing such a model as applied to the process of ground and orbital semiconductor crystallization is partially implemented by the authors in this work. The model development makes it possible to consider the fundamental peculiarities of crystal melt growth processes, in particular, to correctly calculate the HMT both in the core liquid region and in interface vicinity.
2. Mathematical model of melt crystallization The laminar flow in the melt is described by the Navier–Stokes equations in their Boussinesq approximation form: @h ~ r þ V rh ¼ rl rT; ð1Þ @t
r0
! ~ @V ~ rÞ V ~ ¼ rp þ ðV @t ~r0 ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ ~ þ rm rV g;
ð2Þ
In the solidified sections of the calculated area ~ ¼ 0, and Eqs. (1)–(4) will take the form: V r
@h ¼ rl rT; @t
@C ¼ rD rC: @t
ð5Þ
ð6Þ
In Eqs. (1)–(6) r; m; l; D are density, viscosity, thermal conductivity and diffusion coefficients respectively; bT ; bC are thermal expansion and concentration expansion coefficients; h; T; C are deflections of enthalpy, temperature and concentration from certain stationary values h0 ; T 0 ; C 0 ; r0 is the melt density at T ¼ T 0 ; C ¼ C 0 ; and ~ gthe gravity acceleration. For simplicity’s sake h; T; C will subsequently be referred to as enthalpy, temperature and concentration. Eqs. (1)–(6) are supplemented by initial and boundary conditions. Let us break down the region over which the solution is to be computed into meshes and average the enthalpy values in each of the mesh cells. Let us express the average mesh temperature T in terms of enthalpy h: 8 ðh h1 Þ > > Tn þ for hph1 ; > > > C ps ðhÞ < T ¼ Tn for h1 ohoh1 þ x; > > > h h x 1 > > for hXh1 þ x; : Tn þ C pl ðhÞ where C ps ðhÞ is heat capacity of the solid phase of a substance; C pl ðhÞ is the heat capacity for the liquid phase; h1 is the maximum value of enthalpy for the solid phase, x is phase transition latent heat, T n is crystallization temperature. The state of a substance into mesh cells with enthalpy values h1 ohoh1 þ x is intermediate (transient) from thermodynamics viewpoint between solid and liquid. For pure substances T ¼ T 0 ¼ const. For substances containing an impurity T n is dependent on C: T n ¼ T n0 þ mC;
~ ¼ 0; rV
ð3Þ
@C ~ þ V rC ¼ rD rC: @t
ð4Þ
where m is the slope of the solidus line in the phase diagram. Let us use L and S to denote the fraction of liquid and solid phase respectively in the inter-
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
mediate-state region h1 ohoh1 þ x: 8 0 for hph1 ; > >
333
and porosity coefficient K, and is opposed to the ~ in direction: vector V ~ @V ~rÞ V ~ ¼ 1 rp þ n DV ~ þ ðV @t r0 ~ nV K; ð7Þ S where n ¼ rm ; S is square of mesh section normally 0 ~, a is an empirical ‘‘clusterizato velocity vector V tion constant’’ which together with K describes the double-phase medium nature. In physical sense a will depend firstly on the shape and averaged size of clusters, i.e. on the surface/volume ratio of solid particles putting up resistance to the flow motion. When K ¼ 1 (in solid state) Eq. (7) changes into ~ ¼ 0, and when K ¼ 0 it goes over to usual V hydrodynamics Eq. (2). Impurity transfer in meshes where the substance is in transient state will be described by the following Equation [18]: ga ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ~
Let us assume that the quantities r; l; C p ; D are constants for the solid and liquid phases of the substance. Meshes with transient enthalpy values will also have intermediate values of r; l; C p ; D which we will assume to be linear with h. If we use X s to denote a parameter value for the solid phase and X L for the liquid phase, the value of X for the intermediate state can be obtained from the formula: X ¼ X S L þ X L S ; T¼
ðh L xÞ ; Cp
where the parameter X is r; l; C p ; D. Let us assume that the intermediate melt region is filled with clusters, i.e. population of atoms (molecules) associated in solid structures. The size and number of clusters can be arbitrary. The clusters’ heat movement velocity is significantly lower than that for single melt atoms because they have much bigger sizes. The results of moleculardynamic calculations confirm this assumption. For example, the estimation of cluster diffusion coefficient with the number of associated atoms of about 70 in GaP melt given in [2] shows the value D ¼ 8:6 107 cm2 =s, whereas for single phosphorus atoms D ¼ 8:1 105 cm2 =s. The difference is seen consists 2 orders of magnitude. The given value allows the multiatomic cluster movement velocity to be neglected within the model. Then the melt flow in the medium that is filled not tightly with motionless clusters can be presented as a liquid flow in a porous double-phase body [17]. The melt flow in meshes where the substance is in intermediate state will be described using the porous solid approximation with the double-phase (porosity) coefficient K ¼ LS . To be able to make such a description let us make Eq. (2) include a resisting force due to conditions of porosity. This force varies as the melt viscosity, as flow velocity
@C ~ @S þ V rC ¼ rD rC þ ð1 K p ÞC; ð8Þ @t @t where K p o1 is impurity segregation coefficient. Finally, the mathematical model of heat-mass transfer for the whole area, where the crystallization process takes place, includes Eqs. (1), (3), (7), (8) in which the temperature T and values of r; l; C p ; D; L ; S are single-valued and continuous one-variable functions of the enthalpy h. The hallmark of the proposed model is that there is no localized phase transition interface. This boundary is fluctuated (‘‘smudged’’) forming a region of meshes whose enthalpy values refer to an intermediate state of the substance. Nevertheless the crystallization front velocity can be defined as the rate at which the proportion of the solid/liquid phases of the substance changes. This approach allows one to describe impurity segregation in the process of crystallization together with the impurity diffusion–convection transfer in the melt. As far as the hydrodynamics equations are concerned the adhesion condition for the flow velocity is laid down at the boundary between the solid phase and the intermediate state of the substance being crystallized. And the additional resistance to the flow in the intermediate state region is described in terms of porous body
L
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
334
approximation in which the porosity coefficient is determined by the solid-to-liquid phase ratio.
3. Computation algorithm Let us introduce a time grid with variable pitch t and use an implicit discretization scheme. This yields the following equation system: ~ rTðhÞ þ h ¼ h ; raðhÞ rTðhÞ þ V Cpt Cpt ~ þ ðV ~ rÞ V ~þ n DV ¼
ð9Þ
~ ~ V nV þa K t S
~ V 1 rp ðbT ðT T 0 Þ þ bC ðC C 0 ÞÞ~ gþ ; t r0 ð10Þ
~ ¼ 0; rV
ð11Þ
S ð1 K p ÞC þ L C t S ð1 K p ÞC þ L C ; ð12Þ ¼ t ~ where the symbols V ; S ; C refer to the values of those quantities at the preceding point of time. The system of equations (9)–(12) is linearized and solved by the method proposed in [19] and described in [20]. This method consists from the Patankar method [21], exponential transformation method and the conjugate gradient method with preconditioning by the incomplete factorisation method [22]. ~rC þ rDrC þ V
melting point). Meanwhile, the given model does not consider true (molecular) viscosity of the melt (n), but some effective (integral) viscosity (n ) of substance having been in transient state between fluid and solid. Let us call this model of integral viscosity increase one (IVI-model). The cause of integral viscosity increase near the interface, as it was mentioned above, could be the formation of high-drag low-mobile multiatomic clusters resisting the melt flow within the transient area (just is cluster formation or CF-model). In the present work the last term in the Eq. (7) describes the resistive force like as Darcy law for the fluids current in a porous medium. In CF-model the melt viscosity is assumed to be equal to its true value, and is either constant or slightly depending function of temperature. For establishing the quantitative relationship between IVI- and CF-models it is necessary to evaluate the empirical parameter a in the last term of the Eq. (7). For this purpose we have used for a long-time-knowing fundamental effect of compositional melt stratification in hermetic ampoule under the action of weak diametrical gravitational force (see theory e.g. in [24]). The action of stationary axis-perpendicular microacceleration leads to formation of monotonously increasing diametrical concentration profile of dopant, having K p o1, in the growing crystal along the gravity vector direction. Such the effect was observed repeatedly by various authors in the course of space-grown crystal study (see, e.g. [25,26]). As the test basis we have chose the paper of Duffar et al. [26], where GaSb:In crystallization experiments were performed by Bridgman technique both aboard the EURECA orbiter and on-ground. For sketch depicted in Fig. 1 the parametric
4. Empirical parameter evaluation In paper [23] the mathematical model is proposed allowing describe adequately the experimentally revealed anomalous distribution of impurity in low-doped semiconductor crystals grown under conditions of Photon series satellites orbital flight. In the foundation of this model the insight lays on the sharp increasing of melt viscosity in growth interface vicinity within the transient layer limiting by [T þ 10K; T ] isotherms (where T is
Fig. 1. The calculation scheme (indium dopant concentrations are given for initial GaSb charge parts before the melting). Growth rate—0:1 mm=s, initial axial thermal gradient—20 K/ cm, final axial thermal gradient-30 K/cm.
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339 3.1 3.0 2.9
∆C,× 1017at/cc
2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0 10.0
(1/v)⋅(dv*/dT), K-1
Fig. 2. Relationship between diametrical inhomogeneity ratio of In-dopant distribution in GaSb crystal, and integral viscosity parameter in transient boundary layer. Calculation is done for axis-perpendicular stationary microgravity level of 104 g0 , after 106 s since solidification beginning.
335
according both IVI- and CF-models. In the former the increase of integral viscosity was varied up to 100 times within transient area, and in the latter one the empirical clusterization constant a was modified. In both cases the resulting concentration difference (DC) of In impurity was calculated between the edges of crystal diameter. Fig. 2 shows the DC dependence on relative thermal gradient of integral viscosity ð1=nÞ ðdn =dTÞ, obtained by IVI-model. Fig. 3 shows the similar dependence of a, for which both models (IVI- and CF-) give identical diametrical non-uniformity DC in impurity distribution. The good agreement between calculated and experimental results was achieved in the work [23] when ð1=nÞ ðdn =dTÞ ¼ 6–7. These values correspond in Fig. 3 to a ð6:2–7:0Þ 1017 estimation. Assuming as the first approximation that quite conservative CF-model is valid for the most of semiconductors, in the current work we used for orbital crystallization process simulation the a ¼ 7 1017 value.
8.0
5. Model testing
α × 1017
6.0
In order to test the computer code, the numerical simulation was performed of the real crystallization experiments [26] within the frame of proposed CF-model both in space and on-ground. The calculation scheme is given in Fig. 1. Calculation was carried out in 2-D (X,Y) geometry. The values of physical parameters used in calculations are given in Table 1.
4.0
2.0
Table 1 Physical parameters of GaSb used in calculations
0.0 2.0
4.0
6.0 (1/v)⋅(dv*/dT),
8.0
10.0
K-1
Fig. 3. Relationship between the clusterization constant a and integral viscosity parameter in transient boundary layer. Calculation is done for the same conditions as that of Fig. 2.
calculations were made during 106 s at residual gravity vector direction of 104 g0 perpendicularly to growth axis. Calculations were performed
Melting point, K Density of the melt, kg/m3 Density of the solid, kg/m3 Thermal conductivity, Wt/(m K) Thermal capacity, J/(kg K) Kinematic viscosity, m2/s Diffusion coefficient of In in GaSb melt, m2/s Segregation coefficient of In in GaSb Temperature-expansion coefficient, K 1 Prandtl number
T rm rs lm Cp n D
985 6010 5600 8 175 0:38 106 1:2 108
Kp bT Pr
0.15 2:0 104 0.05
ARTICLE IN PRESS 336
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
Fig. 4 shows the calculated and experimental distributions of In isovalent dopant concentration along the reference crystal axis for the conditions of earth gravity. The calculation is made with a ¼ 0 in the Eq. (7). The agreement between the calculated results and experimental data is rather satisfactory. It obviously points to that the melt layer containing clusters has very small width and its influence on the convective heat and masstransfer in earth conditions can be neglected. Quite another matter is crystal grown in microgravity. The good agreement of calculated
and experimental axial distributions of In was achieved in Fig. 5 under zero gravity conditions only with the a ¼ 7 1017 value application. That indicates a really low level of convection and the existence of highly pronounced transient boundary layer during the experiments aboard the EURECA orbiter. In paper [26] the impurity radial distributions had different slopes for different cross-sections.
Fig. 4. Longitudinal In concentration distribution in GaSb for the conditions of earth gravity: 1—measured [26], 2—calculated using the clusterization parameter ratio a ¼ 0.
Fig. 5. Longitudinal In concentration distribution in GaSb for the conditions of zero gravity: 1—measured [26], 2—calculated using the clusterization parameter ratio a ¼ 7 1017 .
Fig. 6. Distribution of In impurity concentration in C 1019 at=cm3 (a) and equal In concentration lines in the longitudinal crystal section (b).
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
337
Fig. 7. Spatial distribution of longitudinal and transverse components of melt flow rate vector and melt flow rate fields. (a)—at the moment of 5 104 s, (b)—at the moment of 105 s.
The authors of [26] suppose that it indicates the change in the direction of residual gravity vector. In order to check this assumption the calculation was made of GaSb crystallization at g ¼ 104 g0 when the residual gravity vector assumed been perpendicular to the growth axis during the first 5 104 s of crystallization, and during the next 5104 s it had the opposite direction. Fig. 6 shows the spatial distributions of In concentration in the solid GaSb being crystallized, as well as concentration isolines obtained in result of this calcula-
tion according the CF-model. The presented data show that isolines firstly are curving towards one side, then are rectified, and finally are curving towards the other side. Thus confirming the model performance at the given assumption on the existence of residual gravitation and precession of its vector during solidification. Fig. 7 shows spatial distributions of flow rate vector components for axes X and Y at the time moments of 5 104 and 105 s in order to illustrate the flow nature in the considered case. These data
ARTICLE IN PRESS 338
V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
testify that the vortex reverses its direction when the direction of residual gravitation vector changes. The maximum value of flow rate module in the first case was equal to 0:34 105 m=s, and in the second case 0:39 105 m=s.
6. Conclusion The model presented in the paper allows the description of melt crystallization process in view of moving interface, convection heat and mass transfer and impurity segregation. For the first time the hydrodynamics model considers the structural pre-crystallization state of the boundary layer in the interface area. The transient region is firstly described as a double-phase medium containing three-dimensional solid phase clusters with a melt together. The presented model does not simulate the quantity and the size of the clusters near the growth interface, and takes into account only the influence of cluster-containing transient area on hydrodynamics of melt and distribution of impurity. The objective of the current work was the estimation of this influence. The calculations show that this influence in earth conditions is insignificant and it can be only detected when the external forces are small, e.g. in condition of microgravity. The specific feature of the model consists in the description of clusterization processes and those of substance crystallization in enthalpy (thermodynamic) variables, thus giving further the possibility to cross-over the results with the data obtained with other independent approaches, i.e. physical chemistry, molecular dynamics, etc. So there is a potential for further development and improvement of model representations. The next step proposes a quantitative description of clusters distribution by their sizes using molecular-dynamic approach and establishing of these data connection with local values of enthalpy. In the suggested version of model the semiempirical consideration was only given for the atmost case of melts structure ordering, i.e. for the formation of solid-state associated particles. Nevertheless the model gives quite adequate description of crystallization system response to
the change of key external conditions—temperature, crystallization rate, magnitude and direction of microacceleration vector, etc. The results of Indoped GaSb crystallization test calculations show a good agreement with the space experimental data. The most pronounced effect from the use of model can be expected for description of crystallization processes performed in the conditions of reduced gravity. That is the case when the boundary segregation phenomena near the crystallization front the most prominently affects the quality of the crystal growing.
Acknowledgements We would like to acknowledge support from the Russian Foundation for Basic Research, grant No. 03-02-16282.
References [1] A.Ya. Gubenko, High purity substances 4 (1994) 110 (in Russian). [2] S.V. Kotov, A.R. Lyutikov, Yu.P. Khukhryansky, et al., Tech. Phys. Lett. 28 (2002) 584. [3] J.R. Chelikowsky, J. Derby, V. Godlewsky, et al., J. Phys.: Condens. Matter 13 (2001) R817. [4] A.G. Ostrogorsky, P.S. Dutta, J. Crystal Growth 197 (1999) 1. [5] A.Ya. Gubenko, Microelectronic Eng. 45 (1999) 161. [6] A.Ya. Gubenko, Crystallogr. Rep. 46 (2001) 81. [7] H. Sasaki, E. Tokizaki, K. Terashima, S. Kimura, Japan J. Appl. Phys. 33 (1994) 3803. [8] S. Kawanishi, H. Sasaki, K. Terashima, S. Kimura, Japan J. Appl. Phys. 34 (1995) L1509. [9] S. Kimura, K. Terashima, J. Crystal Growth 180 (1997) 323. [10] M. Langen, T. Hibiya, M. Eguchi, I. Egry, J. Crystal Growth 186 (1998) 550. [11] W.-K. Rhim, K. Ohsaka, J. Crystal Growth 208 (2000) 313. [12] A.V. Kartavykh, E.S. Kopeliovich, M.G. Milvidskii, V.V. Rakov, Crystallogr. Rep. 43 (1998) 1075. [13] A.V. Kartavykh, E.S. Kopeliovich, M.G. Mil’vidskii, V.V. Rakov, Microgravity Sci. Technol. XII (1999) 16. [14] A.V. Kartavykh, E.S. Kopeliovich, M.G. Mil’vidskii, V.V. Rakov, J. Crystal Growth 205 (1999) 497. [15] R.W. Smith, X. Zhu, M.C. Tunnicliffe, et al., Proceedings of the International Conference ‘‘Spacebound-2000’’.
ARTICLE IN PRESS V. Ginkin et al. / Journal of Crystal Growth 270 (2004) 329–339
[16]
[17]
[18] [19]
Banff, Canada, 2000. p. 96 (also www.space.gc.ca/asc/pdf/ 2000 SMITH1.PDF). R.W. Smith, X. Zhu, M. C. Tunnicliffe, et al., Proceedings of the First International Symposium on Microgravity Research and Application in Physical Sciences and Biotechnology, Sorrento, Italy, Vol. II. 2000. ESA publ. SP-454, p. 887. V. Ginkin, O. Naumenko, M. Zabudko, et al., in: Proceedings of the International Symposium on Advances in Computational Heat Transfer, (CHT-04). Midnatsol Cruise, Norway, 2004, p. 226. V. Timchenko, P.Y.P. Chen, D.G. de Vahl, et al., Int. J. Heat and Mass Transfer 43 (2000) 963. V.P. Ginkin, in: Proceedings of the Second International Symposium on Advances in Computational Heat Transfer, Palm Cove, Australia, 2001, p. 1161.
339
[20] V.P. Ginkin, S.M. Ganina, Math. Model. 15 (2003) 53 (in Russian). [21] S.V. Patankar, D.B. Spalding, J. Heat Mass Transfer 15 (1972) 1787. [22] V.P. Ginkin, A.V. Kulik, O.M. Naumenko, in: Proceedings of the 4th International Conference on Supercomputing in Nuclear Applications (SNA 2000), Tokyo, Japan, 2000, p. 78. [23] V.K. Artemyev, V.I. Folomeev, V.P. Ginkin, et al., J. Crystal Growth 223 (2001) 29. [24] S.A. Nikitin, V.I. Polezhaev, A.I. Fedyushkin, J. Crystal Growth 52 (1981) 471. [25] V.S. Zemskov, V.N. Kubasov, I.N. Belokurova, et al., Apollo-Soyuz Test Project. Summary Science Report. V.1. Washington, NASA edit., 1977, p. 539–553. [26] T. Duffar, M.D. Serrano, C.D. Moore, et al., J. Crystal Growth 192 (1998) 63.