{'fraction_non_alphanumeric': 0.04014738339747916, 'fraction_numerical': 0.012877806759936524, 'mean_word_length': 4.467092142002393, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling Elastic and Plastic Deformations in Non-Equilibrium Processing Using Phase Field Crystals

A continuum field theory approach is presented for modeling elastic and plastic deformation, free surfaces and multiple crystal orientations in non-equilibrium processing phenomena. Many basic properties of the model are calculated analytically and numerical simulations are presented for a number of important applications including, epitaxial growth, material hardness, grain growth, reconstructive phase transitions and crack propagation.


I. INTRODUCTION
Material properties are often controlled by the complex micro-structures that form during non-equilibrium processing. In general terms the dynamics that occur during the processing are controlled by the nature and interaction of the topological defects that delineate the spatial patterns. For example, in spinodal decomposition the topological defects are the surfaces that separate regions of different concentration. The motion of these surfaces is mainly controlled by the local surface curvature and non-local interaction with other surfaces or boundaries through the diffusion of concentration. In contrast, block-copolymer systems form lamellar or striped phases and the topological defects are dislocations in the striped lattice. In this instance the defects interact through longrange elastic fields.
One method for modeling the topological defects is through the use of 'phase fields'. These can be thought of as physically relevant fields (such as concentration, density, magnetization, etc.), or simply as auxiliary fields constructed to produce the correct topological defect motion. In constructing phenomenological models it is often convenient to take the former point of view, since physical insight or empirical knowledge can be used to construct an appropriate mathematical description. In this paper the construction of a phase field model [1] to describe the dynamics of crystal growth that includes elastic and plastic deformations is described. The model differs from other phase field approaches [2,3,4] in that the model is constructed to produce phase fields that are periodic. This is done by introducing a free energy that is a functional of the local time averaged density field, ρ( r, t). In this description the liquid state is represented by a uniform ρ and the crystal state is described a ρ that has the same periodic crystal symmetry as a given crystalline lattice. This description of a crystal has been used in other contexts [5,6], but not previously for describing material processing phenomena. For simplicity this model will be referred to as the Phase Field Crystal model, or PFC model for short.
This approach exploits the fact that many properties of crystals are controlled by elasticity and symmetry. As will be discussed in latter sections, any free energy functional that is minimized by a periodic field naturally includes elastic energy and symmetry properties of the periodic field. Thus any property of a crystal that is determined by symmetry (e.g., relationship between elastic constants, number and type of dislocations, low-angle grain boundary energy, coincident site lattices, etc.) is also automatically incorporated in the PFC model.
The purpose of this paper is to introduce and motivate this new modeling technique, discuss the basic properties of the model and to present several applications to technologically important non-equilibrium phenomena. In remainder of this section a brief introduction to phase field modeling techniques for uniform and periodic field is discussed and related to the study of generic liquid/crystal transitions. In the following section a simplified PFC model is presented and the basic properties of the model are calculated analytically. This includes calculation of the phase diagram, linear elastic constants and the vacancy diffusion constant.
In Section (III) the PFC model is applied to a number of interesting phenomena including, the determination of grain boundary energies, liquid phase epitaxial growth and material hardness. In each of these cases the phenomena are studied in some detail and the results are compared with standard theoretical results. At the end of this section sample simulations of grain growth, crack propagation and a reconstructive phase transition are presented to illustrate the versatility of the PFC model. Finally a summary of the results are presented in Section (IV).

A. Uniform Fields and Elasticity
Many non-equilibrium phenomena that lead to dynamic spatial patterns can be described by fields that are relatively uniform in space, except near interfaces where a rapid change in the field occurs. Classic examples include order/disorder transitions (where the field is the sublattice concentration), spinodal decomposition (where the field is concentration) [7], dendritic growth [8] and eutectics [9]. To a large extent the dynamics of these phenomena are controlled by the motion and interaction of the interfaces. A great deal of work has gone into constructing and solving models that describe both the interfaces ('sharp interface models') and the fields ('phase field models'). Phase field models are constructed by considering symmetries and conservation laws and lead to a relatively small (or generic) set of sharp interface equations [10].
To make matters concrete consider the case of spinodal decomposition in AlZn. If a high temperature homogeneous mixture of Al and Zn atoms is quenched below the critical temperature small Al and Zn rich zones will form and coarsen in time. The order parameter field that describes this phase transition is the concentration field. To describe the phase transition a free energy is postulated (i.e., made up) by consideration of symmetries. For spinodal decomposition the free energy is typically written as follows, where f (φ) is the bulk free energy and must contain two wells one for each phase (i.e., one for Al-rich zones and one for Zn-rich zones). The second term in Eq. (1) takes into account the fact that gradients in the concentration field are energetically unfavorable. This is the term that leads to a surface tension (or energy/length) of domain walls that separate Al and Zn rich zones. The dynamics are postulated to be dissipative and act such that an arbitrary initial condition evolves to a lower energy state. These general ideas lead to the well known equation of motion; where, Γ is a phenomenological constant. The Gaussian random variable, η c , is chosen to recover the correct equilibrium fluctuation spectrum, has zero mean and two point correlation, The variable a is equal to one if φ is a conserved field, such as concentration, and is equal to zero if φ is a nonconserved field, such as sublattice concentration.
A great deal of physics is contained in Eq. (2) and many papers have been devoted to the study of this equation. While the reader is refereed to [7] for details, the only salient points that will be made here is that 1) the gradient term and double well structure of f (φ) in Eq.
(1) lead to a surface separating different phases and 2) the equation of motion of these interfaces is relatively independent of the form of f (φ). For example it is well known [10] that if φ is non-conserved the normal velocity, V n , of the interface is given by; where κ is the local curvature of the interface and A is directly proportional to the free energy difference between the two phases. If φ is conserved then the motion of the interface is described by the following set of equations [10], where µ ≡ δF/δφ is the chemical potential, d o is the capillary length, β is the kinetic under-cooling coefficient,n is a unit vector perpendicular to the interface position, D is the bulk diffusion constant and 0 + and 0 − are positions just ahead and behind the interface respectively. It turns out that Eqs. (4) and (5) always emerge when the bulk free energy contains two wells and the local free energy increases when gradients in the order parameter field are present [10]. In this sense Eqs. (4) and (5) can be thought of as generic or universal features of systems that contain domain walls or surfaces. As will be discussed in the next subsection a different set of generic features arise when the field prefers to be periodic in space. Some generic features of periodic systems are that they naturally contain an elastic energy, are anisotropic and have defects that are topologically identical to those found in crystals. A number of research groups have built these 'periodic features' into phase field models describing uniform fields. This approach has some appealing features, as one can consider mesoscopic length and time scales. But it can involve complicated continuum models. For example, in Refs. [3,4], a continuum phase-field model was constructed to treat the motion of defects, as well as their interaction with moving free surfaces. Although such an approach gives explicit access to the stresses and strains, including the Burger's vector via a ghost field, the interactions between the nonuniform stresses and plasticity are complicated, since the former constitutes a free-boundary problem, while the latter involves singular contributions to the strain, within the continuum formulation.

