Modeling interfacial attachment kinetics during crystal growth from the liquid phase
S. Brandon1, I.G. Rasin1, A. Virozub1, O. Weinstein2
1-Department of Chemical Engineering, Technion, Haifa 32000, Israel
2-Institut für Kristallzüchtung, Max-Born-Str.2, D-12489 Berlin, Germany
We discuss our recent advances in modeling interfacial attachment kinetics during crystal growth from the liquid phase. Two types of systems, in which both transport phenomena and interfacial attachment kinetics are of importance, are addressed. The main focus is on large-scale crystals directionally solidified from their melt; the behavior of such systems is typically dominated by transport phenomena. However, in our case we discuss a class of melt-growth systems in which interface attachment kinetics impact the growth process and quality of the final product. Another type of system additionally (albeit briefly) discussed, concerned with the growth of crystals from supersaturated solutions, is more strongly affected by interfacial attachment kinetics. In both cases modeling efforts have traditionally focused either on large scale analyses of transport phenomena in the crystal growth system or on limited and/or localized theoretical treatments of liquid/crystal interfacial (e.g. stability) phenomena.
The approach presented here is based on rigorous modeling of interfacial attachment kinetics at the liquid/crystal interface, coupled with the analysis of transport phenomena at different length scales depending on the particular problem under consideration. In the case of melt growth the appearance of facets on the liquid/crystal interface is demonstrated, and the sensitivity of this phenomenon to important kinetic parameters is discussed. When analyzing crystal growth from solution, the impact of flow in the liquid phase on interface morphology is addressed. In particular the possibility of designing flow patterns capable of suppressing unwanted morphology is discussed.
1. Introduction
Crystal growth from the liquid phase is an important process both in nature and in synthetic environments. Systems of relevance include the synthesis of pharmaceuticals, the formation of ice, the growth of large-scale semiconductor and dielectric crystals, and the precipitation of salts in an aqueous environment. The physics and chemistry underlying these processes include transport phenomena as well as interfacial attachment kinetics. Modeling of these processes is an important tool used for their analysis, design and control (see e.g. [1]).
Traditionally, large scale modeling of crystal growth from the liquid phase has focused on transport phenomena in the growth system while ignoring or, at best, simplifying kinetics of attachment at the interface between the liquid and solid phases. This is especially true with respect to crystal growth from the melt which is indeed, in many cases, dominated by transport phenomena. However, interfacial attachment kinetics does become an important factor in the case of certain melt-grown materials such as oxides [2] or semiconductor systems grown along singular crystallographic directions (e.g. <111> in silicon [3]). In such cases, the liquid/solid interface typically deviates from the melting point isotherm position while exposing smooth facets supporting a spatially dependent undercooling. Facets along the melt/crystal interface have been associated with inhomogeneous dopant segregation due to the so-called "facet effect" [4], striations due to the step-growth mechanism governing growth on vicinal surfaces [5], twinning related to 2-D nucleation phenomena [6-8], and surface instability due to sudden switching between kinetic mechanisms [3,9].
While it is recognized that solution growth systems are often strongly affected by the combination of interfacial attachment kinetics and transport phenomena [10,11], modeling of these processes in relevant large-scale systems is even less common than in the case of melt-growth processes. The current state of the art includes a few large-scale 3-D models of transport phenomena [12-14] as well as a number of 1-D microscopic models of step-dynamics [15,16,10] applicable to both large and small scale systems. To the best of our knowledge, none of these test methods for directly suppressing step-bunching; this is an unwanted phenomenon in which unit steps on a vicinal surface move one towards another in a manner forming bunches that may eventually lead to the formation of defects such as kinetic striations and liquid inclusions (see. e.g. [17]).
In Section 2 we outline the theory of interface attachment kinetics and transport phenomena. Following this, in Section 3, we briefly review our computational approach. Results, provided in Section 4, summarize two studies of combined transport and interfacial attachment kinetics. First we analyze a specific melt growth system in which we examine the sensitivity of interfacial geometry and temperature distribution to different kinetic growth scenarios and parameters. Next, we briefly discuss simulations of flow modulation as a possible approach to the control of step bunching phenomena during crystal growth from solution. Finally, we end this manuscript with an appropriate conclusions section.
2. Theory
In the following we review the theory underlying both interface-attachment kinetics and heat transport in melt-growth systems. Relevant changes necessary when applying the same approach to solution-growth systems will be briefly given towards the end of this section.
2.1 Interface attachment kinetics
Crystal growth is driven by a difference in chemical potential between the solid and liquid phases, favoring the formation of a crystalline phase. In single component systems, exhibiting a planar interface between the solid and liquid phases, this driving force can be shown to be proportional to the difference between the material's melting point and the actual interfacial temperature (DT º Tmp - T). The growth rate normal to the interface (Vn) is related to DT, where in the most general form it is convenient to use the following relation between growth rate and undercooling:
. (1)
In Eq. (1) the kinetic coefficient (b) may depend both on the undercooling and on the crystallographic orientation of the surface, where q, representing this orientation, is usually defined as the deviation between the surface's orientation and that of a given singular orientation.
The detailed form of b(DT,q) depends primarily on the type of surface considered. Singular surfaces, which are below their thermodynamic roughening temperature, are (theoretically) atomically smooth and can grow only via step generation, mostly due to two dimensional nucleation or the rotation of screw dislocations. Vicinal surfaces advance by the lateral motion of steps, and atomically rough surfaces grow by the direct incorporation of growth-units from the liquid phase. Each of these mechanisms is characterized by a specific functional dependence of the kinetic coefficient on q and DT, given by:
, (2)
where br, bst, C, B, A can be treated as constants (note that bst is the kinetic coefficient associated with motion of steps on a vicinal surface).
Fig. 1. Kinetic coefficient of silicon versus misorientation (from the <111> direction) and undercooling. Both screw-dislocation and 2-D nucleation singular surface growth mechanisms are considered
The consolidation of these four equations into one overall expression for the kinetic coefficient is non-trivial (see [18]). The resultant plot for the kinetic coefficient versus undercooling and the deviation from a specific singular orientation, for the case of silicon, is presented in Fig. 1. Notice the fact that a constant coefficient, characteristic of a rough surface, is observed both for large undercooling values (kinetic roughening) and for large misorientation values (thermodynamic roughening).
2.2 Heat transport
As mentioned above, melt growth is governed by a combination of interfacial attachment kinetics and heat transport. The equations describing heat transport in the melt, crystal and crucible (when relevant) are given by:
, (3)
where i = m,s,c stands for melt, crystalline and crucible phases respectively, while k,r,Cp are the relevant thermal conductivity, density and heat capacity, respectively. The term involving the velocity vector (v) is active only within the molten phase. In addition, note that internal radiative transport, often important within the crystal (and sometimes melt) is not considered in this manuscript. Studying systems which participate (i.e. absorb, emit and possibly scatter thermal radiation) is beyond the scope of the analysis presented here (see e.g. [19] for the relevant state-of-the-art).
When i = m in Eq. (3) we must account for melt convection either by solving an appropriate equation of change of momentum (e.g. the Navier-Stokes equation) together with the continuity equation, or by inserting a pre-defined velocity profile. In all melt-growth results presented here we neglect effects of convection (i.e. v = 0), while in the case of solution growth we assume a pre-defined velocity profile as described in Section 4.2 below.
Solving Eq. (3) requires that we pose an initial condition for the temperature field as well as appropriate boundary conditions along all internal and external boundaries. In the event that the crucible is accounted for in the model, continuity of heat flux as well as continuity of temperature must be applied along all crucible/melt and crucible/crystal internal boundaries. Boundary conditions along external surfaces depend on the system studied and on available data. As discussed in Section 4.1 below, in the melt-growth system considered here it is assumed that the temperature profile along all external boundaries is known; note that in the case described in Section 4.1 the crucible is not modeled. In many other cases (see e.g. [20]) external boundaries are assume to interact via natural convection and enclosure radiation with a known external environment (i.e. the furnace).
Along the melt/crystal interface two conditions must be met. The first, accounting for a balance between heat fluxes normal to the interface and latent heat of solidification released during crystal growth, is given by:
, (4)
where in all cases presented here the heat flux vectors at the melt/crystal interface in the melt (qm) and crystal (qs) phases are simply given by Fourier's law of heat conduction. In addition, in Eq. (4), nms and DHv are, respectively, the unit vector normal to the melt/crystal interface and the volumetric latent heat of solidification. The second condition to be met can be viewed as a modified isotherm condition according to which the relation between the temperature and normal growth rate distributions along the melt/crystal interface must be consistent with the relevant interface attachment kinetics according to Eqs. (1,2).
2.3 Crystal growth from solution
In the case of solution growth the driving force for crystal growth is the degree of supersaturation of the solution, which is directly related to the difference in chemical potential (of the nutrient) between the liquid and solid phases. In the case of a dilute system the supersaturation is given by:
, (5)
where c and ceq are, respectively, the concentration of solute and its equilibrium concentration (at the given system temperature).
In principle the kinetic relations given in Eqs. (1,2) apply in this case with the supersaturation (s) replacing the undercooling (DT) and, obviously, with appropriate changes in dimensions and values of kinetic parameters (A,B,C,br,bst). In this manuscript (i.e. in Section 4.2) we consider a vicinal surface, growing from solution, along which we resolve unit steps whose lateral growth velocity is given by:
(6)
This equation for step velocity is consistent with the expression for normal growth velocity of vicinal surfaces given in Eq, (2) (with s replacing DT).
The equation of change of solute concentration, which we must solve in the case of solution growth, is identical to Eq. (3) without the riCpi pre-multiplier on the left hand side (i.e. riCpi is replaced by unity), with c replacing T, and with the diffusivity of solute (D) replacing the thermal conductivity ki. The interfacial condition (Eq. 4) is replaced by a similar relation accounting for mass balance of solute across the interface (i.e. the normal growth rate is directly proportional to the mass flux of solute from the solution, given by Fick's law while accounting for the convective as well as diffusive flux terms). In addition, in analogy to the case of melt growth, the normal growth rate must be consistent with the local supersaturation and with the relevant kinetic relation. All other boundary conditions as well as the initial condition, which are system specific, are given in Section 4.2 below.
3. Computational approach
The computational approach we use to analyze melt growth is based on the Galerkin Finite Element Method (GFEM) together with an interface motion operator method developed specifically for tracking the melt/crystal interface according to kinetics described by Eqs. (1,2) and depicted in Fig. 1. A full description of the numerical approach is available in [18,20,21]. Here we provide a brief outline of its essentials.
The temperature field in the crystal growth system is expanded in a set of basis functions consistent with mesh used to span the entire computational domain, while the melt/crystal interface position is expanded in a set of basis functions consistent with the mesh used to span the interface. In the sample results given in Section 4.1 we discuss a three-dimensional system in which we use an unstructured mesh involving tetrahedral quadratic elements, i.e. 3-D and 2-D quadratic basis functions are respectively used to approximate the temperature field and melt/crystal interface position.
The approximated temperature field and interface position are inserted into Eq. (3). The time derivative of temperature in this equation is formulated while accounting for mesh velocity [22], and associated derivatives of nodal position and temperature values are discretized using a backward Euler formulation. The resultant discretized version of Eq. (3) is weighted and integrated over the entire computational domain yielding nonlinear algebraic residual equations that are evaluated at each time-step using a the Newton Raphson method. Interface motion is activated based on interfacial undercooling and orientation via Eqs. (1,2) using an interface motion operator approach described in [20]. An iterative procedure ensures a consistent solution of Eqs. (1,2,3,4) and associated boundary conditions at each time step.
The approach we take to solve solution crystal growth problems depends on the length-scale of the system being considered. In the case of the 2-D combined mass transport and vicinal surface growth problem discussed in Section 4.2., we use a phase-field based approach for the resolution of the stepped structure of the surface (see. e.g. [23]). The resultant equation is solved using a standard finite difference method with a uniform 1-D mesh conforming with the lower part of a 2-D mesh applied in the (rectangular-shaped) bulk liquid phase. The Lax-Wendrof method (see .e.g. [24]) is employed in the simulation of the advection-diffusion (Eq. (3) with modifications described in Sections 2.3, 4.2) process in the bulk using an efficient graded mesh.