PHYSICAL REVIEW E 89, 012405 (2014)
Phase-field-crystal simulation of nonequilibrium crystal growth Sai Tang,1,* Yan-Mei Yu,2,† Jincheng Wang,1,‡ Junjie Li,1 Zhijun Wang,1 Yaolin Guo,1 and Yaohe Zhou1 1
State Key Laboratory of Solidification Processing, Northwestern Polytechnical University, Youyi Western Road 127, 710072 Xi’an, China 2 Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China (Received 2 September 2013; revised manuscript received 29 November 2013; published 13 January 2014) By using the phase field crystal model, we simulate the morphological transition of the crystal growth of equilibrium crystal shape, dendrite, and spherical crystal shape. The relationship among growth morphology, velocity, and density distribution is investigated. The competition between interface energy anisotropy and interface kinetic anisotropy gives rise to the pattern selection of dendritic growth in the diffusion controlled regime under low-crystal-growth velocities. The trapping effect in density diffusion suppresses morphological instabilities under high-crystal-growth velocities, resulting in isotropic growth of spherical crystal. Finally, a morphological phase diagram of crystal growth is constructed as function of the phase field crystal model parameters. DOI: 10.1103/PhysRevE.89.012405
PACS number(s): 81.10.Aj, 64.70.dm, 61.50.Ah, 89.75.Kd
I. INTRODUCTION
Pattern formation in crystal growth is a fascinating phenomenon of longstanding interest for scientists [1]. A vast variety of amazing morphologies can be produced from disordered precursors. For example, it is found that crystal growth morphology changes from faceted structure at low undercooling to dendrite at intermediate undercooling, and to spherical crystal at super high undercooling in solidification of silicon [2]. The crystal growth of organic material also exhibits the similar morphological transition of petal shape, dendrite, and isotropic circular shape with supersaturation increasing [3]. The challenge is to understand the spontaneous formation and the growth mechanism of these diversified patterns. A cornerstone of related theoretical works is the recognition that growth patterns may be determined by the interplay between microscopic interfacial dynamics and external macroscopic forces. Ivantsov solution theory [4] gives the macroscopic description of the interface instability controlled by the mass or thermal diffusion, and on the other hand, the Mullins-Sekerka stability analysis [5] and the microscopic solvability theory [6–8] give the microscopic theory on the effect of the interface energy on the interface instability. The role of anisotropy during morphological selection of crystal growth is also insightfully discussed from thermodynamics and kinetics aspects [9–11]. In experiments, the Hele-Shaw cell provides a simple method for pattern formation studying [9,10]. In Hele-Shaw experiments, faceted, dendritic, and dense-branching morphologies are observed, and the morphological transformation driven by applied air pressure and anisotropy of surface tension has shed light on the understanding of the interplay of macroscopic thermodynamic factor and microscopic kinetics effect. For a deep understanding of morphological transformation mechanism of crystal growth, atomic-scale simulation will give a lot of useful information through tracking crystal growth processes directly. However, such work needs to model atomic processes in realistic time scale and therefore
*
[email protected] [email protected] ‡ [email protected] †
1539-3755/2014/89(1)/012405(6)
is very challenging. Recently, the phase field crystal (PFC) model has been widely used to study various interfacial and microstructure evolutions in crystal growth [12–18]. In the PFC model, the atomic scale structure spanning the diffusive time scale of crystal growth can be captured, so that crystal growth can be simulated by considering a microscopic factor such as crystalline anisotropy, and macroscopic growth driving force. Moreover, a large parameter space of crystalline anisotropy and growth driving force can be spanned easily in the PFC model, allowing us to reinvestigate the crystal growth ranging from near-equilibrium to far-from-equilibrium in the frame of the PFC model. In this study, we reproduce a rich variety of crystal growth morphology, including equilibrium crystal shape (ECS), dendrite, and single-crystal and polycrystal with spherical shape. The relationship between growth morphology, velocity, and density distribution is investigated. Three types dendritic patterns are obtained, including interface energy dendrite (IED) with the tip orientation being aligned with a laboratory frame of reference, tip splitting dendrite (TSD), and interface kinetic dendrite (IKD) with the tip orientation twisting 30◦ relative to the laboratory frame of reference, which shows the competition between interface energy anisotropy and interface kinetic anisotropy. Finally, a complete crystal growth morphology is constructed as a function of macroscopic driving force and microscopic crystalline anisotropy in the PFC model. II. PFC MODEL AND SIMULATION DETAILS
We use the two-dimensional (2D) PFC model developed by Elder et al. [12,13]. The PFC model can be derived from classic density functional theory of freezing 2 φ φ4 α + λ q02 + ∇ 2 φ + g , (1) F = d r 2 4 where φ is the time-averaged atom number density measured with respect to a reference liquid state, which is periodic in solid phase and homogeneous in liquid phase, q0 is the first-peak position of structure factor S(k) of liquid phase, being related to the lattice constant of crystal phase. Parameters α and λ deduced from the shape of polynomial approximation of S(k) determine the liquid-phase compressibility, bulk
012405-1
©2014 American Physical Society
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
modulus, and solid-liquid interface properties [19–21]. The dynamic equation is given by Swift-Hohenberg model of pattern formation, δF ∂φ = 2 , ∂τ δφ
crystal growth velocity conveniently. The position of each atom is identified as the local maximum of density field. The atom configuration of crystal growth is visualized in terms of the position of atom by using the visualization software [22].
(2)
where is related to the diffusion coefficient of atoms [13]. By defining x = rq0 , ψ = φ λqg 4 , = − λqα 4 , F˜ = λ2gq 6 F , and 0
0
and corresponding dimensionless dynamic equation reads
(4)
where corresponds to crystalline anisotropy, which is chosen to be 0.1 ∼ 0.5, indicating an increasing crystalline anisotropy [14]. The crystal growth begins with a circular nucleus located at the center of simulation box and surrounded by homogenous liquid. The initial density of the nucleus is given by one-mode approximation which reads √ √ ψ = A[cos(qx) cos(qy/ 3) − cos(2qy/ 3)/2] + ψ¯ 0 , (5) where ψ¯ 0 is the initial density of liquid phase, which determines supersaturation as σ = (ψ¯ 0 − ψ¯ Le )/(ψ¯ Se − ψ¯ Le ), where ψ¯ Se and ψ¯ Le are, respectively, the solid phase and liquid phase lines of the 2D hexagonal-liquid phase coexistence region in the 2D PFC phase diagram [13], and A is the amplitude of density waves, and q is the magnitude √ of principal reciprocal lattice vectors with the value of 3/2. Thus, the lattice constant of hexagonal structure is given by a = 2π/q. The radius of initial nucleus is chosen to be 5a empirically. The dynamic equation is solved through semi-implicit Fourier spectral method. The dimensionless grid size is dx = π/4, and time step is dt = 0.75. The finite size effect is an important factor that may affect the diffusion field of mass and energy in front of solid-liquid interface. Finite-size effect is expected to be negligible for the growth of ECS with very low velocity. For crystal growth with high velocities, finite-size effect is limited by implementing large simulation boxes. The rectangular grid with 2400 × 2400 is used for low-velocity growth and diffusionless growth, while the rectangular grid with 3600 × 3600 is used for dendritic growth. In all cases, the density far from solid-liquid interface front is maintained to be the initial ψ¯ 0 values. In order to obtain crystal growth velocity, we firstly filter the periodic oscillating density field through iterative adjacent averaging procedure, k+1 ψ¯ (i,j ) =
1
k ψ¯ (i±n,j ±m) /9,
A. Crystal growth morphology
0
dimensionless time t = λq06 τ , the dimensionless free energy functional reads as ψ ψ4 F˜ = d x [− + (∇ 2 + 1)2 ]ψ + , (3) 2 4 ∂ψ = ∇ 2 (δ F˜ /δψ), ∂t = ∇ 2 {[− + (∇ 2 + 1)2 ]ψ + ψ 3 },
III. RESULTS AND DISCUSSION
(6)
n,m=0
where k is the iteration time. The position of solid-liquid interface is defined by the isoline of smoothed density with the value of (ψ¯s + ψ¯ l )/2, where ψ¯s and ψ¯ l are the smoothed density of solid phase and liquid phase, respectively. By tracking the position of solid-liquid interface, we obtained
The crystal growth morphologies are summarized in Fig. 1. ¯ are chosen based on the PFC The simulation parameters ( ,ψ) phase diagram, as shown in Fig. 1(s). The increasing and ψ¯ correspond to the increasing crystalline anisotropy and the increasing growth driving force, respectively. For a given value, the closer to the solid-phase line the ψ¯ value is, the more rapid increase in the growth velocity is obtained usually. The crystal growth is simulated for three different values, = 0.15, 0.3, and 0.4. For each case, the crystal growth undergoes four growth stages with the ψ¯ increasing. In the first stage, when ψ¯ 0 is small, the crystal growth demonstrates ECS that changes from faceted hexagonal [Fig. 1(a) for = 0.4] to faceted hexagonal with round corners [Fig. 1(f) for = 0.3], and then to rough spherical morphology [Fig. 1(m) for = 0.15]. During the faceted ECS growth, we can observe that the solid-liquid interface grows in the layer-by-layer mode which exhibits the nucleation of a new layer and subsequent growth by absorbing atoms into kinks on edges. As ψ¯ 0 increases, crystal growth enters into the second stage. In this stage, crystal growth is dominated by dendritic morphology which demonstrates rich morphological variations caused by different ψ¯ 0 and . Shown in Fig. 2 is the convergent dendritic tip velocities of IED and IKD as the time goes by, indicating that the dendritic growth has entered into the steady state. After the growths of IED and IKD reach the steady state, the dendritic tips of IED and IKD keep the constant state while the dendritic arms grow larger and larger. For the case of strong crystalline anisotropy, such as = 0.4, the preferred growth direction of dendrite tip is dominated by interface energy anisotropy, where the principal arms of dendrites [Fig. 1(b) and (c)] are aligned with the laboratory frame as denoted by the dashed line. In terms of atom arrangement (inset A) we can find that the principle arms of dendrite are aligned with the close-packed direction. Such orientation selection is dominated by the interface energy anisotropy inhered in crystalline structure, which is consistent with the surface-tension dendritic growth as observed in Hele-Shaw experiments [9]. Here, we call this type of dendrite as interface energy dendrite (IED). When the interface energy anisotropy is of intermediate magnitude, such as = 0.3, the dendritic growth morphology changes from IED [Fig. 1(g)] to tip splitting dendrite (TSD) [Fig. 1(h) and (i)], and then to the dendrite with principal arms twisting 30◦ in relative to the laboratory frame [Fig. 1(j)]. The principle direction of dendrite tip is aligned with the open-packed direction, as seen from inset B. Such 30◦ angle twist has been observed in experiments of 2D hexagonal lattice crystal growth confined by two slices of parallel glasses, which is attributed to the competition of interface energy anisotropy and interface kinetic anisotropy [3]. As the experimental and theoretical results have proved, the interface
012405-2
PHASE-FIELD-CRYSTAL SIMULATION OF . . .
PHYSICAL REVIEW E 89, 012405 (2014)
FIG. 1. (Color online) The crystal growth morphologies simulated under different ψ¯ 0 and with increasing growth driving force and increasing crystalline anisotropy, respectively. Patterns in (a)–(p) are obtained through performing the simulations denoted by crosses in the PFC phase diagram (s). (q) and (r) are the patterns of (k) and (j) in reciprocal space, respectively. Inset figures A–E at the bottom are the density map of the parts enclosed by square box in (c), (j), (k), (l), and (m), respectively. The 2400 × 2400 simulation grid is used for the growth of ECS, CSC, and PSC, and 3600 × 3600 for dendrite patterns. The crystals shown in (a), (f), (m) contain around 20000 atoms, in (b), (c), (g), (h), (i), (j), (n) around 80000 atoms, and around 35000 atoms in (d), (e), (k), (l), (o), (p).
kinetics factor will become dominant for high crystal growth velocity [3,9,10]. Considering the twisting 30◦ dendrite tip orientation is caused by interface kinetics anisotropy, we call the dendrite as interface kinetic dendrite (IKD) that is similar to the surface-kinetic dendrite observed in Hele-Shaw experiments [9].
FIG. 2. (Color online) The dendritic tip velocity of IED in Fig. 1(g) and IKD in Fig. 1(j) growth versus time.
When the interface energy anisotropy is too weak to sustain the pattern of IED, such as = 0.15, the dendritic growth [Fig. 1(n)] is mainly of the IKD mode. This indicates that the orientation of dendritic growth can be controlled by interface kinetic anisotropy when the interface energy anisotropy is weak enough. The crystal growth enters into the third stage, as ψ¯ 0 increases further. During this stage, being independent of the values, the crystal growth retrieves to the isotropic growth. The grown crystal is of compact spherical shape that we call as compact sphere-shape crystal (CSC), as shown in Fig. 1(d) for = 0.4, Fig. 1(b) for = 0.3, and Fig. 1(o) for = 0.15. The ECS, dendrite, and CSC are of single crystal structure, which is proved by the regular hexagonal diffraction lattice in reciprocal space, as shown in Fig. 1(q). With ψ¯ 0 increasing further, crystal growth enters into the forth stage. As shown in Fig. 1(e) for = 0.4, Fig. (l) for = 0.3, and Fig. 1(p) for = 0.15, the obtained crystal morphology is of polycrystal structure, which is proved by the diffraction circle pattern in reciprocal space, as shown in Fig. 1(r). The polycrystal structure has lost crystalline anisotropy and therefore produces isotropic spherical morphology called as polycrystal sphere-shaped crystal (PSC) here.
012405-3
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
B. Crystal growth velocity
C. Density distribution
Through analyzing crystal growth velocity, we can find some common growth regimes behind different morphologies shown in Fig. 1. As shown in Fig. 3, the low-velocity growth regime A, the intermediate-velocity growth regime B, and the high-velocity growth regime C can be found for all cases of = 0.4,0.3, and 0.15. There also exists a transition zone T between B and C for the cases of = 0.3 and = 0.15. The T regime can be regarded as the signal of the competition of interface energy anisotropy and interface kinetics anisotropy. The relationship between growth morphology and velocity can be found, i.e., ECS grows in A regime, IED and TSD grow in B regime, IKD grows in T regime (partly in B regime), and CSC grows in C regime (PCS also grows in C regime, not shown here). The curves of V versus ψ¯ 0 shown in Fig. 3(a)–(c) are similar to the experimental curves of growth velocity versus increasing temperature undercoolings [23,24] or supersaturation [25]. It is easy to understand that increasing φ¯ means increasing growth driving force in PFC model that is equivalent to undercoolings or supersaturation. The curve of V versus σ is shown in Fig. 3(d). For the diffusion-controlled growth, there exists a relationship V = Aσ α , as predicted by the LKT model [26]. Here, we find that the crystal growths of ECS, IED, TSD, and part of IKD can be described by V = Aσ α with α = 4.45, nevertheless the crystal growths of other IKD and CSC deviate from such relationship. Therefore, we conclude that the ECS, IED and TSD growths are of the diffusion-controlled mode, whereas CSC growth is not diffusion-controlled, and IKD growth lies between the two growth modes.
Figure 4 illustrates the influence of crystal growth velocity on density distribution. The influence of crystal growth velocity on the density diffusion is explained firstly, which shows that the diffusion layer of the low-velocity growth [IED in Fig. 4(a)] is much thicker than that of the high growth velocity [CSC in Fig. 4(b)]. Next, the influence of growth velocity on interface structure is explained in terms of interface thickness that is defined as the distance from solid surface to liquid front. For the cases of = 0.4 and 0.3, there is a sudden increase of the interface thickness, as marked by the transition point M in Fig. 4(c). We can observe that the faceted solidliquid interface becomes diffusive after M. As we suppose, the diffusive interface leaves more possibilities for atoms to attach, which is advantageous to the large crystal growth velocity. Thus, the sudden increase of interface thickness illustrates the relationship between interface structure and crystal growth velocity. In Fig. 4(d) and (e), the density distribution is illustrated in terms of solid-phase density ψ¯ S , liquid-phase density ψ¯ L , and reduced partition coefficient k = ψ¯ S /ψ¯ L at solid-liquid interface with increasing V . For ECS and IED growths, ψ¯ S and ψ¯ L are in accordance with the values of the solid-phase and liquid-phase lines in PFC phase diagram, respectively, for example, ψ¯ S = −0.315 and ψ¯ L = −0.366 at point N as denoted by the black arrow in Fig. 4(d) and (e), which are close
FIG. 3. (Color online) The crystal growth velocity V versus initial liquid density ψ¯ 0 under three different values in (a)–(c). Shown in (d) is V versus the supersaturation σ , corresponding to (b), where the solid line is the fitting curves of V = Aσ α with α = 4.45.
FIG. 4. (Color online) The density distribution at the solid-liquid interface. (a) and (b) The smoothed density profile along the direction normal to solid-liquid interface, wherein the solid, dashed, dotted, and dot-dashed lines correspond to the advancing position, (c) the thickness of the solid-liquid interface versus the crystal growth velocity V , (d) the solid-phase and liquid-phase densities ψ¯ S (smoothed value) and ψ¯ L at the solid-liquid interface versus V , wherein ψ¯ 0 is the corresponding initial liquid density, (e) the reduced solid-liquid partition coefficient k versus V , as fitted by the Aziz’s model (red line) [27].
012405-4
PHASE-FIELD-CRYSTAL SIMULATION OF . . .
PHYSICAL REVIEW E 89, 012405 (2014)
to the equilibrium values, ψ¯ Se = −0.312 and ψ¯ Le = −0.367, as predicted by the PFC phase diagram [13]. Therefore, we suppose that the growths of ECS and IED are in the under equilibrium or near-equilibrium state. However, for TSD and IKD, ψ¯ S and ψ¯ L begin to deviate from the values of the solid-phase and liquid-phase lines, respectively. The increasing crystal growth velocity causes increasing k. At the point M when the interface thickness increases suddenly, ψ¯ S approaches the initial liquid phase density ψ¯ 0 [denoted by the blue arrow in Fig. 4(d)]. This is a signal of diffusionless growth mode, i.e., the solid-liquid interface advances so fast that density-field of the front of liquid phase has insufficient time to diffuse. The insufficient diffusion ahead of advancing solid-liquid interface causes a trapping phenomenon under high growth velocities. As shown in Fig. 4(e), the tendency of k → 1 arises up when growth velocity V increases up to the region of CSC, where ψ¯ S approaches ψ¯ L . The relationship of k vs. V is fitted in terms of k = (k0e + V /VD )/(1 + V /VD ), wherein equilibrium solid-liquid partition coefficient k0e = 0.8260 ± 0.0012 and atom diffusion velocity VD = 0.0573 ± 0.0020. The obtained relationship of k vs. V in the PFC simulation is consistent with the solute trapping phenomena in solidification of alloy, as described by the Aziz’s model [27]. We have noticed that the Aziz’s model may not work under very high velocities. In the Aziz’s model, the complete trapping limit given by k(V ) = 1 is approached asymptotically as V → ∞. In contrast, the work of Sobolev [28] predicts that there is an abrupt change of k at a finite velocity. In other words, complete trapping (k = 1) occurs at a well-defined velocity. The experimental evidence supports that the Sobolev description of solute trapping is accurate, for example, splat cooling experiments of Al-Mg performed by Galenko and Herlach [29] show a change from a eutectic to supersaturated solid solution at a finite velocity. Very recently, Humadi and coworkers investigated the solute trapping in PFC model [30]. Their results show that incorporation of a higher order time derivatives in the PFC model can lead to complete trapping occurring at a finite velocity, while the PFC model with only diffusive dynamics (i.e., first-order time derivative) can lead to the solute trapping as the Aziz’s model predicted. Since in this work, we use the PFC model with only diffusive dynamics, the consistence of our relationship of k vs. V with the Aziz’s model is reasonable. The similar trapping phenomenon has been found in a previous PFC simulation, which is resulted from the trapping of vacancy-diffusion in PFC solid [17]. Our result is in accordance with the previous PFC simulations [17,30]. The trapping effect in density diffusion exerts significant effect on crystal growth morphology in the PFC simulation. Generally, the morphological instability, such as branching, dendritic, and cellular interface, is mediated through mass diffusion field in terms of a ‘thumb’ rule of diffusion-controlled growth [31]. The trapping effect inhibits density redistribution and segregation, and hence the morphology instability is suppressed greatly. This contributes to the stable solid-liquid interface of CSC and PSC growth. Despite of the same spherical shape, the CSC and PSC are of different growth mechanisms with the circular ECS growth in Fig. 1(m). The former is attributed to the stability induced by trapping effect,
while the latter is due to the weak crystalline anisotropy determined by = 0.15. D. Morphological diagram
Finally, various crystal growth morphologies are summarized as function of σ and , as shown in Fig. 5. ECS appears in the low σ regime (for example, σ < 0.6 at the bottom of the axis). The anisotropic ECS with faceted interface structure and the isotropic ECS with rough interface structure are found in the upper panel with > 0.25 and the lower panel with < 0.25, respectively. When σ is of intermediate values (for example, 0.6 < σ < 2.8 at the bottom of the axis), the crystal growth mode turns to be dendritic growth, where the specific dendrite pattern is further determined by . In the upper panel of > 0.25, the dendritic growth is dominated by interface energy anisotropy and hence IED is obtained, while in the lower panel of < 0.25, the dendritic growth is dominated by interface kinetics anisotropy and hence IKD is obtained. The TSD region is between IED and IKD regions. As σ increases further (for example, σ > 2.8 at the bottom of the axis), crystal growth enters into diffusionless mode, where CSC and PSC are obtained. Figure 5 illustrates how the morphological selection of crystal growth is determined by growth driving force and crystalline anisotropy. The similar morphological diagram has also been obtained in the Hele-Shaw experiment for the patterns of vapor-liquid interface, which is composed by the faceted shape (region I), surface-tension dendrite (region II), tip-splitting dendrite (region III), and kinetic dendrite (IV) [9]. Our PFC morphological diagram elucidates more morphologies because the growth driving force and crystalline anisotropy span broader ranges. When the crystalline anisotropy (i.e., the y-axis variable) is weaker, as shown in the left-lower-half panel ( < 0.25, σ < 0.2) of Fig. 5, crystal growth directly enters into the IKD region after non-faceted ECS growth region with driving force increasing, and therefore no IED and TSD is obtained. Besides, the PFC simulations explore the morphological diagram when the growth driving forces
FIG. 5. (Color online) The phase-diagram of crystal growth morphology as function of PFC parameter and supersaturation σ that describe the crystalline anisotropy and the growth driving force, respectively.
012405-5
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
is very large, demonstrating CSC grows at high velocity and PSC at higher velocity. Larger parameter spaces in the PFC model allow us to obtain a full picture of morphological selection map of the crystal growth from near-equilibrium to far-from-equilibrium. Crystal growth under very strong anisotropy, = 0.75, is studied in previous PFC simulations, and shows faceted ECS and IED growth under diffusion-controlled mode and compact hexagon-shaped crystal growth under diffusionless mode. Here, a large simulation space with ranging from 0.1 to 0.5 is explored. More PFC growth morphologies are obtained in the range, such as TSD, IKD, and PSC, enabling us to summarize the dendritic growth at intermediate velocities and the spherical crystal growth at high velocities systematically.
CSC at high velocities, and PSC at extremely high velocities. Furthermore, the competition between interface energy anisotropy and the interface kinetics anisotropy causes the orientation of dendritic arms twisting from the close-packed direction to the open-packed direction, i.e., from IED to IKD. The TSD appears under the interplay of interface energy dominated growth and interface kinetics dominated growth. At high velocities, the trapping effect in density diffusion suppresses interface instability, which is responsible for the isotropic growth of CSC and PSC. Such a broad range of PFC crystal growth morphologies reflects the underling relationship between the morphology, velocity and density diffusion of the morphological transition during crystal growth. ACKNOWLEDGMENTS
In summary, our PFC simulations present a full picture of the velocity-dependent crystal growth morphology from near-equilibrium to far-equilibrium regime, including an ECS crystal at low velocities, the dendrite at intermediate velocities,
This work has been supported by Nature Science Foundation of China (Grant Nos. 10974228 and 51071128), National Basic Research Program of China (Grant No. 2011CB610401), and EU FP7, PHASEFIELD under no. 247504, Excellent Doctorate Foundation of Northwestern Polytechnical University.
[1] J. S. Langer, Rev. Mod. Phys. 52, 1 (1980). [2] R. P. Liu, T. Volkmann, and D. M. Herlach, Acta Mater. 49, 439 (2001). [3] P. Oswald and P. Pieranski, Smectic and Columnar Liquid Crystals (Taylor and Francis, New York, 2006). [4] G. P. Ivantsov, Dokl. Akad. Nauk. SSSR 58, 567 (1947). [5] W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 35, 444 (1964). [6] A. Barbieri and J. S. Langer, Phys. Rev. A 39, 5314 (1989). [7] D. Kessler, J. Koplik, and H. Levine, Adv. Phys. 37, 255 (1988). [8] M. Ben Amar and E. Brener, Phys. Rev. Lett. 71, 589 (1993). [9] E. Ben-Jacob and P. Garik, Nature (London) 343, 523 (1990). [10] E. Ben-Jacob, N. Goldenfeld, J. S. Langer, and G. Sch¨on, Phys. Rev. Lett. 51, 1930 (1983). [11] F. Liu and N. Goldenfeld, Phys. Rev. A 42, 895 (1990). [12] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys. Rev. Lett. 88, 245701 (2002). [13] K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004). [14] R. Backofen, A. Ratz, and A. Voigt, Philos. Mag. Lett. 87, 813 (2007). [15] S. Tang, R. Backofen, J. Wang, Y. Zhou, A. Voigt, and Y. M. Yu, J. Cryst. Growth 334, 146 (2011). [16] S. Tang, Z. J. Wang, Y. L. Guo, J. C. Wang, Y. M. Yu, and Y. H. Zhou, Acta Mater. 60, 5501 (2012).
[17] G. Tegze, G. I. T´oth, and L. Gr´an´asy, Phys. Rev. Lett. 106, 195502 (2011). [18] G. Tegze, L. Gr´an´asy, G. I. T´oth, J. F. Douglas, and T. Pusztai, Soft Matter 7, 1789 (2011). [19] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phys. Rev. B 75, 064107 (2007). [20] K.-A. Wu and A. Karma, Phys. Rev. B 76, 184107 (2007). [21] A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, Phys. Rev. E 80, 031602 (2009). [22] J. Li, Model. Simul. Mater. Sci. Eng. 11, 173 (2003). [23] R. Willnecker, D. M. Herlach, and B. Feuerbacher, Phys. Rev. Lett. 62, 2707 (1989). [24] C. Panofen and D. M. Herlach, Appl. Phys. Lett. 88, 171913 (2006). [25] F. C. Frank, Contemp. Phys. 23, 3 (1982). [26] Lipton, W. Kurz and R. Trivedi, Acta Metall. 35, 957 (1987). [27] M. J. Aziz, J. Appl. Phys. 53, 1158 (1982). [28] S. L. Sobolev, Phys. Rev. E 55, 6845 (1997). [29] P. K. Galenko and D. M. Herlach, Phys. Rev. Lett. 96, 150602 (2006). [30] H. Humadi, N. Ofori-Opoku, N. Provatas, and J. J. Hoyt, JOM 65, 1103 (2013). [31] R. Trivedi and W. Kurz, Acta Metall. 34, 1663 (1986).
IV. SUMMARY
012405-6
Phase-field-crystal simulation of nonequilibrium crystal growth Sai Tang,1,* Yan-Mei Yu,2,† Jincheng Wang,1,‡ Junjie Li,1 Zhijun Wang,1 Yaolin Guo,1 and Yaohe Zhou1 1
State Key Laboratory of Solidification Processing, Northwestern Polytechnical University, Youyi Western Road 127, 710072 Xi’an, China 2 Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China (Received 2 September 2013; revised manuscript received 29 November 2013; published 13 January 2014) By using the phase field crystal model, we simulate the morphological transition of the crystal growth of equilibrium crystal shape, dendrite, and spherical crystal shape. The relationship among growth morphology, velocity, and density distribution is investigated. The competition between interface energy anisotropy and interface kinetic anisotropy gives rise to the pattern selection of dendritic growth in the diffusion controlled regime under low-crystal-growth velocities. The trapping effect in density diffusion suppresses morphological instabilities under high-crystal-growth velocities, resulting in isotropic growth of spherical crystal. Finally, a morphological phase diagram of crystal growth is constructed as function of the phase field crystal model parameters. DOI: 10.1103/PhysRevE.89.012405
PACS number(s): 81.10.Aj, 64.70.dm, 61.50.Ah, 89.75.Kd
I. INTRODUCTION
Pattern formation in crystal growth is a fascinating phenomenon of longstanding interest for scientists [1]. A vast variety of amazing morphologies can be produced from disordered precursors. For example, it is found that crystal growth morphology changes from faceted structure at low undercooling to dendrite at intermediate undercooling, and to spherical crystal at super high undercooling in solidification of silicon [2]. The crystal growth of organic material also exhibits the similar morphological transition of petal shape, dendrite, and isotropic circular shape with supersaturation increasing [3]. The challenge is to understand the spontaneous formation and the growth mechanism of these diversified patterns. A cornerstone of related theoretical works is the recognition that growth patterns may be determined by the interplay between microscopic interfacial dynamics and external macroscopic forces. Ivantsov solution theory [4] gives the macroscopic description of the interface instability controlled by the mass or thermal diffusion, and on the other hand, the Mullins-Sekerka stability analysis [5] and the microscopic solvability theory [6–8] give the microscopic theory on the effect of the interface energy on the interface instability. The role of anisotropy during morphological selection of crystal growth is also insightfully discussed from thermodynamics and kinetics aspects [9–11]. In experiments, the Hele-Shaw cell provides a simple method for pattern formation studying [9,10]. In Hele-Shaw experiments, faceted, dendritic, and dense-branching morphologies are observed, and the morphological transformation driven by applied air pressure and anisotropy of surface tension has shed light on the understanding of the interplay of macroscopic thermodynamic factor and microscopic kinetics effect. For a deep understanding of morphological transformation mechanism of crystal growth, atomic-scale simulation will give a lot of useful information through tracking crystal growth processes directly. However, such work needs to model atomic processes in realistic time scale and therefore
*
[email protected] [email protected] ‡ [email protected] †
1539-3755/2014/89(1)/012405(6)
is very challenging. Recently, the phase field crystal (PFC) model has been widely used to study various interfacial and microstructure evolutions in crystal growth [12–18]. In the PFC model, the atomic scale structure spanning the diffusive time scale of crystal growth can be captured, so that crystal growth can be simulated by considering a microscopic factor such as crystalline anisotropy, and macroscopic growth driving force. Moreover, a large parameter space of crystalline anisotropy and growth driving force can be spanned easily in the PFC model, allowing us to reinvestigate the crystal growth ranging from near-equilibrium to far-from-equilibrium in the frame of the PFC model. In this study, we reproduce a rich variety of crystal growth morphology, including equilibrium crystal shape (ECS), dendrite, and single-crystal and polycrystal with spherical shape. The relationship between growth morphology, velocity, and density distribution is investigated. Three types dendritic patterns are obtained, including interface energy dendrite (IED) with the tip orientation being aligned with a laboratory frame of reference, tip splitting dendrite (TSD), and interface kinetic dendrite (IKD) with the tip orientation twisting 30◦ relative to the laboratory frame of reference, which shows the competition between interface energy anisotropy and interface kinetic anisotropy. Finally, a complete crystal growth morphology is constructed as a function of macroscopic driving force and microscopic crystalline anisotropy in the PFC model. II. PFC MODEL AND SIMULATION DETAILS
We use the two-dimensional (2D) PFC model developed by Elder et al. [12,13]. The PFC model can be derived from classic density functional theory of freezing 2 φ φ4 α + λ q02 + ∇ 2 φ + g , (1) F = d r 2 4 where φ is the time-averaged atom number density measured with respect to a reference liquid state, which is periodic in solid phase and homogeneous in liquid phase, q0 is the first-peak position of structure factor S(k) of liquid phase, being related to the lattice constant of crystal phase. Parameters α and λ deduced from the shape of polynomial approximation of S(k) determine the liquid-phase compressibility, bulk
012405-1
©2014 American Physical Society
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
modulus, and solid-liquid interface properties [19–21]. The dynamic equation is given by Swift-Hohenberg model of pattern formation, δF ∂φ = 2 , ∂τ δφ
crystal growth velocity conveniently. The position of each atom is identified as the local maximum of density field. The atom configuration of crystal growth is visualized in terms of the position of atom by using the visualization software [22].
(2)
where is related to the diffusion coefficient of atoms [13]. By defining x = rq0 , ψ = φ λqg 4 , = − λqα 4 , F˜ = λ2gq 6 F , and 0
0
and corresponding dimensionless dynamic equation reads
(4)
where corresponds to crystalline anisotropy, which is chosen to be 0.1 ∼ 0.5, indicating an increasing crystalline anisotropy [14]. The crystal growth begins with a circular nucleus located at the center of simulation box and surrounded by homogenous liquid. The initial density of the nucleus is given by one-mode approximation which reads √ √ ψ = A[cos(qx) cos(qy/ 3) − cos(2qy/ 3)/2] + ψ¯ 0 , (5) where ψ¯ 0 is the initial density of liquid phase, which determines supersaturation as σ = (ψ¯ 0 − ψ¯ Le )/(ψ¯ Se − ψ¯ Le ), where ψ¯ Se and ψ¯ Le are, respectively, the solid phase and liquid phase lines of the 2D hexagonal-liquid phase coexistence region in the 2D PFC phase diagram [13], and A is the amplitude of density waves, and q is the magnitude √ of principal reciprocal lattice vectors with the value of 3/2. Thus, the lattice constant of hexagonal structure is given by a = 2π/q. The radius of initial nucleus is chosen to be 5a empirically. The dynamic equation is solved through semi-implicit Fourier spectral method. The dimensionless grid size is dx = π/4, and time step is dt = 0.75. The finite size effect is an important factor that may affect the diffusion field of mass and energy in front of solid-liquid interface. Finite-size effect is expected to be negligible for the growth of ECS with very low velocity. For crystal growth with high velocities, finite-size effect is limited by implementing large simulation boxes. The rectangular grid with 2400 × 2400 is used for low-velocity growth and diffusionless growth, while the rectangular grid with 3600 × 3600 is used for dendritic growth. In all cases, the density far from solid-liquid interface front is maintained to be the initial ψ¯ 0 values. In order to obtain crystal growth velocity, we firstly filter the periodic oscillating density field through iterative adjacent averaging procedure, k+1 ψ¯ (i,j ) =
1
k ψ¯ (i±n,j ±m) /9,
A. Crystal growth morphology
0
dimensionless time t = λq06 τ , the dimensionless free energy functional reads as ψ ψ4 F˜ = d x [− + (∇ 2 + 1)2 ]ψ + , (3) 2 4 ∂ψ = ∇ 2 (δ F˜ /δψ), ∂t = ∇ 2 {[− + (∇ 2 + 1)2 ]ψ + ψ 3 },
III. RESULTS AND DISCUSSION
(6)
n,m=0
where k is the iteration time. The position of solid-liquid interface is defined by the isoline of smoothed density with the value of (ψ¯s + ψ¯ l )/2, where ψ¯s and ψ¯ l are the smoothed density of solid phase and liquid phase, respectively. By tracking the position of solid-liquid interface, we obtained
The crystal growth morphologies are summarized in Fig. 1. ¯ are chosen based on the PFC The simulation parameters ( ,ψ) phase diagram, as shown in Fig. 1(s). The increasing and ψ¯ correspond to the increasing crystalline anisotropy and the increasing growth driving force, respectively. For a given value, the closer to the solid-phase line the ψ¯ value is, the more rapid increase in the growth velocity is obtained usually. The crystal growth is simulated for three different values, = 0.15, 0.3, and 0.4. For each case, the crystal growth undergoes four growth stages with the ψ¯ increasing. In the first stage, when ψ¯ 0 is small, the crystal growth demonstrates ECS that changes from faceted hexagonal [Fig. 1(a) for = 0.4] to faceted hexagonal with round corners [Fig. 1(f) for = 0.3], and then to rough spherical morphology [Fig. 1(m) for = 0.15]. During the faceted ECS growth, we can observe that the solid-liquid interface grows in the layer-by-layer mode which exhibits the nucleation of a new layer and subsequent growth by absorbing atoms into kinks on edges. As ψ¯ 0 increases, crystal growth enters into the second stage. In this stage, crystal growth is dominated by dendritic morphology which demonstrates rich morphological variations caused by different ψ¯ 0 and . Shown in Fig. 2 is the convergent dendritic tip velocities of IED and IKD as the time goes by, indicating that the dendritic growth has entered into the steady state. After the growths of IED and IKD reach the steady state, the dendritic tips of IED and IKD keep the constant state while the dendritic arms grow larger and larger. For the case of strong crystalline anisotropy, such as = 0.4, the preferred growth direction of dendrite tip is dominated by interface energy anisotropy, where the principal arms of dendrites [Fig. 1(b) and (c)] are aligned with the laboratory frame as denoted by the dashed line. In terms of atom arrangement (inset A) we can find that the principle arms of dendrite are aligned with the close-packed direction. Such orientation selection is dominated by the interface energy anisotropy inhered in crystalline structure, which is consistent with the surface-tension dendritic growth as observed in Hele-Shaw experiments [9]. Here, we call this type of dendrite as interface energy dendrite (IED). When the interface energy anisotropy is of intermediate magnitude, such as = 0.3, the dendritic growth morphology changes from IED [Fig. 1(g)] to tip splitting dendrite (TSD) [Fig. 1(h) and (i)], and then to the dendrite with principal arms twisting 30◦ in relative to the laboratory frame [Fig. 1(j)]. The principle direction of dendrite tip is aligned with the open-packed direction, as seen from inset B. Such 30◦ angle twist has been observed in experiments of 2D hexagonal lattice crystal growth confined by two slices of parallel glasses, which is attributed to the competition of interface energy anisotropy and interface kinetic anisotropy [3]. As the experimental and theoretical results have proved, the interface
012405-2
PHASE-FIELD-CRYSTAL SIMULATION OF . . .
PHYSICAL REVIEW E 89, 012405 (2014)
FIG. 1. (Color online) The crystal growth morphologies simulated under different ψ¯ 0 and with increasing growth driving force and increasing crystalline anisotropy, respectively. Patterns in (a)–(p) are obtained through performing the simulations denoted by crosses in the PFC phase diagram (s). (q) and (r) are the patterns of (k) and (j) in reciprocal space, respectively. Inset figures A–E at the bottom are the density map of the parts enclosed by square box in (c), (j), (k), (l), and (m), respectively. The 2400 × 2400 simulation grid is used for the growth of ECS, CSC, and PSC, and 3600 × 3600 for dendrite patterns. The crystals shown in (a), (f), (m) contain around 20000 atoms, in (b), (c), (g), (h), (i), (j), (n) around 80000 atoms, and around 35000 atoms in (d), (e), (k), (l), (o), (p).
kinetics factor will become dominant for high crystal growth velocity [3,9,10]. Considering the twisting 30◦ dendrite tip orientation is caused by interface kinetics anisotropy, we call the dendrite as interface kinetic dendrite (IKD) that is similar to the surface-kinetic dendrite observed in Hele-Shaw experiments [9].
FIG. 2. (Color online) The dendritic tip velocity of IED in Fig. 1(g) and IKD in Fig. 1(j) growth versus time.
When the interface energy anisotropy is too weak to sustain the pattern of IED, such as = 0.15, the dendritic growth [Fig. 1(n)] is mainly of the IKD mode. This indicates that the orientation of dendritic growth can be controlled by interface kinetic anisotropy when the interface energy anisotropy is weak enough. The crystal growth enters into the third stage, as ψ¯ 0 increases further. During this stage, being independent of the values, the crystal growth retrieves to the isotropic growth. The grown crystal is of compact spherical shape that we call as compact sphere-shape crystal (CSC), as shown in Fig. 1(d) for = 0.4, Fig. 1(b) for = 0.3, and Fig. 1(o) for = 0.15. The ECS, dendrite, and CSC are of single crystal structure, which is proved by the regular hexagonal diffraction lattice in reciprocal space, as shown in Fig. 1(q). With ψ¯ 0 increasing further, crystal growth enters into the forth stage. As shown in Fig. 1(e) for = 0.4, Fig. (l) for = 0.3, and Fig. 1(p) for = 0.15, the obtained crystal morphology is of polycrystal structure, which is proved by the diffraction circle pattern in reciprocal space, as shown in Fig. 1(r). The polycrystal structure has lost crystalline anisotropy and therefore produces isotropic spherical morphology called as polycrystal sphere-shaped crystal (PSC) here.
012405-3
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
B. Crystal growth velocity
C. Density distribution
Through analyzing crystal growth velocity, we can find some common growth regimes behind different morphologies shown in Fig. 1. As shown in Fig. 3, the low-velocity growth regime A, the intermediate-velocity growth regime B, and the high-velocity growth regime C can be found for all cases of = 0.4,0.3, and 0.15. There also exists a transition zone T between B and C for the cases of = 0.3 and = 0.15. The T regime can be regarded as the signal of the competition of interface energy anisotropy and interface kinetics anisotropy. The relationship between growth morphology and velocity can be found, i.e., ECS grows in A regime, IED and TSD grow in B regime, IKD grows in T regime (partly in B regime), and CSC grows in C regime (PCS also grows in C regime, not shown here). The curves of V versus ψ¯ 0 shown in Fig. 3(a)–(c) are similar to the experimental curves of growth velocity versus increasing temperature undercoolings [23,24] or supersaturation [25]. It is easy to understand that increasing φ¯ means increasing growth driving force in PFC model that is equivalent to undercoolings or supersaturation. The curve of V versus σ is shown in Fig. 3(d). For the diffusion-controlled growth, there exists a relationship V = Aσ α , as predicted by the LKT model [26]. Here, we find that the crystal growths of ECS, IED, TSD, and part of IKD can be described by V = Aσ α with α = 4.45, nevertheless the crystal growths of other IKD and CSC deviate from such relationship. Therefore, we conclude that the ECS, IED and TSD growths are of the diffusion-controlled mode, whereas CSC growth is not diffusion-controlled, and IKD growth lies between the two growth modes.
Figure 4 illustrates the influence of crystal growth velocity on density distribution. The influence of crystal growth velocity on the density diffusion is explained firstly, which shows that the diffusion layer of the low-velocity growth [IED in Fig. 4(a)] is much thicker than that of the high growth velocity [CSC in Fig. 4(b)]. Next, the influence of growth velocity on interface structure is explained in terms of interface thickness that is defined as the distance from solid surface to liquid front. For the cases of = 0.4 and 0.3, there is a sudden increase of the interface thickness, as marked by the transition point M in Fig. 4(c). We can observe that the faceted solidliquid interface becomes diffusive after M. As we suppose, the diffusive interface leaves more possibilities for atoms to attach, which is advantageous to the large crystal growth velocity. Thus, the sudden increase of interface thickness illustrates the relationship between interface structure and crystal growth velocity. In Fig. 4(d) and (e), the density distribution is illustrated in terms of solid-phase density ψ¯ S , liquid-phase density ψ¯ L , and reduced partition coefficient k = ψ¯ S /ψ¯ L at solid-liquid interface with increasing V . For ECS and IED growths, ψ¯ S and ψ¯ L are in accordance with the values of the solid-phase and liquid-phase lines in PFC phase diagram, respectively, for example, ψ¯ S = −0.315 and ψ¯ L = −0.366 at point N as denoted by the black arrow in Fig. 4(d) and (e), which are close
FIG. 3. (Color online) The crystal growth velocity V versus initial liquid density ψ¯ 0 under three different values in (a)–(c). Shown in (d) is V versus the supersaturation σ , corresponding to (b), where the solid line is the fitting curves of V = Aσ α with α = 4.45.
FIG. 4. (Color online) The density distribution at the solid-liquid interface. (a) and (b) The smoothed density profile along the direction normal to solid-liquid interface, wherein the solid, dashed, dotted, and dot-dashed lines correspond to the advancing position, (c) the thickness of the solid-liquid interface versus the crystal growth velocity V , (d) the solid-phase and liquid-phase densities ψ¯ S (smoothed value) and ψ¯ L at the solid-liquid interface versus V , wherein ψ¯ 0 is the corresponding initial liquid density, (e) the reduced solid-liquid partition coefficient k versus V , as fitted by the Aziz’s model (red line) [27].
012405-4
PHASE-FIELD-CRYSTAL SIMULATION OF . . .
PHYSICAL REVIEW E 89, 012405 (2014)
to the equilibrium values, ψ¯ Se = −0.312 and ψ¯ Le = −0.367, as predicted by the PFC phase diagram [13]. Therefore, we suppose that the growths of ECS and IED are in the under equilibrium or near-equilibrium state. However, for TSD and IKD, ψ¯ S and ψ¯ L begin to deviate from the values of the solid-phase and liquid-phase lines, respectively. The increasing crystal growth velocity causes increasing k. At the point M when the interface thickness increases suddenly, ψ¯ S approaches the initial liquid phase density ψ¯ 0 [denoted by the blue arrow in Fig. 4(d)]. This is a signal of diffusionless growth mode, i.e., the solid-liquid interface advances so fast that density-field of the front of liquid phase has insufficient time to diffuse. The insufficient diffusion ahead of advancing solid-liquid interface causes a trapping phenomenon under high growth velocities. As shown in Fig. 4(e), the tendency of k → 1 arises up when growth velocity V increases up to the region of CSC, where ψ¯ S approaches ψ¯ L . The relationship of k vs. V is fitted in terms of k = (k0e + V /VD )/(1 + V /VD ), wherein equilibrium solid-liquid partition coefficient k0e = 0.8260 ± 0.0012 and atom diffusion velocity VD = 0.0573 ± 0.0020. The obtained relationship of k vs. V in the PFC simulation is consistent with the solute trapping phenomena in solidification of alloy, as described by the Aziz’s model [27]. We have noticed that the Aziz’s model may not work under very high velocities. In the Aziz’s model, the complete trapping limit given by k(V ) = 1 is approached asymptotically as V → ∞. In contrast, the work of Sobolev [28] predicts that there is an abrupt change of k at a finite velocity. In other words, complete trapping (k = 1) occurs at a well-defined velocity. The experimental evidence supports that the Sobolev description of solute trapping is accurate, for example, splat cooling experiments of Al-Mg performed by Galenko and Herlach [29] show a change from a eutectic to supersaturated solid solution at a finite velocity. Very recently, Humadi and coworkers investigated the solute trapping in PFC model [30]. Their results show that incorporation of a higher order time derivatives in the PFC model can lead to complete trapping occurring at a finite velocity, while the PFC model with only diffusive dynamics (i.e., first-order time derivative) can lead to the solute trapping as the Aziz’s model predicted. Since in this work, we use the PFC model with only diffusive dynamics, the consistence of our relationship of k vs. V with the Aziz’s model is reasonable. The similar trapping phenomenon has been found in a previous PFC simulation, which is resulted from the trapping of vacancy-diffusion in PFC solid [17]. Our result is in accordance with the previous PFC simulations [17,30]. The trapping effect in density diffusion exerts significant effect on crystal growth morphology in the PFC simulation. Generally, the morphological instability, such as branching, dendritic, and cellular interface, is mediated through mass diffusion field in terms of a ‘thumb’ rule of diffusion-controlled growth [31]. The trapping effect inhibits density redistribution and segregation, and hence the morphology instability is suppressed greatly. This contributes to the stable solid-liquid interface of CSC and PSC growth. Despite of the same spherical shape, the CSC and PSC are of different growth mechanisms with the circular ECS growth in Fig. 1(m). The former is attributed to the stability induced by trapping effect,
while the latter is due to the weak crystalline anisotropy determined by = 0.15. D. Morphological diagram
Finally, various crystal growth morphologies are summarized as function of σ and , as shown in Fig. 5. ECS appears in the low σ regime (for example, σ < 0.6 at the bottom of the axis). The anisotropic ECS with faceted interface structure and the isotropic ECS with rough interface structure are found in the upper panel with > 0.25 and the lower panel with < 0.25, respectively. When σ is of intermediate values (for example, 0.6 < σ < 2.8 at the bottom of the axis), the crystal growth mode turns to be dendritic growth, where the specific dendrite pattern is further determined by . In the upper panel of > 0.25, the dendritic growth is dominated by interface energy anisotropy and hence IED is obtained, while in the lower panel of < 0.25, the dendritic growth is dominated by interface kinetics anisotropy and hence IKD is obtained. The TSD region is between IED and IKD regions. As σ increases further (for example, σ > 2.8 at the bottom of the axis), crystal growth enters into diffusionless mode, where CSC and PSC are obtained. Figure 5 illustrates how the morphological selection of crystal growth is determined by growth driving force and crystalline anisotropy. The similar morphological diagram has also been obtained in the Hele-Shaw experiment for the patterns of vapor-liquid interface, which is composed by the faceted shape (region I), surface-tension dendrite (region II), tip-splitting dendrite (region III), and kinetic dendrite (IV) [9]. Our PFC morphological diagram elucidates more morphologies because the growth driving force and crystalline anisotropy span broader ranges. When the crystalline anisotropy (i.e., the y-axis variable) is weaker, as shown in the left-lower-half panel ( < 0.25, σ < 0.2) of Fig. 5, crystal growth directly enters into the IKD region after non-faceted ECS growth region with driving force increasing, and therefore no IED and TSD is obtained. Besides, the PFC simulations explore the morphological diagram when the growth driving forces
FIG. 5. (Color online) The phase-diagram of crystal growth morphology as function of PFC parameter and supersaturation σ that describe the crystalline anisotropy and the growth driving force, respectively.
012405-5
TANG, YU, WANG, LI, WANG, GUO, AND ZHOU
PHYSICAL REVIEW E 89, 012405 (2014)
is very large, demonstrating CSC grows at high velocity and PSC at higher velocity. Larger parameter spaces in the PFC model allow us to obtain a full picture of morphological selection map of the crystal growth from near-equilibrium to far-from-equilibrium. Crystal growth under very strong anisotropy, = 0.75, is studied in previous PFC simulations, and shows faceted ECS and IED growth under diffusion-controlled mode and compact hexagon-shaped crystal growth under diffusionless mode. Here, a large simulation space with ranging from 0.1 to 0.5 is explored. More PFC growth morphologies are obtained in the range, such as TSD, IKD, and PSC, enabling us to summarize the dendritic growth at intermediate velocities and the spherical crystal growth at high velocities systematically.
CSC at high velocities, and PSC at extremely high velocities. Furthermore, the competition between interface energy anisotropy and the interface kinetics anisotropy causes the orientation of dendritic arms twisting from the close-packed direction to the open-packed direction, i.e., from IED to IKD. The TSD appears under the interplay of interface energy dominated growth and interface kinetics dominated growth. At high velocities, the trapping effect in density diffusion suppresses interface instability, which is responsible for the isotropic growth of CSC and PSC. Such a broad range of PFC crystal growth morphologies reflects the underling relationship between the morphology, velocity and density diffusion of the morphological transition during crystal growth. ACKNOWLEDGMENTS
In summary, our PFC simulations present a full picture of the velocity-dependent crystal growth morphology from near-equilibrium to far-equilibrium regime, including an ECS crystal at low velocities, the dendrite at intermediate velocities,
This work has been supported by Nature Science Foundation of China (Grant Nos. 10974228 and 51071128), National Basic Research Program of China (Grant No. 2011CB610401), and EU FP7, PHASEFIELD under no. 247504, Excellent Doctorate Foundation of Northwestern Polytechnical University.
[1] J. S. Langer, Rev. Mod. Phys. 52, 1 (1980). [2] R. P. Liu, T. Volkmann, and D. M. Herlach, Acta Mater. 49, 439 (2001). [3] P. Oswald and P. Pieranski, Smectic and Columnar Liquid Crystals (Taylor and Francis, New York, 2006). [4] G. P. Ivantsov, Dokl. Akad. Nauk. SSSR 58, 567 (1947). [5] W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 35, 444 (1964). [6] A. Barbieri and J. S. Langer, Phys. Rev. A 39, 5314 (1989). [7] D. Kessler, J. Koplik, and H. Levine, Adv. Phys. 37, 255 (1988). [8] M. Ben Amar and E. Brener, Phys. Rev. Lett. 71, 589 (1993). [9] E. Ben-Jacob and P. Garik, Nature (London) 343, 523 (1990). [10] E. Ben-Jacob, N. Goldenfeld, J. S. Langer, and G. Sch¨on, Phys. Rev. Lett. 51, 1930 (1983). [11] F. Liu and N. Goldenfeld, Phys. Rev. A 42, 895 (1990). [12] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys. Rev. Lett. 88, 245701 (2002). [13] K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004). [14] R. Backofen, A. Ratz, and A. Voigt, Philos. Mag. Lett. 87, 813 (2007). [15] S. Tang, R. Backofen, J. Wang, Y. Zhou, A. Voigt, and Y. M. Yu, J. Cryst. Growth 334, 146 (2011). [16] S. Tang, Z. J. Wang, Y. L. Guo, J. C. Wang, Y. M. Yu, and Y. H. Zhou, Acta Mater. 60, 5501 (2012).
[17] G. Tegze, G. I. T´oth, and L. Gr´an´asy, Phys. Rev. Lett. 106, 195502 (2011). [18] G. Tegze, L. Gr´an´asy, G. I. T´oth, J. F. Douglas, and T. Pusztai, Soft Matter 7, 1789 (2011). [19] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phys. Rev. B 75, 064107 (2007). [20] K.-A. Wu and A. Karma, Phys. Rev. B 76, 184107 (2007). [21] A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, Phys. Rev. E 80, 031602 (2009). [22] J. Li, Model. Simul. Mater. Sci. Eng. 11, 173 (2003). [23] R. Willnecker, D. M. Herlach, and B. Feuerbacher, Phys. Rev. Lett. 62, 2707 (1989). [24] C. Panofen and D. M. Herlach, Appl. Phys. Lett. 88, 171913 (2006). [25] F. C. Frank, Contemp. Phys. 23, 3 (1982). [26] Lipton, W. Kurz and R. Trivedi, Acta Metall. 35, 957 (1987). [27] M. J. Aziz, J. Appl. Phys. 53, 1158 (1982). [28] S. L. Sobolev, Phys. Rev. E 55, 6845 (1997). [29] P. K. Galenko and D. M. Herlach, Phys. Rev. Lett. 96, 150602 (2006). [30] H. Humadi, N. Ofori-Opoku, N. Provatas, and J. J. Hoyt, JOM 65, 1103 (2013). [31] R. Trivedi and W. Kurz, Acta Metall. 34, 1663 (1986).
IV. SUMMARY
012405-6