B. Periodic Systems
In many physical systems periodic structures emerge.
In addition many convective instabilities [14], such as Rayleigh-Bénard convection and a Margonoli instability lead to periodic structures (although it is not always possible to describe such systems using a free energy functional). To construct a free energy functional for periodic systems it is important to make the somewhat trivial observation that unlike uniform systems, these systems are minimized by spatial structures that contain spatial gradients. This simple observation implies that in a lowest order gradient expansion the coefficient of | ∇φ| 2 in the free energy (see Eq. (1)) is negative. By itself this term would lead to infinite gradients in φ so that the next order term in the gradient expansion must to be included (i.e., |∇ 2 φ| 2 ). In addition to these two terms a bulk free energy with two wells is also needed, so that a generic free energy functional that gives rise to periodic structures can be written, where K and a o are phenomenological constants. Insight into the influence of the gradient energy terms can be obtained by considering a solution for φ of the form φ = A sin(2πx/a). For this particular functional form for φ the free energy becomes, where ∆a ≡ a − a o . At this level of simplification it can be seen that the free energy per unit length is minimized when a = a o , or a o is the equilibrium periodicity of the system. Perhaps more importantly it highlights the fact the energy can be written in a Hooke's law form (i.e., E = E o + (k∆a) 2 ) that is so common in elastic phenomena. Thus a generic feature of periodic systems is that for small perturbation (e.g., compressions or expansions) away from the equilibrium they behave elastically. This feature of periodic systems will be exploited to develop models for crystal systems in the next section.

C. Liquid/Solid systems
In a liquid/solid transition the obvious field of interest is the density field since it is significantly different in the liquid and solid phases. More precisely the density is relatively homogeneous in the liquid phase and spatially periodic (i.e., crystalline) in the solid phase. The free energy functional can then be approximated as; (8) where f and G are to be determined and δρ is the deviation of the density from the average density (ρ). Under constant volume conditions δρ is a conserved field, so that the dynamics are given by; where η is a Gaussian random variable with zero mean and two point correlation, To determine the precise functional form of the operator G(∇ 2 ) it is useful to consider a simple liquid since δρ is small and f can be expanded to lowest order in δρ, i.e., where f (i) ≡ (∂ i f /∂δρ i ) δρ=0 . In this limit Eq. (9) takes the form which can be easily solved to give, where, q is the wavevector,ω q ≡ f (2) + G(q 2 ), δρ is the Fourier transform of δρ, i.e., δρ( q, t) ≡ d re i q· r δρ( r, t)/(2π) d .
and d is the dimension of space. The structure factor, S(q, t) ≡ |δρ| 2 , is then; In a liquid system the density is stable with respect to fluctuations which implies thatω q > 0. The equilibrium liquid state structure factor, S eq liq (q) = S(q, ∞), then becomes, This simple calculation indicates the phase field method can model a liquid state if the functionω q is replaced with k B T /S eq liq (q), or A typical liquid state structure factor and the correspondingω q is shown in Fig. (1). Thus G(∇ 2 ) can be obtained for any pure material through Eq. (17).
In the solid state the density is unstable to the formation of a periodic structure (i.e., to forming a crystalline solid phase) and thusω q must go negative for certain values of q. This instability is taken into account by the temperature dependence of f (2) , i.e., Thus when T > T m , w q is positive and the density is uniform. When T < T m , w q is negative (for some values of q) and the density is unstable to the formation of a periodic structure. To properly describe this state, higher order terms in δρ must be included in the expansion of f (δρ), since δρ is no longer small. Before discussing the properties of a specific choice for f (δρ) it is worth pointing out some generic elastic features of such a model. As illustrated in the Section (I B) a free energy that is minimized by a periodic structure has 'elastic' properties. The elastic constants of the system can be obtained by formally expanding around an equilibrium state in the strain tensor. If the equilibrium state is defined to be δρ o ( r) and the displacement field is u, then δρ can be written δρ( r) = δρ o ( r + u) + ǫ, where ǫ will always be chosen to minimize the free energy. Expanding to lowest order in the strain tensor gives, where C ij,kl is the elastic constant given by, In Eq. (20) the Einstein summation convection is used and u ij represents the usual components of the strain tensor, i.e., and the subscript o indicates that the derivative are evaluated at δρ = δρ o (i.e., δρ o = u ij = 0). While Eq. (20) is somewhat formal and difficult to use for a specific model it does highlight several important features. Eq. (20) shows that the elastic constants are simply related to the curvature of the free energy along given strain directions Perhaps more importantly Eq. (20) shows that C ij,kl is proportional to H which is a function of the equilibrium density field δρ o . Thus if the free energy is written such that F is minimized by a δρ o that is cubic, tetragonal, hexagonal, etc., then C ij,kl will automatically contain the symmetry requirements of that particular system. In other words the elastic constants will always satisfy any symmetry requirement for a particular crystal symmetry since C ij,kl is directly proportional to a function that has the correct symmetry. This also applies to the type or kind of defects or dislocations that can occur in any particular crystal system, since such deformation are determined by symmetry alone.
In the next section a very simple model of a liquid/crystal transition will be presented and discussed in some detail. This model is constructed by providing the simplest possible approximation for f (δρ) that will lead to towards a transition from a uniform density state (i.e., a liquid) to a periodic density state (i.e., a crystal).

