# Phase-field-crystal model for fcc ordering

## Abstract

We develop and analyze a two-mode phase-field-crystal model to describe fcc ordering. The model is formulated by coupling two different sets of crystal density waves corresponding to and reciprocal lattice vectors, which are chosen to form triads so as to produce a simple free-energy landscape with coexistence of crystal and liquid phases. The feasibility of the approach is demonstrated with numerical examples of polycrystalline and (111) twin growth. We use a two-mode amplitude expansion to characterize analytically the free-energy landscape of the model, identifying parameter ranges where fcc is stable or metastable with respect to bcc. In addition, we derive analytical expressions for the elastic constants for both fcc and bcc. Those expressions show that a non-vanishing amplitude of [200] density waves is essential to obtain mechanically stable fcc crystals with a non-vanishing tetragonal shear modulus . We determine the model parameters for specific materials by fitting the peak liquid structure factor properties and solid density wave amplitudes following the approach developed for bcc [K.-A. Wu and A. Karma, Phys. Rev. B 76, 184107 (2007)]. This procedure yields reasonable predictions of elastic constants for both bcc Fe and fcc Ni using input parameters from molecular dynamics simulations. The application of the model to two-dimensional square lattices is also briefly examined.

###### pacs:

61.72.Mm,68.08.De,68.08.-p,81.16.Rf## I Introduction and summary

The phase-field-crystal (PFC) method has emerged as an attractive computational approach to simulate the evolution of crystalline patterns (1); (2); (3); (4); (5); (6). By resolving the crystal density field, it naturally incorporates defects and elastic interactions arising from localized and large scale distortions of this field, respectively. Moreover, this method can in principle be used to simulate microstructural evolution on diffusive time scales that are much longer than typical time scales accessible by molecular dynamics (MD) simulations.

Like classical density function theory (DFT), the PFC method is based on representing the free-energy of a material by a functional of its density (7); (8); (9); (11); (10); (12); (13). However, classical DFT and the PFC method use different functionals to achieve different goals. Classical DFT seeks a physically realistic mean-field description of the crystal density field to reproduce quantitatively as accurately as possible the properties of a material. Since is sharply peaked around mean atomic positions, this generally requires a very large number of terms in the traditional expansion of the number density as a sum of density waves

(1) |

where each represents a different reciprocal lattice vector (RLV) in this unrestricted sum. In contrast, by using a considerably simplified density functional, the PFC method essentially restricts this sum to a much smaller set of reciprocal lattice vectors in order to simulate efficiently the evolution of the crystal field on length and time scales as large as possible. Recent studies have shown that, despite this loss of realism, PFC models are able to reproduce quantitatively certain key properties that influence microstructural evolution such as the crystal-melt interfacial free-energy (14); (15), the bulk modulus (15), and grain-boundary energies (15), which have been computed for the test case of pure Fe.

Despite this progress, the PFC method has only been developed for a small set of crystal structures. The original formulation of Elder et al. (1); (2) uses the same free-energy functional as the Swift-Hohenberg model of pattern formation (16); (17) of the form

(2) |

with the free-energy density

(3) |

where represents the crystal density field. This one-mode model essentially truncates the sum (1) to one set of RLVs with equal magnitude since higher modes have much smaller amplitude. As a result, it favors crystal structures for which the principal RLVs can form “triads” (i.e. closed triangles), which include hexagonal and body-centered-cubic (bcc) ordering in two and three dimensions, respectively. Aside from favoring those structures, triad interactions are essential for solid-liquid coexistence. This is because in a weakly nonlinear expansion of the bulk free-energy density of the form, (with for all principal RLVs), triads contribute a cubic term with a negative coefficient . Since and are both positive, this cubic term is responsible for the existence of a free-energy barrier between the two minima of corresponding to liquid () and solid ().

In this paper, we use a “two-mode” phase-field-crystal model to model face-centered-cubic (fcc) structures, which has the free-energy density

(4) |

This model truncates the sum (1) to two sets of RLVs with magnitude and , respectively, where the first set corresponds in general to the principal RLVs and the second to some other set with larger wavevector magnitude; all other RLVs have much smaller amplitude. This construct provides more flexibility to form triad interactions by combining RLVs from those two sets, and hence to describe other crystal structures. We demonstrate this here for fcc ordering, which is obtained by choosing the sets and to correspond to (principal set) and RLVs, respectively, with . While there is in principle freedom in the choice of the second set for a given structure, we have chosen this set such that is as small as possible, as desired for computational efficiency.

