Amplitude expansion of the binary phase field crystal model

Amplitude expansion of the binary phase field crystal model

K. R. Elder Department of Physics, Oakland University, Rochester, MI 48309    Zhi-Feng Huang Department of Physics and Astronomy, Wayne State University, Detroit, MI 48201    Nikolas Provatas Department of Materials Science and Engineering and Brockhouse Institute for Materials Research, McMaster University, Hamilton, ON, Canada L8S-4L7
August 5, 2019

Amplitude representations of a binary phase field crystal model are developed for a two dimensional triangular lattice and three dimensional BCC and FCC crystal structures. The relationship between these amplitude equations and the standard phase field models for binary alloy solidification with elasticity are derived, providing an explicit connection between phase field crystal and phase field models. Sample simulations of solute migration at grain boundaries, eutectic solidification and quantum dot formation on nano-membranes are also presented.

81.10.Aj, 81.15.Aa, 81.30.-t, 61.72.Bb

I Introduction

The development and analysis of continuum field theories to model complex non-equilibrium structures or spatial patterns has made a tremendous impact in many areas of research in condensed matter and materials physics. A central idea in the development of such models is the recognition that the patterns are controlled by the type and interaction of defects that define the patterns. For example in spinodal decomposition surfaces between different atomic species interact through diffusion in the bulk and control the phase segregation process. In block copolymer systems Harrison et al. (2000) disclinations interact through elastic fields and control the ordering of lamellar phases. While it is desirable that models of these processes be derived from some fundamental atomic theory, they are frequently phenomenologically proposed. Classic phenomenological models include the Ginzburg-Landau theory of superconductivity Ginzburg and Landau (1950) and the Cahn-Hilliard-Cook theory of phase segregation Cahn and Hilliard (1958); Cook (1970).

Several years ago a phase field model of crystallization was phenomenologically proposed by exploiting the properties of free energy functionals that are minimized by periodic fields. In crystallization this field is interpreted as the atomic number density (), which is uniform in a liquid phase and is typically periodic in the solid phase. By incorporating elasticity, dislocations and multiple crystal orientations, such functionals naturally incorporate the type and interaction of the defects that control many crystallization phenomena. This so-called phase field crystal (PFC) model Elder et al. (2002); Elder and Grant (2004) has been used to study glass formation Berry et al. (2008a), climb and glide dynamics Berry et al. (2006), pre-melting at grain boundaries Berry et al. (2008b); Mellenthin et al. (2008), epitaxial growth Huang and Elder (2008); Yu et al. (2009), commensurate/incommensurate transitions Achim et al. (2006); Ramos et al. (2008), sliding friction phenomena Achim et al. (2009) and the yield strength of polycrystals Elder et al. (2002); Elder and Grant (2004); Hirouchi et al. (2009). More recently a simple binary phase field crystal modeled was developed Elder et al. (2007) that couples the features of the PFC model of a pure materials with a concentration field so that eutectic growth, spinodal decomposition and dendritic growth can also be studied. This model can be linked with classical density functional theory (CDFT) (and the parameters entering CDFT), although the approximations needed to connect PFC with CDFT are quite drastic. As shown in recent studies on Fe Jaatinen et al. (2009) and colloidal systems van Teeffelen et al. (2009), CDFT predicts that is very sharply peaked in space (at atomic lattice positions) while the PFC solutions are almost sinusoidal in space. Nevertheless these same studies indicate that the parameters entering PFC models can be adjusted to match experimental quantities.

While the periodic structure of PFC models is essential for describing elasticity and plasticity, it is very inconvenient for numerical calculations. For example PFC simulations typically require (where is dimension) spatial grid points per atomic lattice site. Obviously this limits the method to relatively small systems, although several new computational algorithms have been developed that can significantly extend the applicability of both pure Backofen et al. (2007); Cheng and Warren (2008); Wise et al. (2009); Hu et al. (2009) and binary PFC models Tegze et al. (2009). To alleviate this limitation an amplitude expansion of the PFC model was developed by Goldenfeld et al. Goldenfeld et al. (2006a); Goldenfeld et al. (2005); Goldenfeld et al. (2006b). To understand the idea behind such expansions it is useful to consider a one-dimensional equilibrium state of the form , where the amplitude, , is zero in the liquid and finite in the solid state. While the field varies rapidly in space, on a length scale set by , the amplitude is uniform in crystalline regions and only varies near dislocations and liquid solid surfaces. Deformations of the crystal lattice can be represented by spatial variations in the phase of the amplitude. Using this amplitude representation, Athreya et al. Athreya et al. (2007) were able to apply adaptive mesh refinement to simulate grain growth on micron scales while simultaneously resolving atomic scale structures at interfaces. This remarkable achievement suggests that the development of amplitude expansions is very promising for computational materials research. More recently this expansion has been extended to include spatial variations in the average number density in two and three dimensions Yeon et al. (2009).