II. SIMPLE PFC MODEL: BASIC PROPERTIES
In this section perhaps the simplest possible periodic model of a liquid/crystal transition will be presented. Several basic features of this model will be approximated analytically in the next few subsections. This includes calculations of the phase diagram, the elastic constants and the vacancy diffusion constant.

A. Model
In the preceding section it was shown that a particular material can be modeled by incorporating the two point correlation function into the free energy through Eq. (17). It was also argued that the basic physical features of elasticity are naturally incorporated by any free energy that is minimized by a spatially periodic function. In this section the simplest possible free energy that produces periodic structures will be examined in detail. This free energy can be constructed by fitting the following functional form for G, to the first order peak in an experimental structure factor. As an example such a fit is shown for argon in Fig. (1). At this level of simplification the minimal free energy functional is given by; In principle other non-linear terms (such as ρ 3 ) can be included in the expansion but retaining only ρ 4 simplifies calculations. The dynamics of ρ are then described by the following equation, For convenience it useful to rewrite the free energy in dimensionless units, i.e., In dimensionless units the free energy becomes, The dimensionless equation of motion becomes, (27) and (28) describe a material with specific elastic properties. In the next few sections the properties of this 'material' will be discussed in detail. As will be shown, some of the properties can be adjusted to match a given experimental system and others cannot be matched without changing the functional form of the free energy. For example the periodicity (or lattice constant) can be adjusted since all lengths have been scaled with q o . The bulk modulus can be also be easily adjusted since the free energy has been scaled with λ, u and q o . On the other hand this free energy will always produce a triangular lattice in two dimensions [5,6]. To obtain a square lattice a different choice of the non-linear terms must be made. This is the most difficult feature to vary as there are no systematic methods (known to the authors) for determining which functional form will produce which crystal symmetry. Cubic symmetry can be obtained by replacing ψ 4 with |∇ψ| 4 [16,17].
In the next few subsections the properties of this free energy and some minor extensions will be considered in one and two dimensions. The three dimensional case will be discussed in a future paper.

B. One dimension
In one dimension the free energy described by Eq. (26) is minimized by a periodic function for small values of ψ o and a constant for large values. To determine the properties of the periodic state it is useful to make a one mode approximation, i.e., ψ ≈ A sin(qx) + ψ o , which is valid in the small r limit. Substitution of this function into Eq. (26) gives, whereω q is the fourier transform of ω(∇ 2 ), i.e.,ω q = r + (1 − q 2 ) 2 . Minimizing Eq. (29) with respect to q gives the selected wavevector, q * = 1. Minimizing F with respect to A gives, . This solution is only meaningful if A is real, since the density is a real field. This implies that periodic solutions only exist when r < −3ψ 2 o , sinceω q * = r. The minimum free energy density is then; Equation (30) represents the free energy density of a periodic solution in the one mode approximation. To determine the phase diagram this energy must be compared to that for a constant state (i.e., the state for which ψ c = ψ o ) which is, To obtain the equilibrium states the Maxwell equal area construction rule must be satisfied, i.e., where and µ = ∂F /∂ψ o . Using these conditions it is straightforward to show that for r > −1/4 a periodic state is selected for |ψ o | < −r/3 and a constant state is selected when |ψ o | > −r/3. For r < −1/4, there can exist a coexistence of periodic and constant states. If the constant and periodic states are considered to be a liquid and crystal respectively then this simple free energy allows for the coexistence of a liquid and crystal, which implies a free surface. The entire phase diagram is shown in Fig. (2). It is also relatively easy to calculate the elastic energy in the one mode approximation. If a ≡ 2π/q is defined as the one-dimensional lattice parameter, then the F can be written, where u ≡ (a − a o )/a o is the strain and K is the bulk modulus and is equal to; or for the particular dispersion relationship used here, The existence of such a Hooke's law relationships is automatic when a periodic state is selected since F always increases when the wavelength deviates from the equilibrium wavelength.
C. Two dimensions

