# The displacement field associated with the freezing of a melt and its role in determining crystal growth kinetics

See allHide authors and affiliations

Edited by Pablo G. Debenedetti, Princeton University, Princeton, NJ, and approved January 4, 2020 (received for review September 11, 2019)

## Significance

We demonstrate that an accurate estimation of the displacements associated with the transformation of liquid into crystal is necessary to explain the striking variations in the temperature dependence of the addition rate of liquid atoms to the growing crystal interface during freezing. An assignment algorithm, adapted from operations theory, is shown to provide a good estimate of these atomic displacements. As the assignment algorithm requires only initial and target locations, it is applicable to all forms of structural transformations. In resolving a fundamental feature of the kinetics of freezing, a phenomenon of central importance to material fabrication, this paper also provides the tools to open lines of research into the kinetics of structural transformation.

## Abstract

The atomic displacements associated with the freezing of metals and salts are calculated by treating crystal growth as an assignment problem through the use of an optimal transport algorithm. Converting these displacements into timescales based on the dynamics of the bulk liquid, we show that we can predict the activation energy for crystal growth rates, including activation energies significantly smaller than those for atomic diffusion in the liquid. The exception to this success, pure metals that freeze into face-centered cubic crystals with little to no activation energy, are discussed. The atomic displacements generated by the assignment algorithm allows us to quantify the key roles of crystal structure and liquid caging length in determining the temperature dependence of crystal growth kinetics.

Crystal growth from the melt is a commonplace phenomenon, one fundamental to material fabrication, whose kinetics depends on extraordinary sequences of cooperative events. In this paper we shall demonstrate how a map of this atomic choreography, the atomic displacement field, can be obtained from the structure of the interface alone and show how this displacement field allows us a conceptually deeper and, quantitatively, more accurate treatment of crystal growth kinetics. While the use of a displacement field that maps particles from liquid to crystal is novel, such maps are standard in the description of solid–solid transformations (1) where the existence of well-defined initial and final structures makes defining the associated transformation displacements a simple task. If the required displacements are less than a crystal layer spacing, the transition is classified as diffusionless or displacive*,* otherwise it is described as diffusive (1). In freezing and melting transitions, the large excursions of atoms in a liquid obscure those movements essential for the transformation into the crystal phase. For this reason, freezing has not been previously analyzed in terms of an atomic displacement field. The goal of this paper is to present an explicit definition and determination of the displacement field for the freezing of a liquid and to explain how to integrate this information into the theory of growth kinetics.

The “classical” theory of crystal growth, as expressed in the theory of Wilson (2) and Frenkel (3) and the solid-on-solid model (4), considers crystallization to take place via the incoherent deposition of particles from an otherwise unstructured liquid into available sites on the growing crystal surface. The timescale for this deposition is assumed to be set by a liquid transport coefficient, diffusion (3), or viscosity (2). This last assumption corresponds to an implicit proposition that crystal growth is diffusive in nature. This classical picture may be appropriate for precipitation from a dilute solution (5) but it ignores, in the context of growth from the melt, the strong local correlations that govern the structure and dynamics of the liquid state and the diffuse structure of the crystal–liquid interface (6, 7). These omissions are addressed, to a degree, in a Landau–Ginzburg treatment of freezing (8⇓–10). Eschewing reference to individual particles, this approach treats the crystal as a liquid with periodic density waves and the dynamics of crystal growth is now governed by the same timescale that governs density fluctuations in the liquid at the Bragg wave vector (8). While the density wave representation makes sense for the construction of an equilibrium description of the freezing transition, it is not clear that the dynamics of ordering can be equally well described without reference to the motions of individual atoms. Our goal in this paper is to explore an alternative account of the process of freezing based on an extension of the existing formalism for solid–solid transformation.

Previously (11), having demonstrated that the crystal growth of a number of metals that form face-centered cubic (FCC) crystals is nondiffusive, we argued that the absence of activated control arose as a consequence of the crystalline order that characterized the local ground state of the liquid at the crystal–liquid interface. The evidence of the existence of these ordered interfacial ground states is the propagation of the crystal interface during energy minimization. An example of this behavior is shown in Fig. 1 along with an example of a system for which there is no advance of the interface. These interfacial ordered ground states are in sharp contrast to those of the bulk liquid that exhibit no such order.