In addition to greatly increasing computational efficiency, amplitude representations of PFC models can also be exploited to establish a link between PFC type models and traditional phase field models. This link provides insight into the specific terms that enter the bulk free energy and gradient energy coefficients of traditional phase field models Wu et al. (2006); Wu and A.Karma (2007); Majaniemi and Provatas (2009); Provatas and Majaniemi (2009). Since the relationship between the parameters that enter phase field models and sharp interface models are well established Karma and Rappel (1998); Elder et al. (2001); Echebarria et al. (2004), the relationship between parameters in PFC models to sharp interface models can then be established.

In this paper, amplitude expansions are developed for triangular (2d), BCC and FCC crystal symmetries. The method of multiple scales expansion methods employed by Yeon et al. Yeon et al. (2009) is used. In the small deformation limit the expansions are shown to reproduce standard phase field models of solidification and eutectic growth which incorporate elasticity and solute segregation effects. Sample simulations of grain boundary segregation, eutectic solidification and quantum dot growth on nano-membranes are also presented to illustrate the flexibility of the amplitude equations.

Ii Binary Phase Field Crystal Model

As discussed above, the binary-alloy PFC model developed recently Elder et al. (2007) can incorporate the important features of solidification, phase segregation and solute expansion in alloy systems, in addition to the elasticity, plasticity, and multiple crystal orientations that characterize the crystalline state. For an alloy consisting of and atoms the model can be written in the case of equal atomic mobilities of the constituents as


where , () is mobility of the () species and is time. The field is the dimensionless number density difference given by , where , and are the atomic number densities of atoms, atoms and a reference liquid respectively. The field , plays the role of a concentration field. The free energy is given by


where and and are constants. The parameters , and depend on and in the simplest non-trivial case can be set to , and . Details of the parameters entering this model are described in Ref. Elder et al. (2007). Here is the solute expansion coefficient and to further simply calculations it will be assumed to be small. In the small limit the free energy can be rewritten as




, , and . The equations of motion are then,


where a time scale has been adopted for the rescaling.

In the following section the equations of motion will be developed for the slowly varying amplitudes that describe various crystalline systems. More formally, for any given periodic structure the density field can be expanded as


which separates the “slow”-scale complex amplitudes and average density field from the underlying small-scale crystalline structure that is characterized by , where are the principle reciprocal lattice vectors and are integers. The corresponding solutions for are then relatively smooth and can be approximated using only a few amplitudes. In what follows, model equations will be developed for only the lowest order amplitudes that are needed to reconstruct a given crystal symmetry and defect structures relevant for controlling elastic and plastic effects in solidification and impurity segregation.

Recently many efforts have been devoted to developing amplitude expansion for various physical systems. The central assumption of these approaches is that the amplitudes vary on scales much larger that the short (or “fast”) atomic spacing scale. Formally a small parameter can then be introduced that represents the ratio of the two scales and an expansion in this variable can be performed. To apply this analysis to Eqs. (5) and (6), both the amplitudes () and concentration field are assumed to be slow variables. For technical details of this multiple-scale analysis the reader can refer to Yeon et al. Yeon et al. (2009) and the references therein. For simplicity in this paper the average atomic density (i.e., ) will be assumed as constant and zero, since miscibility gaps between liquid and solid phase can be accounted for (to some extent) by a miscibility in . Moreover, noise dynamics will be neglected here. More complete analysis involving dynamic variation of , noise effects, as well as the general case of different atomic mobilities will be presented elsewhere Huang et al. (2009a).

Iii Amplitude expansion for triangular symmetry