Phase Diagram
In two dimensions F is minimized by three distinct solutions for ψ. These solutions are periodic in either zero dimensions (i.e., a constant), one dimension (i.e., stripes) or two dimensions (i.e., triangular distributions of drops or 'particles'). The free energy density for the constant and stripe solutions are identical to the periodic and constant solution discussed in the preceding section. The two dimensional solution can be written in the general form, where, G ≡ n b 1 + m b 2 and the vectors b 1 and b 2 are reciprocal lattice vectors. For a triangular lattice the reciprocal lattice vectors can be written, where a is the distance between nearest neighbor local maxima of ψ (which corresponds to the atomic positions). In analogy with the one-dimensional calculations presented (see Sec. (II B)) a one mode approximation will be made to evaluate the phase diagram and elastic constants. In a two dimensional triangular system a one mode approximation corresponds to retaining all fourier components that have the same length. More precisely the lowest order harmonics consist of all (n, m) pairs such that the vector G has length 2π/(a √ 3/2). This set of vectors includes (n, m) = (±1, 0), (0, ±1), (1, −1) and (−1, 1). Furthermore since ψ is a real function the fourier coefficients must satisfy the following relationship, a n,m = a −n,m = a n,−m . In addition, by symmetry, a ±1,0 = a 0,±1 = a 1,−1 = a −1,1 . Taking these considerations into account it is easy to show that in the lowest order harmonic expansion for a triangular solution for ψ can be represented by; where A t is an unknown constant and q t = 2π/a. Substituting Eq. (37) into Eq. (23), minimizing with respect to A t and q t gives, where, q t = √ 3/2, and S is a unit area. The accuracy of this one mode approximation was tested by comparison with a direct numerical calculation for a range of r's, using 'Method I' as described in Appendix (VI). The time step (∆t) and grid size (∆x) were 0.0075 and π/4 respectively and a periodic grid of a maximum size of 512∆x × 512∆x [18] was used. A comparison of the analytic and numerical solutions are shown in Fig. (3) for a variety of values of r (ψ o was set to be √ −r/2). The approximate solution is quite close to the numerical one and becomes exact in the limit r → 0. The analytic results can in principle be systematically improved by including more harmonics in the expansion. To determine the phase diagram in two dimensions the free energy of the triangular state (i.e., Eq. (38) must be compared with the free energy of a striped state (i.e., Eq. (30)) and a constant state (i.e., Eq. (31)). In addition since ψ is a conserved field Maxwell's equal area construction must be used to determine the coexistence regions. The phase diagram arising from these calculations is shown in Fig. (4). While this figure does not look like a typical liquid/solid phase diagram in the densitytemperature plane it can be superimposed onto a portion of an experimental phase diagram. As an example the PFC phase diagram is superimposed onto the argon phase diagram in Fig. (5).