In a second publication (12), the mapping from liquid to crystal associated with energy minimization into the ordered interfacial ground states was used to provide the displacement maps for the FCC metals and a molten salt, NaCl. An example of the resulting distribution of displacements for Cu is shown in Fig. 2*A*. The time required for atoms in the bulk liquid to move the median value of the displacement field (Fig. 2*B*) was found to provide a reasonable prediction of the observed activation energy for crystal growth.

The existence of ordered interfacial ground states––the observation that formed the basis of the analysis in refs. 11 and 12—is not the norm in crystal growth. In ref. 12 Fe and ZnS were shown to not order at the interface during energy minimization. The (111) face of a Lennard-Jones crystal is another example (11). We propose that the majority of materials do not exhibit ordering on minimization. Clearly, we cannot simulate every system. Our claim is based on the generic experimental observations of activated control of crystal growth (13) for materials other than pure metals and some simple salts. Activation indicates an energy barrier for the ordering process at the interface and, hence, strongly implies the inaccessibility of the crystal from the liquid by any minimization algorithm that imposes a continuous decrease in the potential energy.

## Transformation Displacement Fields as a Problem of Optimal Transport

If energy minimization does not provide a mapping from liquid to crystal at the interface, we must find another way of defining the transformation displacement field for an arbitrary crystal growth problem. In general, the displacement field corresponds to an assignment mapping where each liquid atom is assigned to a crystal site. The total number of distinct possible ways of assigning *N*_{L} liquid atoms to *N*_{C} crystal sites (both regarded as distinguishable by virtue of their locations) is *N*_{C} *> N*_{L} (14). The displacement fields for most of these possible assignments will be of no physical significance. We propose that the set of displacements of interest is that with the minimum average displacement. While this minimal field may not represent the path taken by any given ordering trajectory, it does provide a lower bound on the characteristic transformation displacement. This means that if this minimal length exceeds the threshold length obtained from the liquid cage length, then we can be confident that the associated freezing transition is diffusional in character.

The problem of assigning atoms to crystal sites to minimize the mean displacement is an example of a problem well known in linear programming where it is called the optimal assignment problem or, more generally, the optimal transport problem (15). Transport problems occur frequently in operations research––matching workers to tasks, examinations to time slots, etc. A widely used solution is the Hungarian Algorithm (HA) described by Kuhn in 1955 (16). We have implemented the HA for the crystal growth problem as follows. 1) All calculations have been carried out with a particle-to-site ratio, *N*_{L}*/N*_{C} *∼0.6*. 2) The reported displacements are the components in the interfacial plane only. 3) The assignment is applied to atoms in the liquid that lies within two lattice spacings of the interface position. (A detailed explanation of these choices in the implementation of the HA is provided in *SI Appendix*.) An example of the displacement field for an interfacial plane of the Cu crystal–liquid interface is shown in Fig. 3. We note the presence of collective motions arising solely from the particle exclusion associated with the requirement of single occupancy of crystal sites.

## Assignment Displacements and the Temperature Dependence of Crystal Growth Rates