In two dimensions the free energy given in Eq. (3) is minimized by a triangular lattice. The corresponding principle reciprocal lattice vectors are given by


where is the equilibrium wave number. To construct the minimal model of a triangular lattice only the lowest order reciprocal lattice vectors are needed, which correspond to , and . Following the standard methods of multiple-scale expansion Cross and Hohenberg (1993); Yeon et al. (2009) the following equations of motion can be derived




To further simply calculations a long wavelength approximation will be made such that . Unfortunately a similar long wavelength approximation can’t be made for (e.g., replacing by ) as then the equations would not be rotationally invariant. In this limit the equations of motion become


These dynamics can alternatively be written in a form of Model C type (in the Halperin-Hohenberg classification scheme Hohenberg and Halperin (1977)), i.e.,




with referring to the complex conjugate.

iii.1 Small deformation limit

To gain insight into the above results and make connection with traditional phase field models it is useful to rewrite in the form . In this case the magnitude of distinguishes between liquid () and solid phases (), while is the displacement vector introduced in continuum elasticity theory Landau and Lifsitz (1999) to describe displacement of atoms from a perfectly ordered crystal lattice. Substituting this expression into Eq. (13) and taking the long wavelength limit (done by retaining derivatives only to second order) gives,


where is the linear strain tensor. The first set of terms (defined by the brackets) is remarkably similar to standard phase field models developed for eutectic and dendritic solidification Elder et al. (1994, 2000). The polynomial in gives a first order transition from a liquid () to solid phase (). The polynomial in is the typical ‘’ free energy used in the Cahn-Hilliard-Cook type models. The coupling term (note: is negative) can lead to phase segregation at low temperatures when becomes large. The second and third set of terms account for surface and linear elastic energy respectively. From the form of the third term it is straightforward Landau and Lifsitz (1999) to derive the elastic constants in dimensionless units, i.e., and .

Finally the last set of terms couples the concentration field to the liquid/solid order parameter when the atomic species have a different size (i.e., ). The term, , implies preferential phase segregation to liquid/solid surface, grain boundaries and dislocations (i.e., regions in which varies in space). Dynamically this term is responsible for solute migration at grain boundaries and solute drag. The last term implies a coupling between strain and concentration as should be expected when the atomic species have different sizes. It may appear unusual that the free energy functional depends on the sign of the displacement gradients (via ); however this sign determines whether there is a local compression or expansion of the lattice which would favor a specific atomic species based on the sign of solute expansion coefficient .

iii.2 Equilibrium phase diagram

The equilibrium phase diagram can be evaluated by considering and constant and a bulk compression to account for solute expansion, such that . In this limit the free energy per unit area () becomes


Minimizing with respect to gives . As expected the contraction/expansion of the lattice is controlled by . Substituting for leads to


which is then minimized with respect to , yielding


Substituting Eq. (17) into Eq. (16) yields a free energy per unit area that is only a function of . This result can be used to construct the phase diagram in an analogous manner as was done in Ref. Elder et al. (2007). Two example phase diagrams and the corresponding model parameters required to obtain each are show in Fig. 1.

Figure 1: Sample phase diagram for two dimensional triangular system with parameter set . Also and in (a) and (b) respectively. In each panel the filled regions correspond to regions of phase coexistence.

The derivation of amplitude equations and the related simplifications presented above can be readily extended to three dimensional systems. In Appendix A and B the relevant equations and sample phase diagrams for both BCC and FCC symmetries are presented.

Iv Applications

To further examine the above amplitude equations, numerical simulations of the model described by Eqs. (12) and (13) were undertaken and the results are shown in Figs. 26. In Fig. 2 simulations were conducted to study the coupling between composition and topological defects (dislocations and grain boundaries) in binary alloys with components of different atomic sizes, i.e., nonzero solute expansion coefficient . For this study a symmetric tilt grain boundary between two grains with a misorientation angle of was created by dynamically evolving an initial configuration of two perfect crystals separated by a layer of supercooled liquid. As time evolves the liquid solidifies and a grain boundary spontaneously forms. The dislocation cores that comprise the grain boundary interact with the different atomic species or alloy composition. As shown in Fig. 2 the larger (smaller) solute atoms preferentially segregates around the dislocation cores in regions of tensile (compressive) strain.