Elastic Energy
The elastic properties of the two dimensional triangular state can be obtained by considering the energy costs for deforming the equilibrium state. The free energy density associated with bulk, shear and deviatoric deformations can be calculated by considering modified forms of Eq. (37), i.e., ψ t (x/(1 + ζ), y/(1 + ζ)) (bulk), ψ t (x+ζy, y) (shear) and ψ t (x(1+ζ), y(1−ζ) (deviatoric). In such calculations ζ represents the dimensionless deformation, q t = √ 3/2 and A t is obtained by minimizing F . The results of these calculations are given below; where These results can be used to determine the elastic constants by noting that for a two dimensional system [5,19], The elastic constants are then These results are consistent with the symmetries of a two dimensional triangular system, i.e., C 11 = C 12 +2 C 44 . In two dimensions this implies a bulk modulus of B = α/2, a shear modules of µ = α/4, a Poisson's ratio of σ = 1/3, and a two dimensional (i.e., Y 2 = 4Bµ/(B + µ)) Young's modulus of Y 2 = 2α/3. Numerical simulations were conducted (using the parameters and numerical technique discussed in the previous section) to test the validity of these approximations for the bulk modulus. The results, shown in Fig. (3), indicate that the approximation is quite good in the small r limit.
These calculations highlight the strengths and limitations of the simplistic model described by Eq. (23). On the positive side the model contains all the expected elastic properties (with the correct symmetries) and the elastic constants can be approximated analytically within a one mode analysis. One the negative side, the model as written, can only describe a system where C 11 = 3 C 12 . Thus parameters in the free energy can be chosen to produce any C 11 , but C 12 cannot be varied independently.
To obtain more flexibility a term ψ 3 could be added to the free energy.

Dynamics
The relatively simple dynamical equation for ψ (i.e., Eq. (28)) can describe a large number of physical phenomena depending on the initial conditions and boundary conditions. To illustrate this versatility it is useful to consider the growth a crystalline phase from a supercooled liquid, since this phenomena simultaneously involves the motion of liquid/crystal interfaces and grain boundaries separating crystals of different orientations. Numerical simulations were conducted using the 'Method I' as described in Appendix (VI A). The parameters for these simulations were, (r, ψ o , D, ∆x, ∆t) = (−1/4, 0.285, 10 −9 , π/4, 0.0075) on a system of size 512∆x × 512∆x with periodic boundary conditions. The initial condition consisted of large random gaussian fluctuations (amplitude 0.1) covering (10x10) grid points in three locations in the simulation cell. As shown in Fig. (6) the initial state evolves into three crystallites each with a different orientation and a well defined liquid/crystal interface. The excess energy of the liquid/crystal interfaces is highlighted in Fig. (6d) where the local free energy density is plotted.
As time evolves the crystallites impinge and form grain boundaries. As can be seen in Fig. (6) the nature of the grain boundary between grains (1) and (3) is significantly different from the boundary between grains (2) and grain (1) (or (3)). The reason is that the orientation of grains (1) and (3) is quite close but significantly different from (2). The low angle grain boundary consists of dislocations separated by large distances, while the high angle grain boundary consists of many dislocations piled together. A more detailed discussion of the grain boundaries will be given in Section (III A). Even this small sample simulation illustrates the flexibility and power of the PFC technique. This simulation incorporates, the heterogeneous nucleation of crystallites, crystallites with triangular symmetry and elastic constants, crystallites of multi orientations, the motion of liquid/crystal interfaces and the creation and motion of grain boundaries. While all these features are incorporated in standard microscopic simulations (e.g., molecular dynamics) the time scales of these simulations are much longer than could be achieved using microscopic models.
One fundamental time scale in the PFC model is the diffusion time. To envision mass diffusion in the PFC model it convenient to consider a perfect equilibrium (ψ t ) configuration with one 'particle' missing. At the atomic level this would correspond to a vacancy in the lattice. Phonon vibrations would occasionally cause neighboring atoms to hop into the vacancy and eventually the vacancy would diffuse throughout the lattice. In the PFC model the time scales associated with lattice vibrations are effectively integrated out and all that is left is long time mass diffusion. In this instance the density at the missing spot will gradually increase as the density at neighboring sites slowly decreases. Numerical simulations of this process are shown in Fig. (7) using Method I (see App. (VI A)) with the parameters (r, ψ o , D, ∆x, ∆t) = (1/4, 1/4, 0, π/4, 0.0075). To highlight diffusion of the vacancy, the difference between ψ( r, t) and a perfect equilibrium state (ψ t ) is plotted in this Fig (7). The diffusion constant in this system can be obtained by a simple linear stability analysis, or Bloch-Floquet analysis, around an equilibrium state. To begin the analysis the equation of motion for ψ is linearized around ψ t , i.e., ψ = ψ t ( r) + δψ( r, t). To first order in δψ Eq. (28) becomes, where φ t = ψ t − ψ o (see Eq. (35)). The perturbation, δψ, is then expanded as follows, where,ω ≡ r To solve Eq. (46) a finite number of modes are chosen and the eigenvalues are determined. Using the modes corresponding to the reciprocal lattice vectors in the one mode approximation ((m, n) = (±1, 0), (0, ±1), (1, −1), (−1, 1) and the (0, 0) mode gives four eigenvectors that are always negative and thus irrelevant and three eigenvalues that have the form −DQ 2 . The smallest D arises from b 0,0 mode and can be determined analytically if only this mode is used (the other eigenvalues correspond to D ≈ 3, 9). Since this is the smallest D it determines the diffusion constant in the lattice. The solution is The accuracy of Eq. (47) was tested by numerically measuring the diffusion constant for the simulations shown in Fig. (7). In this calculation the envelope of profile of δψ was fit to a gaussian (Ae −r 2 /2σ 2 ) and the standard deviation (σ) was measured. The diffusion constant can be obtained by noting that the solution of a diffusion equation (i.e., ∂C/∂t = D∇ 2 C) is C ∝ e −r 2 /4Dt , i.e., σ 2 = Dt/2. In Fig. (8) σ 2 is plotted as a function of time and the slope of this curve gives, D ≈ 1.22. This is quite close the value predicted by Eq. (47) which is 1.25.

III. SIMPLE PFC MODEL: APPLICATIONS
In this section several application of the PFC model will be considered that highlight the flexibility of the model. In Sec. (III A) the energy of a grain boundary separating two grains of different orientation is considered. The results are compared with the Read-Shockley equation [20] and shown to agree quite well for small orientational mismatch. This calculation, in part, provides evidence that the interaction between dislocations is correctly captured by the PFC model, since the grain boundary energy contains a term that is due to the elastic field set up by a line of dislocations. In Sec. (III B) the technologically important process of liquid phase epitaxial growth is considered. Numerical simulations are conducted as a function of mismatch strain and show how the model naturally produces the buckling instability and nucleation of dislocations. In Sec. (III C) the yield strength of poly-(nano-) crystalline materials is examined. This is a phenomena that requires many of the features of the PFC model (i.e., multi orientations, elastic and plastic deformations, grain boundaries) that are difficult to incorporate in standard uniform phase field models. The yield strength is examined as a function of grain size and the reverse Hall-Petch effect is observed. Finally some very preliminary numerical simulations are presented in Sec. (III D) to demonstrate the versatility of the technique. This section includes simulations of grain growth, crack propagation and reconstructive phase transitions.

A. Grain Boundary Energy
The free energy density of a boundary between two grains that differ in orientation is largely controlled by geometry. In a finite size two-dimensional system the parameters that control this energy are the orientational mismatch, θ and an offset distance ∆, as shown in Fig.  (9). For small θ, θ controls the number of dislocations per unit length and ∆ controls the average core energy. For an infinite grain boundary ∆ becomes irrelevant, unless the distance between dislocation is an integer number of lattice constants (and the integer is relatively small). Nevertheless it is straightforward to determine a lower bound on the grain boundary energy in the small θ limit, by directly relating the dislocation density to θ and assuming that the dislocation cores can always find the minimum energy location. The later assumption restricts the calculation to providing a lower bound on the grain boundary energy.
For small θ, Read and Shockley [20] were able to derive an expression for the grain boundary energy, assuming the dislocation core energy was a constant independent of geometry. In two dimensions the energy/length of the grain boundary is [5], where b is the magnitude of the Burger's vector, a is the size of the dislocation core, d is the distance between dislocations, Y 2 is the two dimensional Young's modulus and E core is the energy/length of the dislocation core.
To estimate the minimum core energy it is convenient to assume the core energy is proportional to the size of the core [5], i.e., E core = Ba 2 , where B is an unknown constant. The total energy/length then becomes To obtain a lower bound on F/L the unknown parameter B is chosen to minimize F/L, i.e., B is chosen to satisfy, d(F/L)/da = 0, which gives Ba 2 = b 2 Y 2 /16πd. Thus the free energy per unit length is, Furthermore, from purely geometrical considerations, the distance between dislocations is d = a/ tan(θ), where θ is the orientational mismatch. Finally in the small angle limit (tan(θ) ≈ θ) Eq. (50) reduces to, where the dislocation core size b was assumed to be equal to the lattice constant a.
To examine the validity of Eq. (51) the grain boundary energy was measured as a function of angle. In these simulations numerical Method I (see App. (VI A)) was used with the parameter set (r, ψ o , D, ∆x, ∆t) = (−4/15, 1/5, 0, π/4, 0.01). The initial condition was constructed as follows. On a periodic grid of size L x × L y , a triangular solution (i.e., Eq. (37)) for ψ was constructed in one orientation between 0 < x < L x /4 and 3L x /4 < x < L x . In the center of the simulation (i.e., L x /4 < x < 3L x /4) a triangular solution of a different orientation was constructed. A small slab of supercooled liquid was placed between the two crystals so as not to influence the nature of the grain boundary that emerged. The systems were then evolved for a time of t = 10, 000, after which the grain boundary energy was measured. Small portions of sample configurations are shown in Fig. (10) for θ = 5.8 o and θ = 34.2 o (the grain boundary energy is symmetric around 30 o ). As expected the Read-Shockley description of a grain boundary is consistent with the small angle configuration. In contrast the large angle grain boundary is much more complicated and harder to identify individual dislocations.
The measured grain boundary energy is compared with Eq. (51) in Fig. (11). As expected Eq. (51) provides an adequate description for small angles but not for large angles. The Read-Shockley equation does fit the measured result for all θ reasonably well if coefficients that enter the equation are adjusted, as has been observed in experiment [21,22]. This fit is shown in Fig. (12).
dimension since another degree of freedom exists. This degree of freedom can be visualized by considering taking one of the crystals shown in Fig. (9) and rotating it out of the page. The extra degree of free can lead to interesting phenomena, such as coincident site lattices that significantly alter the grain boundary energy. The PFC model should provide an excellent tool for studying such phenomena since it is purely a geometrical effect that is naturally incorporated in the PFC approach.

B. Liquid Phase Epitaxial Growth
Liquid phase epitaxial growth is a common industrial method used to grow thin films that are coherent with a substrate. The properties of such films depend on the structural integrity of the film. Unfortunately flat defectfree heteroepitaxial films of appreciable thickness are often difficult to grow due to morphological instabilities induced by the anisotropic strain arising from the mismatch between film and substrate lattice constants. Consequently, there has been a tremendous amount of scientific effort devoted to understanding the morphological stability of epitaxially grown films [1,3,4,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43].
The stability and resulting structural properties of epitaxial films are often compromised by at least two distinct processes which reduce the anisotropic strain. In one process, small mounds or ridges form as the surface buckles or corrugates to reduce the overall strain in the film. This instability to buckling can be predicted by considering the linear stability of an anisotropically strained film as done by Asaro and Tiller [23] and Grin-  Figures a), b), c) and d) correspond to times t=150, 300, 450 and 600 respectively. The grey-scale is proportional to the local density (i.e., ψ) in the film and the liquid. The substrate is highlighted by a darker grey background. To highlight nucleated dislocations, small white dots were placed on atoms near the two dislocation cores that appear in this configuration.
feld [24]. The initial length scale of the buckling is determined by a competition between the reduction in overall elastic energy which prefers mounds and surface tension and gravity both of which favor a flat interface. Another mechanism that reduces strain is the nucleation of misfit dislocations which can occur when the energy of a dislocation loop is comparable with the elastic energy of the strained film. Matthews and Blakeslee [42] and many others [27,28,29,30,31,32] have used various arguments to provide an expression for the critical height at which a flat epitaxially grown film will nucleate misfit dislocations.
The two mechanisms are often considered separately but it is clear that surface buckling can strongly influence the nucleation of misfit dislocations. Typically, as the film begins to grow, it will deform coherently by the Asaro-Tiller-Grinfeld instability. This leads initially to a roughly sinusoidal film thickness with a periodicity close to the most unstable mode in a linear analysis. As time increases, the sinusoidal pattern grows in amplitude and develops cusps or local regions of high curvature [33,34,35,36] with a periodicity similar to that of the initial instability although some coarsening may occur [3,35,36]. Eventually, the stress at the cusps become too large and a periodic array of misfit dislocations appear which reduces the roughness of the film. These dislocations eventually climb to the film/substrate interface.
The purpose of this section is to illustrate how the PFC method can be exploited to examine surface buckling and dislocation nucleation in liquid phase epitaxial growth. Modeling this process requires a slight modification of the model to incorporate a substrate that has a different lattice constant than the growing film. This can be accomplished by changing the operator ω given in Eq. (27) to be, where the parameter q controls the lattice constant of the growing film and is set to one in the substrate. To incorporate a constant mass flux the field ψ was fixed to be, ψ ℓ at a constant distance (L = 100∆x) above the film. Numerical simulations were conducted using Method I (see App. (VI A)) for the parameters (r, ψ ℓ , ∆x, ∆t) = (−1/4, .282, .785, 0.0075). The width of the film grown was L x = 8192∆x, corresponding to a width of roughly 900 particles. The initial condition was such that eight layers of substrate atoms resided at the bottom of the simulation cell with a supercooled (r = −1/4,ψ ℓ = .282) liquid above it. A small portion of a simulation is shown in Fig. (13), for q = 0.93. As can be seen in this figure, and in Fig. (14) the film initially grows in a uniform fashion before becoming unstable to a buckling or mounding instability. The film then nucleates dislocations in the valleys where the stress is the largest. After the dislocations nucleate the liquid/film interface grows in a more regular fashion. To highlight the local elastic energy, the free energy is plotted in Fig. (14). As can be seen in this figure, elastic energy builds up in the valleys during the buckling instability and is released when dislocations appear. The behavior of the liquid/film interface was monitored by calculating the average interface height and width. Both quantities are plotted in Fig. (15). The data shown in this figure is representative of all simulations conducted at different mismatch strains, but the precise details varied from run to run. In all cases the width initially fluctuates around a * /2 (where a * is the thickness of a film layer) during the 'step by step' growth. The average width then increases during buckling and decreases when dislocations nucleate. While these quantities are difficult to measure in situ there is experimental evidence for this behavior in SiGe/Si heterostructures [43].
Assigning a value to the critical height, H c at which dislocations nucleate is very subjective. Typically a first wave of dislocations is nucleated at a density that is determined by the buckling instability. Since this is not the correct density to reduce the strain to zero a subsequent buckling and dislocation occurs above the first wave. To complicate manners the nucleated dislocations climb towards the substrate/film interface. To illustrate these points the dynamics of a sample distribution of defects is shown as function of height in Fig. (16). As can be seen in this figure the first wave of dislocations appears roughly between a film height of six and thirteen layers. The critical height, as defined in the preceding paragraph, was calculated as a function of mismatch strain, ǫ = (a f ilm − a substrate )/a substrate . The equilibrium lattice constant in the film a f ilm was obtained by assuming it was directly proportional to 1/q (note, in the one mode approximation, a = 2π/( (3)q/2)) and determining the constant of proportionality by interpolating to where the critical height diverges. The numerical data was compared with the functional form proposed by Matthews and Blakeslee [42] , i.e., in Fig. (17). This comparison indicates that the data is consistent with a linear relationship between ǫ and [1 + log(H c /a * )]/(H c /a * ), where the constant of proportionality depends on whether a compressive or tensile load is applied to the substrate.