The intrinsic crystal growth rate is the rate at which the crystal front propagates when unrestrained by heat or concentration diffusion. It is the growth rate controlled by the kinetics of microscopic ordering alone. The temperature dependence of these rates allows us to empirically distinguish diffusive from displacive transformations based on whether the associated activation energy is equal to or less than that for diffusion in the bulk liquid. We shall apply the standard resolution (17) of the growth rate *v*(*T*) into a thermodynamic term and a kinetic coefficient *k*(*T*) using*μ* is the chemical potential difference and *k*_{B} is the Boltzmann constant. Estimating the chemical potential difference as *T*_{m} being the melting point and *h* being the enthalpy per particle of the respective phase), Eq. **1** allows us to extract the kinetic coefficient from the calculated rates *v*(*T*). In this paper we shall analyze the crystal growth kinetics of the following model systems. As examples of a metal and a salt whose crystal growth shows little to no activation barrier, we include Cu and NaCl. As examples of a metal and a salt whose crystal growth shows significant activation (and no sign of an ordered interfacial ground state), we include Fe and ZnS. Finally, we include a family of binary Lennard-Jones equimolar mixtures that freeze into a cubic CsCl crystal. This mixture is based on a model due to Kob and Andersen (KA) (18) that has been studied extensively in the context of glass formation. These mixtures have been included as their crystal growth kinetics provide a continuous range of activation energies through a small adjustment of their interaction (as described below). The metals are modeled using an embedded atom model potential due to Foiles et al. (19) and Mendelev et al. (20), NaCl is modeled using the Tosi–Fumi potential (21), and ZnS is modeled using a potential due to Grünwald et al. (22). The binary Lennard-Jones mixtures A_{50}B_{50} are modeled with the following interaction potential:*i* and *j* and finite distance at which the potential (first) vanishes, respectively, *v*(*T*) were calculated using molecular-dynamics simulations from the LAMMPS (23) set of algorithms. As described previously (11, 12), a uniform temperature is maintained throughout the run and the density is allowed to adjust through the inclusion of free liquid. These details, along with the crystal growth rates *v*(*T*), are provided in *SI Appendix*. For FCC, we have calculated growth along the <111> direction, and for body-centered cubic (BCC) structure, the growth direction investigated is along the <110> direction. The growth directions of the rocksalt crystal for NaCl and the wurtzite for ZnS are both along the <100> direction. In the case of the CsCl crystal, we have considered growth along the <100> direction. Details of the simulation protocol along with the tabulated *v*(*T*) vs. *T* data for all of the models are provided in *SI Appendix*.

In Table 1, we present the activation energy *E*_{a}, obtained empirically by fitting the rate coefficient *k*(*T*) [obtained from *v*(*T*) using Eq. **1**)] to an Arrhenius function, *k*(*T*_{m}) is the kinetics rate and *E*_{a} is the energy activation energy. For comparison, we also present the activation energy *E*_{D} for self-diffusion in the homogeneous liquid, obtained by an Arrhenius fit to the self-diffusion coefficient *D* of the form *d* is translated into a time *τ* using the relation *P*(*τ*) we calculate the average rate *SI Appendix*). The resulting values for

We find, as shown in Fig. 5, that the estimate of the activation energy *E*_{a}/*E*_{D} and, hence, the degree of activated control of the crystal growth kinetics. A clear exception to this success is Cu where the assignment analysis results in a significant overestimate in the barrier to growth. It is likely that this represents a general limitation of the assigned displacements. When *E*_{a} <<1.0, the dynamics of growth is dominated by collective vibrations about the crystalline ground state in the interface which is not adequately described by simply assigning atoms to sites based on minimizing the mean displacement. A more appropriate approach for this limiting behavior is to use the characteristic vibrational frequency of the bulk crystal as proposed in ref. 11.

The general success of the freezing displacement fields, as generated by the assignment analysis, in predicting the degree of activated control of crystal growth kinetics represents the major result of this paper. In using the liquid MSD data to map lengths to timescales we are following the same reasoning used in the earlier theories of the kinetic coefficient for crystal growth. Where we deviate from these earlier approaches is in providing a more nuanced measure of the displacement lengths associated with ordering a liquid. The Wilson–Frenkel theory of crystal growth assumes that *E*_{a}/*E*_{D} = 1.0. The marked deviations from this “classical” theory in Fig. 5 are the direct result of identifying length scales small enough to access dynamics that show little activated control. While the lengths we find using the optimal transport analysis may differ by only a small amount from those assumed in the earlier studies (11, 12), the difference can be physically significant and give rise to very large differences in both the temperature dependence of crystal growth rates as well as its absolute magnitude.

## Mapping to Alternate Crystal Structures: Polymorphs and Substitutionally Disordered Crystals

We have reported that Cu freezes into an FCC crystal by a diffusionless process but the freezing of Fe to a BCC structure is diffusional. What is the origin of the difference? Is it some property of the liquid or does the difference lie with the different crystal structures? The optimal transport analysis provides a useful means of directly addressing question such as this since it allows us to select any target crystal structure we choose. If Fe was to freeze into the FCC structure, as opposed to BCC, how would the distribution of displacements and, hence, the nature of the transformation change? The assignment analysis allows us to address this question since we are free to select the target crystal to map the liquid particles to when generating the displacement distribution. The value of

Referring to Table 2, we find that NaCl does exhibit a significantly smaller value of *E*_{a}*/E*_{D} > 0.4, does not depend on the crystal structure but is a consequence of the liquid. It appears that when the barrier to growth is particularly small (as in the case of NaCl), the choice of crystals may play a role in influencing the temperature dependence of growth. The results for the random crystals in Table 2 showed little variation for the values of