Figure 2: Solute segregation in symmetric grain boundaries with misorientation are shown for (a) and (b). The left and right panels correspond to and respectively. In the corner insets of left panels the density field is reconstructed from the amplitudes for the boxed region. In the right panels the dark (light) color corresponds to the larger (smaller) of the atomic species. The parameters for (a) and (b) correspond to Fig. 1a and Fig. 1b at respectively.

Simulation results for eutectic solidification and phase separation are presented in Fig. 3, where three small grains of different orientation are heterogeneously nucleated for the parameters used in Fig. 1b and for a “temperature” below the eutectic point. As the grains grow the lamellar concentration bands form within the grains as a result of phase separation. The relatively large lattice mismatch (roughly in equilibrium, due to a finite solute expansion coefficient) between the lamella results in the spontaneous nucleation of dislocation at the lamellar interfaces. When the grains impinge on one another (i.e., coarsening occurs) additional dislocations form, leading to complex patterns as shown in Figs. 3b and 3c. All these simulations indicate that the amplitude model established here can simultaneously describe the complex evolution of liquid/solid interfaces, grain boundaries, dislocations and interfaces between regions of different compositions.

Figure 3: Eutectic solidification for parameters in Fig. 1b at and . Panels (a), (b) and (c) correspond to dimensionless times 30000, 60000 and 105000 respectively. From left to right the columns correspond to reconstructed from (for the boxed region), , and the local free energy density. Dislocations are most easily identified as small black dots in the local free energy density.

Numerical simulations were also conducted to study islands or quantum dots formation on thin freely standing films or solid nanomembranes. Recent experiments of Ge on Si membranes Huang et al. (2009b); Kim-Lee et al. (2009) have suggested that growth on such nanomembranes strongly influences the maximum size that the strained islands can form coherently (i.e., without dislocations) and can lead to correlation or self-ordering of multiple islands. To examine this phenomenon, simulations were set up such that islands of one material were grown on a thin free standing membrane of another material, by exploiting a eutectic phase diagram (such as that shown in Fig. 1b). To initiate growth of the islands a small crystal region at was constructed coherently on top of a thin membrane at in a supercooled liquid at , with set to be . Sample simulation results are shown in Figs. 4 and 5. As can be seen in these figures the dot grows coherently with the membrane until reaching some critical size at which dislocations are nucleated at the liquid/dot/membrane junction. Perhaps more interestingly by comparing Figs. 4 and 5 it is apparent that for the thinner membrane the dots can grow to a larger size before nucleating dislocations. This result occurs as the strained quantum dots partially relax by straining the substrate membrane and thinner membranes are easier to deform than thicker ones. Such mechanism has been proved to play an important role in engineering the self assembly of thin film nanostructures such as quantum dots Huang et al. (2009b).

Figure 4: Quantum dot formation on a two atomic layer thick nanomembrane at times , and (for (a) to (c) respectively). The columns from left to right correspond to , and the local free energy density. The parameters used for this simulation are from Fig. 1a except and .

Figure 5: Quantum dot formation on a five atomic layer thick nanomembrane at times , and (for (a) to (c) respectively). The parameters used are identical to those in Fig. 4.

Figure 6: Correlated quantum dot formation on a nanomembrane at times and (for (a) to (f) respectively). The parameters used are identical to those in Fig. 4 except for a higher liquid supersaturation.

Another interesting consequence of the membrane deformation is that it locally creates favorable and unfavorable positions for the nucleation of other islands/dots, with an alternative sequence on the two sides of the membrane, as can be seen in Fig. 6 which was conducted at a higher liquid supersaturation (i.e., ) or growth rate. Once a dot is formed, e.g., on the top of the membrane, it is preferable for the next ones to nucleate under the edges rather than directly underneath the top quantum dot. After these bottom-side dots are formed they in turn create preferential regions for a new set of dots to nucleate on the top (above their edges) and the process repeats. In such a fashion the process leads to the correlated dot growth as observed in experiments of Ge or SiGe islands on Si membranes Huang et al. (2009b); Kim-Lee et al. (2009). This phenomenon can in principle be exploited to create periodic strained nanostructures that in turn produce periodically modulated band gaps.