C. Material Hardness
It is well known that mechanical properties of materials depend crucially on the microstructure and grain size. For example, Hall and Petch [44,45] calculated that for large grain sizes, the yield strength of a material is inversely proportional to the square of the average grain radius. This result is due to the pileup of dislocations at grain boundaries and has been verified in many materials including Fe alloys [46,47,48], Ni [49], Ni-P alloys [50], Cu [51] and Pd [51]. However, for very small grain sizes the Hall-Petch relationship must break down, since the yield strength cannot diverge. Experimentally it is found that materials "soften" at very small grain sizes, such that the yield strength begins to decrease when the grain sizes become of the order of tens of nanometers. This 'inverse' Hall-Petch behavior has been observed in in Ni-P alloys [50], Cu and Pd [51] and molecular dynamics experiments [53,54]. Determining the precise the crossover length scale and mechanisms of material breakdown has become increasingly important in technological processes as interest in nano-crystalline materials (and nano-technology in general) increases.
The purpose of this section is to demonstrate how the PFC approach can be used to study the influence of grain size on material strength. In these simulations a poly-crystalline sample was created by heterogenous nucleation (see Sec. (III D 1) for details) in a system with periodic boundary conditions in both the x and y directions. A small coexisting liquid boundary of width 200∆x was included on either side of the sample. To apply a strain the particles near the liquid/crystal boundary (i.e., within a distance of 16∆x) were 'pulled' by coupling these particles to a moving field that fixed the particle positions. Initially the system was equilibrated for a total time of 4000 (2000 before the field was applied and 2000 after). An increasing strain was modeled by moving the field every so many time steps in such a manner that the size of the polycrystal increased by 2∆x. To facilitate relaxation, ψ was extrapolated to the new size after every movement of the external field. The parameters of the simulations to follow were (r, ψ sol , ψ liq , L x , L y , ∆x, ∆t) = (−0.3, 0.312, 2048∆x, 2048∆x, 0.377, 0.79, 0.05) and the pseudo-spectral numerical method described in App. (VI B) was used.
A sample initial configuration is shown in Fig. (18). This particular sample contains approximately 100 grains with an average grain radius of 35 particles. As can be seen in this figure there exists a large variety (i.e., distribution of mismatch orientations) of grain boundaries as would exist in a realistic poly-crystalline sample. The same configuration is shown after it has been stretched in the x direction in Figs. (19) and (20) corresponding to strains of 1.9% and 7.8% respectively.
As the poly-crystalline sample is pulled the total free energy was monitored and used to calculate the stress, i.e., stress ≡ d 2 F/dζ 2 , where ζ is the relative change in the width of the crystal. Stress-strain curves are shown in Fig. (21) as a function of grain size and strain rate. In all cases the stress is initially a linear of function of strain until plastic deformation occurs and the slope of the stress-strain curve decreases. In Fig. (21a) the influence of strain rate is examined for the initial configuration shown in Fig. (18). It is clear from this figure that strain rate plays an strong role in determining the maximum stress that a sample can reach, or the yield stress, as has been observed in experiments [55]. The yield strength increases as the strain rate increases as would be expected. The influence of grain size on the stress-strain relationship is shown in Fig. (21b) for four grain sizes. The initial slope of the stress-strain curve (which will be denoted Y 0 in what follows) increases with increasing grain size as does the maximum stress, or yield stress, sustained by the sample. The yield strength and elastic moduli (Y 0 ) are plotted as a function of inverse grain size in Figs. (22) and (23) respectively for several strain rates. For each strain rate the yield stress is seen to be inversely proportional to the square root of the average grain size, however the constant of proportionality decreases with decreasing strain rate. Thus the PFC approach is able to reproduce the inverse Hall-Petch effect or the softening of nano-crystalline materials.
It would be interesting to observe the cross-over to the normal Hall-Petch effect where the yield stress decreases with increasing grain size. However, it is important to note that the initial conditions in these simulations was set up to explicitly remove the Hall-Petch mechanism, i.e., each nano-crystal was defect free. In addition thermal fluctuations were not included in the simulations. Nevertheless it is unclear whether or not a crossover may occur, due to the fact that low angle grain boundaries may act as sources of movable dislocations. Further study of this interesting phenomena for larger grain sizes would be of great interest.