The form (4) reduces in the limit to the free-energy density introduced by Lifshitz and Petrich (18) as a generalization of the Swift-Hohenberg model to describe two-dimensional quasiperiodic patterns observed in Faraday wave experiments, which result from the superposition of two frequencies. Although formulated primarily to describe those patterns, this model was also shown to describe other patterns, including regular square crystal lattices in two dimensions with the choice , which couples and RLVs.

The present introduction of the parameter in the form (4) provides the additional flexibility to change the relative stability of different crystal structures. This is because in the limit , this form reduces formally to the original Swift-Hohenberg form (3) after a simple rescaling of the parameters. Hence, as is increased the contribution of the second -mode becomes less significant in comparison to the first -mode. Consequently, as is increased from zero, the crystal structure favored by the two-mode interaction becomes metastable with respect to the one-mode structure. This added capability to model the coexistence of two different crystal structures, in addition to the coexistence of each structure with a liquid, should prove useful to model a wide range of phase transformations with a PFC approach.

In the next section, we scale the parameters of the model to write the free-energy functional in a dimensionless form with only three parameters: , which is the standard PFC model parameter analogous to temperature that controls the size of the solid-liquid coexistence regions as a function of density, , whose value is generally determined by the choice of crystal structure, and controls the relative stability of the two-mode and one-mode structures (fcc and bcc, respectively). In this section, we also use a standard common tangent construction to compute the phase-diagram in the plane of density and for an illustrative choice of . The phase diagram exhibits regions of bcc-liquid and fcc-liquid coexistence for small and large epsilon, respectively. The size of the fcc-liquid coexistence region depends generally on . For where Eq. (4) reduces to the free-energy density of Lifshitz and Petrich (18), the analog phase-diagram only exhibits fcc-liquid coexistence, so that a finite is necessary for the phase diagram to exhibit both bcc-liquid and fcc-liquid coexistence. We demonstrate the feasibility of the approach with some simulations of polycrystalline growth and (111) twin growth. A numerical computation of the (111) twin boundary energy for parameters of Ni is given in an appendix. The ability to model twin growth is important for solidification modeling since twins can dramatically alter both eutectic (19); (20) and dendritic (21) microstructures.

In section III, we carry out an amplitude expansion of the bulk free-energy density in the small limit. This expansion exploits the property that, with the scaling , the amplitudes of the and density waves scale as, and , respectively, while the density difference between solid and liquid scales . Therefore, this density difference can be neglected in the small limit and the bulk free-energy density can be expressed solely in terms of those amplitudes. As required for solid-liquid coexistence, the free-energy density has minima in the (,) plane corresponding to liquid () and fcc solid (finite and values that depend on ). By comparing this form to the free-energy density for a single amplitude of bcc density waves (corresponding to RLVs), we identify different regions of relative fcc and bcc stability, which explains the phase diagram computed in section II.

In section IV, we discuss how to determine the two-mode PFC model parameters to relate them quantitatively to different materials. We follow essentially the same approach developed by two of the authors for the standard PFC one-mode model for bcc ordering (14). For bcc, the parameters were completely determined by fitting three parameters: (i) the peak value of the liquid structure factor, , where , (ii) the second derivative of the fourier transform of the direct correlation function at this peak, , and (iii) the solid density wave amplitude . For fcc, all the parameters except are determined by the same fit, where . (The shape of the structure factor at is not realistically modeled given the limited number of model parameters.) then determines the ratio of the and solid amplitudes, which can be varied to alter the relative stability of fcc and bcc.

In section V, we derive analytical expressions for the three independent elastic constants of a cubic material, , , and , for both the standard one-mode PFC model (3) and the present two-mode model (4). We use a brute force approach that consists of calculating to quadratic order the change of solid free-energy density, modeled by a one- or two-mode approximation for bcc and fcc, respectively, due to small dilation or shear transformations of the unit cell. We have checked that we obtain identical expressions to those derived recently by Spatschek and Karma for general lattices using an amplitude equation framework (22), which provides a non-trivial self-consistent test of our calculations. For the one-mode bcc model (3), the elastic constants are

(5) |

where . For the two-mode fcc model (4),

(6) |

and

(7) |

where and for simplicity.