Finally it should be noted that in the dynamics described by Eq. (12) the diffusion constant for the concentration field () is similar in magnitude in the liquid and solid phases. However it is typically the case that diffusion or atomic mobility in liquids is much larger than in solids. This difference can be incorporated by using the following dynamic equation for ,


where the mobility now depends on . The following functional form,


would for example change from the mobility of the crystal () at large to the mobility in the liquid phase () at small or zero . The parameter would control how quickly the changes from liquid to solid phases. In addition, while the binary PFC model presented in Eq. (1) assumes both atomic species have the same mobility, the concentration dependence of could also be added in the above ad hoc method suggested for the phase dependence of atomic mobility.

V Discussions and Conclusions

In this paper amplitude equations have been derived from the binary phase field crystal model for a two dimensional triangular lattice and for three dimensional BCC and FCC structures. Furthermore, the connection to standard phase field models has been established, and for small deformations the results have been shown to recover linear continuum elasticity theory and reconstruct the equilibrium phase diagrams for binary alloy systems. Sample simulations of the amplitude equations have shown that this relatively simple model can effectively model many complex phenomena and the emergent microstructures that arise, and reveal the underlying mechanisms. While these amplitude equations were derived from an atomistic model, they can in themselves be regarded as phase field models with complex order parameters. One advantage of this amplitude description is that the liquid and solid phases are easily distinguished by a relative uniform quantity , and the coupling and interaction between structure (i.e., amplitudes) and concentration can be well identified. In this regard it may be possible to extend the equations to naturally incorporate other uniform fields, such as magnetization, polarization, temperature, etc. and simultaneously include elastic and plastic deformation in polycrystalline samples.

Vi Acknowledgements

K.R.E. acknowledges support from NSF under Grant No. DMR-0906676. Z.-F.H. acknowledges support from NSF under Grant No. CAREER DMR-0845264. N. P. acknowledges support from the National Science and Engineering Research Council of Canada

Appendix A Amplitude Equations for BCC symmetry

For a BCC lattice, the principle reciprocal lattice vectors are


where is the equilibrium wave number. In a one mode approximation the above principle reciprocal lattice vectors need to be combined with the following vectors, , and . Thus the BCC structure can be represented in the usual manner, i.e.,


Repeating the calculations presented in section (III) (with the same level of approximations) gives the following complex amplitude equations:


with equations of motion for and obtained by cyclic permutations on the groups (1,2,3) and (4,5,6) in Eq. (22), while equations for and can be obtained by the similar cyclic permutations of Eq. (23). The corresponding concentration equation is given by

Figure 7: Sample phase diagram for three dimensional BCC system with same parameter set and notation as in Fig. 1 except in (a).

Once again, Eqs. (22)-(LABEL:BCCpsi) can be rewritten in a ‘model C’ type form, i.e.,




As discussed in Sec. III.1 it is interesting to replace by , yielding


Similarly the elastic constants can be derived as and . Furthermore, minimizing with respect to gives, leads to the following result for the free energy per unit volume,


which can be minimized with respect to to obtain


The corresponding phase diagram for this model is shown in Fig. 7.

Appendix B Amplitude Equations for FCC symmetry

Recently Wu Wu (2006) has shown that the basic phase field crystal model can be extended to model a FCC lattice by including an extra length scale. Replacing the operator , where for an FCC lattice and the extra factor of is introduced just for convenience, gives a revised alloy PFC free energy,


Here the parameter controls the symmetry of the phase such that for large (small) a BCC (FCC) structure is favored. The following calculations focus on the FCC phase in the limit of . Furthermore, to examine the influence of solute expansion the parameters and can be set to be and . In the small limit the free energy functional becomes,


showing the same form as Eq. (3) but with different operators


The corresponding equations of motion are also governed by Eqs. (5) and (6), with and given above.

For the FCC symmetry the principle reciprocal lattice vectors are


To construct an FCC crystal the following reciprocal lattice vectors are also required,


Unlike the triangular and BCC symmetries, the FCC lattice requires at minimum two set of vectors of different lengths, i.e., with length and with length . The density field is then expanded in the usual fashion, i.e.,


Following the standard procedure, the amplitude equations can be derived as