D. Other Phenomena
There are many phenomena that the PFC method can be used to explore. To illustrate this a few small simulations were conducted to examine a number of interesting phenomena of current interest. In the next few sections some preliminary results are shown for grain growth, crack propagation and reconstructive phase transitions.

Grain Growth
When a liquid is supercooled just below the melting temperature small crystallites can nucleate homogeneously or heterogeneously. The crystallites will grow and impinge on neighboring crystallites forming grain boundaries. Depending on the temperature and average concentration the final state (i.e., in the infinite time limit) may be a single crystal or a coexistence of liquid and crystal phase since there exists a miscibility gap in density for some regions of the phase diagrams. For deep temperatures quenches the liquid is unstable to the formation of a solid phase and initially an amorphous sample is created very rapidly which will evolve into a poly-crystalline sample and eventually become a single crystal (in the infinite time limit). All these phenomena can be studied with the simple PFC model considered in this paper.
In this section the PFC model is used to examine the heterogenous nucleation of a poly-crystalline sample from a supercooled liquid state. A simulation containing fifty initial seeds (or nucleation sites) was conducted. The initial seeds were identical to those described in Sec. (II C 3) as were all other relevant parameters. The results of the simulations are shown in Fig. (24). Comparison of Figs. (24b) and (24c) shows that there is a wide distribution of grain boundaries each with a different density of dislocations (which appear as black dots in the figure). Comparison of (24c) with later configurations indicates that the low angle grain boundaries disappear much more rapidly than the large angle ones. The simple reason is that it is easy for one or two dislocations to glide in such a manner as to reduce the overall energy (this is usually accompanied with some grain rotation). The simulation was run for up to a time of t = 50, 000 (or approximately 1,200 diffusion times) and contained approximately 15,000 particles. The simulation took roughly 70 hours of cpu on a single alpha chip processor (xp1000).