Using values of and density wave amplitudes from molecular dynamics simulations for parameters of bcc Fe and fcc Ni, we find that the above expressions give reasonable estimates of elastic constants (e.g., GPa for one-mode bcc PFC model compared to GPa in MD Fe and GPa for the two-mode fcc PFC model compared to GPa in MD Ni). The predicted values generally tend to be lower than the constants computed from MD simulations, but such discrepancies are to be expected given the PFC models are based on one or two modes.

The analytical predictions for the elastic constants allow us to draw two important general conclusions pertaining to the development of PFC models for different crystal structures and to the method used to determine the parameters of those models.

The first conclusion, which follows directly from Eqs. (6) and (7), is that the presence of the second mode, which corresponds to [200] density waves, is essential to obtain a physically meaningful set of elastic constants for fcc. Without this second mode (), Eqs. (6) and (7), predict that . This implies that the tetragonal shear modulus vanishes, and that the system is mechanically marginally stable. Of course, these analytical expressions for the elastic constants neglects the contributions of higher modes that are present in a full solution of the PFC equations. However, those higher modes are generally small for the small values of corresponding to Fe and Ni parameters. Therefore, the contributions of those modes will generally be small and will not change qualitatively this picture.

While it is in principle possible to select energetically different crystal structures in the PFC model with the addition of other nonlinearities in the free-energy density (such as and ) (23), this approach will be of limited applicability for crystal structures like fcc where one mode does not suffice to produce the correct elastic properties. This is also true for simple cubic lattices and two-dimensional square lattices. The latter are briefly examined in section VI by coupling and density waves.

The second conclusion, which is general, is that the elastic constants are uniquely determined once the phase-field model parameters have been fitted to the peak liquid structure properties, which fixes , and the solid density wave amplitudes, as in the approach of Wu and Karma (14) summarized above. This also fixes the value of the elastic bulk modulus

(8) |

In general, the bulk modulus can also be defined from the thermodynamic relation

(9) |

where is the total free-energy, is the volume, and is the number density. The second equality in the last equation can in principle be used to compute the bulk modulus directly from the PFC solid free-energy curve ( versus ), without computing the elastic constants. For a perfect crystal without vacancy, Eqs. (8) and (9) should in principle predict the same bulk modulus. However, the two definitions can give different predictions for the PFC model because the number of atoms per peak of the crystal density field is not constrained to unity. While the average number of atoms per peak will also differ from unity in a real crystal with vacancies, thereby altering the open-system elastic constants (24), the vacancy concentration is generally very small even at melting. How to meaningfully relate the predictions of Eqs. (8) and (9) for the bulk modulus is unclear in the PFC approach that, by construct, does not use a realistic description of the crystal density field, and also does not model vacancy formation explicitly.

Despite these limitations of the PFC approach, Jaatinen et al. (15) have recently proposed a modified one-mode PFC model to remedy the fact that, for the standard one-mode PFC model with the free-energy density (3), the bulk modulus predicted by Eq. (9) is several times smaller than the experimental value for parameters of bcc Fe. Their model yields a value of the bulk modulus computed through Eq. (9) that is in better agreement with experiment and also gives an improved prediction of the density difference between solid and liquid. It gives similar predictions of crystal-melt interfacial free-energies for bcc Fe as obtained previously by Wu and Karma using the standard one-mode model (14).

In the light of Eq. (5), it is apparent that any one-mode model that fits the correct peak structure factor properties and solid density wave amplitudes should predict the same elastic constants. This is consistent with the fact that Eq. (5) predicts a shear modulus GPa for the standard one-mode model of bcc Fe, which is reasonably close to the value GPa estimated by Jaatinen et al. (15) from numerical shearing experiments in their model for similar input parameters.

Since elastic constants are a major determinant of grain boundary energies and long-range interactions between crystal defects, reproducing those constants, and hence the bulk modulus predicted by Eq. (8), appears essential for modeling microstructural evolution. Also requiring that Eq. (9) predicts the correct bulk modulus using the solid free-energy curve may appear desirable. However, the motivation for doing so in the context of simple PFC models is somewhat less clear given the lack of realism of the crystal density field and the fact that Eqs. (5)-(7) predict reasonable values of the elastic constants. In fact, any one- or two-mode model with the same peak liquid structure factor properties and density wave amplitudes will predict essentially the same elastic constants associated with the free-energy cost of lattice distortions, and also the same interfacial energies as can be inferred from amplitude equations (14). Since those elastic constants and interfacial energies are the quantities that matter most for modeling microstructural evolution in a PFC context, we have not found it necessary to formulate the two-mode PFC model in such a way that the bulk modulus is also correctly predicted from the solid free-energy curve using Eq. (9). Accordingly, we follow essentially the same approach outlined in Ref. (14) for determining the PFC model parameters.