## Conclusions

We have demonstrated that crystallization from the melt can occur via crystal growth with an activation energy that can take on a value

The exception to the success of the assignment analysis are those materials that exhibit very small or no activation barrier to crystal growth. The pure FCC-forming metals belong to this general group. As shown in Fig. 5 and Table 1, the assignment analysis significantly overestimates *E*_{a} for Cu. Previously (11), we have argued that, in the absence of a barrier, the process of crystallization more closely resembles a complex collective phonon mode than the incoherent motion that characterizes liquid dynamics. It is not surprising, then, that the simple minimization condition that defines the displacements generated by the HA would be a poor estimate of the coherent motions involved in ordering when the local potential energy minimum is crystalline. For these materials, the characteristic frequency of the crystal vibrations provides a better estimate of the magnitude of the ordering coefficient *k*.

Previous studies have observed anisotropy in crystal growth rates. In the case of the Lennard-Jones (LJ) FCC crystal, the (111) surface exhibits activated growth while the (100) surface grows without any apparent barrier (25, 26). Applying the assignment analysis to these two surfaces, we find only a small decrease in the predicted activation energy for growth estimated for LJ(100) surface relative to that the LJ(111) surface, 1.6*ε* and 1.8*ε*, respectively. This result for LJ(100) surface is consistent with the failure of the assignment analysis, already noted, when applied to barrierless growth such as that of Cu. We leave a more extensive study of the origin of anisotropy in crystal growth kinetics for future study.

Our introduction of transformation displacements and the HA for their calculation opens up lines of study of crystallization kinetics based on measures of configurational proximity (i.e., the magnitude of the transformation displacements) between a disordered and ordered structure. As demonstrated here, a researcher has complete freedom in the choice of the target structure. Polymorphs, substitutional disorder, or orientational disorder (in the case of molecules) could be considered in exploring how different aspects of structure influence the characteristic displacements atoms or molecules must undergo in a phase transformation. The optimal transport analysis can also be extended to transformation from one disordered structure to another in bulk liquids, quantifying the kinetic accessibility of different distinct configurations in a liquid. A number of these questions are the subject of current research.

### Data Availability Statement.

All the data supporting the findings of this study are available within this paper and the *SI Appendix*.

## Acknowledgments

P.H. would like to gratefully acknowledge a helpful conversation with Toby Hudson concerning the assignment problem and Paul Madden for discussions on solid-state transformations. This work was supported by a Discovery Grant from the Australian Research Council.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: peter.harrowell{at}sydney.edu.au.

Author contributions: P.H. designed research; G.S. and A.H. performed research; G.S. and P.H. analyzed data; and G.S. and P.H. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1915806117/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- D. A. Porter,
- K. E. Easterling,
- M. Y. Sharif

- ↵
- H. A. Wilson

- ↵
- J. Frenkel

- ↵
- G. H. Gilmer,
- H. J. Leamy,
- K. A. Jackson

- ↵
- A. Pimpinelli

- ↵
- ↵
- J. Zhou et al

- ↵
- L. V. Mikheev,
- A. A. Chernov

- ↵
- K. A. Wu,
- C. H. Wang,
- J. J. Hoyt,
- A. Karma

- ↵
- K. A. Wu,
- S. C. Lin,
- A. Karma

- ↵
- G. Sun,
- J. Xu,
- P. Harrowell

- ↵
- A. Hawken,
- G. Sun,
- P. Harrowell

- ↵
- ↵
- G. E. Martin

- ↵
- C. Villani

- ↵
- ↵
- ↵
- ↵
- ↵
- M. I. Mendelev,
- S. Han,
- D. J. Ackland,
- D. Y. Sun,
- M. Asta

- ↵
- F. Fumi,
- M. Tosi

- ↵
- ↵
- S. J. Plimpton

- ↵
- J. Timmermans

- ↵
- E. Burke,
- J. Q. Broughton,
- G. H. Gilmer

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Chemistry