Crack Propagation
The PFC model can be used to study the propagation of a crack in ductile (but not brittle) material. To illustrate this phenomena a preliminary simulation was conducted on a periodic system of size (4096∆x, 1024∆x) for the parameters (r, ψ o , ∆x, ∆t) = (−1.0, 0.49, π/3, .05). Initially a defect free crystal was set up in the simulation cell that had no strain in the x direction and a 10% strain in the y direction. A notch of size 20∆x × 10∆x was cut out of the center of the simulation cell and replaced with a coexisting liquid (ψ = 0.79). The notch provides a nucleating cite for a crack to start propagating. A sample simulation is shown in Fig. (25).

Reconstructive Phase Transitions
The simple PFC model can be used to study a phase transition from a state with square symmetry to one with triangular symmetry. In the model described by Eq. (26) a state with square symmetry is metastable, i.e., a state with square symmetry will remain unchanged unless boundary conduction or fluctuations are present. Bound- ary conduction or fluctuations allow for the nucleation of a lower energy state which in this particular model is the state of triangular symmetry discussed in Sec. (II C 1). A small simulation was performed to illustrate this phenomena. In this simulation a crystallite with square symmetry coexisting with a liquid was created as an initial condition. The parameters for this simulation were (r, ψ liq , ψ sol , ∆x, ∆t) = (1.0, 0.68, 0.52, 1.0, .02). The simulations depicted in Fig. (26), show the spontaneous transition from square structure to the triangular structure. Two variants of the triangular structure (differing by a rotation of 30 degrees) form in the new phase as highlighted in Fig. (26d).
A better method for studying this phenomena is to create a free energy that contains both square and triangular symmetry equilibrium states. This can be done by including a | ∇ψ| 4 term (which favors square symmetry) in the free energy. This is, unfortunately not the most convenient term for numerical simulations. A better approach is to simply couple two fields in the appropriate manner as was done in an earlier publication [1]. In either case an initial poly-crystalline state can be created of one crystal symmetry.

IV. SUMMARY
The purpose of this was paper was to introduce the PFC method of studying non-equilibrium phenomena involving elastic and plastic deformations and then to show how the technique can be applied to many phenomena. Those phenomena included epitaxial growth, material hardness, grain growth, reconstructive phase transitions, crack propagation, and spinodal decomposition. In the future, we intend to extend this model to study these phenomena in three dimensions. VI. APPENDIX: NUMERICAL METHODS Equation (28) was numerically solved using the two different methods work as described below. In what follows the subscripts n, i and j are integers that corresponds to the number of time steps and distance along the x and y directions of the the lattice respectively. Time and space units are recovered by the simple relations, t = n∆t, x = i∆x and y = j∆x.

A. Method I
In method I an Euler discretization scheme was used for the time derivative and the 'spherical laplacian' approximation was used to calculate all Laplacians. For this method the discrete dynamics read, ψ n+1,i,j = ψ n,i,j + ∇ 2 µ n,i,j , where µ n,i,j is the chemical potential given by; µ n,i,j = (r + (1 + ∇ 2 ) 2 )ψ n,i,j − ψ 3 n,i,j .

B. Method II
In method II an Euler algorithm was again used for the time step, except that a simplifying assumption was made to evaluate (r + (1 + ∇ 2 ) 2 )ψ n,i,j ) in Fourier space. In this approach the fourier transform of ψ n,i,j was numerically calculated then multiplied by w(q) and then an inverse fourier transform was numerically evaluated to obtain an approximation to (r + (1 + ∇ 2 ) 2 )ψ n,i,j ). If w(q) is chosen to be w(q) = r + (1 − q 2 ) 2 then, to within numerical accuracy, there is no approximation. In this work w(q) was chosen to be r + (1 − q 2 ) 2 if w(q) < −2.5 and w(q) = −2.5 otherwise. Thus w(q) is identical to the exact result for wavevectors close to q = 1, i.e., the wavelengths of interest. The advantage of introducing a large wavevector cutoff is that the most numerically unstable modes arise from the largest negative values of w(q). This allows the use of much larger time steps. Other than this approximation the method is identical to Method I.