## Ii Phase-field crystal model

### ii.1 Basic equations and scalings

The PFC equations have the standard form for conserved dynamics

(10) |

where is the free-energy functional defined by Eq. (2) with the free-energy densities given by Eqs. (3) and (4) for the one- and two-mode models, respectively. To minimize the number of parameters, it is useful to rewrite the equations in dimensionless form. For the two-mode model, we define the dimensionless parameters

(11) |

(12) |

(13) |

where we set for fcc (i.e. equal to the ratio of the magnitudes of the and RLVs). We also define the dimensionless variables

(14) |

(15) |

(16) |

(17) |

Substituting the above definitions into Eqs. (2) and (4) yields the dimensionless form

(18) |

with the free-energy functional

(19) |

and free-energy density

(20) |

where we have dropped the prime symbol on the dimensionless spatial coordinate vector for brevity. Even though most of the paper focuses on the two-mode model, we also compute in section V the elastic constants for the standard one-mode PFC model. For this model, we use the same scaling as in Ref. (14) with the parameter

(21) |

and dimensionless variables

(22) |

(23) |

(24) |

where is defined by Eq. (14). Substituting the above forms into Eqs. (2) and (3) yields (after dropping the prime symbol on ) the dimensionless form of the one-mode PFC equations (18) and (19) with

(25) |

### ii.2 Phase diagram

The phase diagram of the two-mode PFC model is obtained by computing the free-energy density as a function of the mean density in solid and liquid, denoted by , and , respectively, and then using a standard common tangent construction to obtain equilibrium values of in solid () and liquid ().

Since the density is constant in the liquid, is obtained directly from Eq. (20)

(26) |

For small , the solid free-energy density can be well approximated by only considering the contribution of the and RLVs. Accordingly, the crystal density field is expanded in the form

(27) | |||||

where we have used the fact that all density waves have the same amplitude in the crystal ( and ) and the magnitude of the principal RLVs are unity in our dimensionless units so (). The parameters and are solved by substituting Eq. (27) into Eqs. (19) and (20) and by minimizing the resulting free-energy with respect to and . This minimization yields the solid free energy density

(28) | |||||

where and are themselves functions of . The coexistence densities and are computed numerically using the standard common tangent construction, which consists of equating the chemical potentials and grand potentials of the two phases. It is also necessary to compute the solid free-energy curve for bcc since the latter can have a lower free-energy than fcc for some regions of the phase diagram. The bcc free-energy density was obtained by expanding the crystal density field using a one-mode approximation, which only involves RLVs as in Ref. (14), and substituting this expansion into the two-mode model defined by Eqs. (19) and (20).

An example of the phase diagram for is shown in Fig. 1, where we also show for completeness the hexagonal and stripe phases. As desired, we obtain a large range of fcc-liquid coexistence. For small , however, bcc becomes favored over fcc. A common tangent construction using fcc and bcc free-energy curves shows that the density range of bcc-fcc coexistence is extremely narrow for small values of and cannot be resolved on the scale of Fig. 1. As will be explained later in section III.3, the range of where bcc is favored depends on the value of . In the limit , the two-mode model reduces to the standard one-mode model after a simple rescaling of parameters, which can be easily seen by comparing Eqs. (20) and (25). Hence, increasing reduces the contribution of the second mode. Conversely, reducing increases the contribution of this mode and tends to favor the fcc structure, which extends to smaller for smaller . In the extreme case where , the region of fcc-liquid coexistence extends all the way to vanishingly small as shown in Fig. 2.

### ii.3 Numerical examples

We now demonstrate the feasibility of the model with some numerical examples of fcc polycrystalline growth and (111) twin growth. The PFC conserved dynamics governed by Eq. (18) with the free-energy defined by Eqs. (19) and (20) was solved using the semi-implicit pseudo-spectral scheme given by Eq. (A2) in Appendix A of Ref. (25). We used the parameters and obtained from our fit of pure Ni presented later in section IV, together with the grid spacing , which determines the number of Fourier modes, and the time step . For this value of and , the computations presented in the next section show that the size of the solid-liquid coexistence region is extremely small, i.e. is two orders of magnitude smaller than as can already be seen from the phase diagram in Fig. 2, and .

The first example in Fig. 3 shows the growth of small fcc crystallites of different orientations for a value of that is well inside the stable fcc-solid region of the phase diagram. The crystallites grow as expected until they collide to form grain boundaries. The second example in Fig. 4 shows a (111) twin crystal for a value of at coexistence and for a system size chosen such that a twin crystal with two stacking faults fits perfectly the periodic boundary conditions in all directions without any liquid present. A computation of the excess free-energy of this twin boundary given in the appendix to this paper yields a value of approximately 30 mJ/m that falls within the range of values typically reported in the literature for fcc metals. Fig. 5 then shows the growth of the same twin crystal in a supercooled liquid for a much larger system with .

## Iii Amplitude equations

### iii.1 Scalings

In this section, we analyze in more detail the properties of the model by expanding the free-energy in terms of the amplitudes of density waves. In the one-mode bcc case analyzed in Ref. (14), a similar expansion exploited the fact that the amplitude of density waves scales as in the small limit. In the present case, the expansion is rendered more difficult by the presence of two different sets of density waves with amplitudes and corresponding to and RLVs, respectively. Therefore, it is not a priori obvious how and should scale in the small limit. If is kept constant, the bcc structure turns out to always be favored in the small limit as apparent in the phase diagram of Fig. 1. Consequently, a small amplitude expansion that captures the fcc structure cannot be carried out at fixed . However, if is decreased proportionally to by imposing the additional scaling , both and scale as , thereby making a rigorous expansion possible. This expansion may seem artificial since the phase diagram of Fig. 1 is computed at fixed . However, as we show below, the results of this expansion can be used to understand the small structure of the phase diagram, in particular the relative stability of fcc and bcc.

To demonstrate the feasibility of this expansion, we first analyze fcc-liquid coexistence for small with the scaling . The equilibrium densities are calculated using the common tangent construction described in the previous section. To make the dependence of the coexistence densities on explicit, we make a log-log plot of the mean coexistence density versus for three different values of . The results in Fig. 6 show that the mean coexistence density scales as . Next in Fig. 7, we show a log-log plot of the density difference between solid and liquid versus for the same three values of . The results show that . Together, these two log-log plots show that, in the small limit, the two-mode PFC model exhibits a weak first-order freezing transition where the size of the solid-liquid coexistence region is at the order of that is much smaller than the mean value of the density .

### iii.2 Free-energy functional

The above scalings suggest that we can expand the crystal density field in powers of as

(29) |

and expand accordingly the average densities

(30) |

and

(31) |

in the liquid and solid, respectively. The numerically determined scaling relations and then imply that

(32) |

and

(33) |

Next, to carry out the amplitude expansion, we start from the equilibrium equation , where is the equilibrium value of the chemical potential. With defined by Eqs. (19) and (20), we obtain

(34) |

We substitute the small expansion of the density field (29) into Eq. (34) and collect terms with the same power of . We find at the order

(35) |

which has the solution

(36) |

where the summations are over and RLVs, respectively, and , in our scaled units. At order , we obtain

(37) |

which has the solution

(38) |

and collecting the terms at order yields

(39) |

Since gives a vanishing contribution for all density waves associated with sets and , all remaining terms and must balance each other in order for a solution of Eq. (III.2) to exist. For example, the condition that the coefficients of balance each other yields

(40) |

and requiring that the coefficients of balance each other yields in turn

(41) |

The above solvability condition must be satisfied independently for each reciprocal lattice vector. This yields a set of fourteen coupled amplitude equations that are straightforward to obtain. From those amplitude equations, it is useful to express the free-energy of the system measured from its constant value in the liquid, defined as the difference , as a functional of the density wave amplitudes and . This quantity can be expressed solely in terms of the amplitudes of density waves owing to the property that the size of the coexistence region () is much smaller than the mean density () in the small limit. Since the amplitudes are not conserved order parameters, the equilibrium state simply corresponds to a minimum of this free-energy. Hence the amplitude equations must be equivalent to

(42) |

and

(43) |

For the case where the amplitudes are spatially uniform, we obtain the free-energy density

(44) | |||||

In general, the above free-energy density is a multi-variate function of the fourteen amplitudes and of different density waves, and cannot be represented graphically in a simple way. However, the free-energy barrier between solid and liquid can be made explicit by assuming that all the and all the density waves have the same amplitude (i.e., and , respectively). In this isotropic approximation (see (14) for the bcc analog), the free-energy density becomes

(45) |

This expression can also be obtained by evaluating directly the difference between the solid and liquid free-energy densities, , with and given by Eqs. (26) and (28), respectively, and the substitutions