Thermodynamic Phase Diagram

Fig. 1 shows the temperature–pressure phase diagram of the embedded atom method (EAM) model–Cu system. The melt line $$ T_m\left(P\right)$$ has been obtained from two-phase (solid–liquid) coexistence simulations in the isobaric–isoenthalpic (NPH) ensemble (464748–49). Details are described in SI Appendix. It is found that the equilibrium solid phase at $$ T_m$$ is fcc for pressures $$ P<71.6$$ GPa, hcp in the interval $$ 71.6<P<85$$ GPa, and bcc for $$ P>85$$ GPa. The phase diagram contains two triple points: one at 71.6 GPa and 3,320 K, where fcc, hcp, and liquid coexist, and one at 85 GPa and 3,598 K, where hcp, bcc, and liquid coexist. A third three-phase coexistence point exists at 79.5 GPa and 3,478 K, but it does not appear in the phase diagram shown in Fig. 1 because it falls inside the stability range of the hcp phase. At this pressure, the hcp phase melts at the slightly higher temperature of 3,487 K.

The solid–solid coexistence lines shown in Fig. 1 are drawn via integration of their slopes, given by the Clausius–Clapeyron relation, starting from the respective triple points. A detailed discussion of these calculations, as well as the latent heats and volumes of melting, is in SI Appendix.

Phase Map of Ultrarapid Solidification

In contrast to the thermodynamic melt line, which delineates the boundary between the liquid and solid phases in the limit of infinitely slow change of thermodynamic variables, we explore in this section the maximal effect of kinetics, by determining the start temperature $$ T_{\text{HS}}$$ for freezing with vanishing nucleation barrier. We have measured $$ T_{\text{HS}}\left(P\right)$$ through nonequilibrium molecular-dynamics (NEMD) simulations of ultrarapid cooling at a rate of 100 K/ns (Fig. 1). We have tracked the solid-phase evolution during these NEMD simulations using the adaptive cutoff common neighbor analysis (a-CNA) (5051–52) method, which allows for identification of the hcp, as well as the fcc and the bcc phases. Fig. 1 depicts by color the solid phases observed at each pressure, with blue depicting the bcc phase and green depicting the close-packed fcc/hcp phases. At 80 GPa and above, pure bcc nucleation is observed, and the bcc nuclei grow and dominate the simulation box. Similarly, at pressures below 40 GPa, pure fcc nucleation is observed, and the close-packed nuclei grow and dominate the simulation box. In the pressure region between 40 and 70 GPa, polymorphic solidification takes place in the simulations. The difference in the observed polymorphic evolution kinetics as pressure is increased from 40 to 70 GPa is quite noteworthy. At 40 GPa, small bcc as well as close-packed (fcc/hcp) nuclei are formed, but the bcc nuclei quickly transform to close-packed phases and vanish. At 60 GPa, pure bcc nuclei do form and grow. However, transformation to close-packed phases during growth is also observed (Fig. 2). In this scenario, at first, small bcc nuclei form. During their growth stage, at or near the bcc solid–liquid interface, transformation to close-packed phases occurs, which starts competing with the original bcc phase as the nucleus grows.

Fig. 2.

A snapshot of an NEMD simulation during solidification at 60 GPa. Particles have been classified as belonging to different crystal phases via the a-CNA method (5051–52): bcc (blue), fcc (green), and hcp (red). Particles identified as liquid have been removed.

It should be noted that $$ T_{\text{H}S}$$ at all pressures is clearly outside the region of thermodynamic stability of the bcc phase. Nevertheless, at pressures 80 GPa and above, not only nucleation occurs in the bcc phase, but the bcc nuclei grow to fill the simulation box. Hence, the solidification process kinetically stabilizes the metastable bcc phase. At pressures below 80 GPa, mixed close-packed and bcc nucleation as well as bcc to close-packed transformation of the growing nuclei can be observed. Hence, the NEMD simulations provide proof for polymorphic solidification kinetics near the triple point of the phase diagram. In the following, we analyze the physical origin of this phenomenon and discuss the conditions for growth of macroscopic-sized metastable bcc crystals, observable in laboratory experiments.

SCL Basins and Thermodynamics of Multiphase Metastable Equilibria

In this section, we discuss the prevalence of solid–liquid metastable equilibria. Multiphase simulations in the NPH ensemble constitute an effective tool for exploring the cores of the MSCL subspaces. Fig. 3 show examples of compact solid-phase nuclei that coexist with the liquid in the pressure range 20 to 100 GPa. They range from single-phase bcc, fcc, or hcp clusters to multiphase mixed fcc/hcp nuclei. These basins are distinct and weakly coupled to each other. It should be noted that the bulk bcc phase becomes thermodynamically stable above 85 GPa. Nevertheless, bcc clusters containing at least several thousand particles can be found stable in MD simulations with durations in the range of 5 to 10 ns, down to 20 GPa. This indicates that phase transformation of these clusters is clearly an activated process with a substantial barrier, which only increases with cluster size since nucleation of a new solid phase involves misfit strain that is easier to accommodate in smaller clusters interfacing the isotropic liquid. Hence, these solid–liquid equilibria constitute the core regions of their respective MSCL subspaces. It should however be noted that most shallow SCL basins may lack metastable subspaces. As is laid out in some detail in SI Appendix, we have not been able to stabilize mixed bcc/close-packed clusters in coexistence with the melt.

Fig. 3.

A–D show four different crystal-phase clusters that can coexist in metastable equilibria with the liquid phase. Color coding is based on phase designation of each particle through the a-CNA method: fcc (green), bcc (blue), and hcp (red).

The study of metastable-phase cluster–liquid equilibria (Fig. 3) reveals many surprising structural features. For instance, in contrast to previous works in Lennard–Jones systems (29, 30, 36), where the fcc clusters were found to interface the liquid through a bcc shell, negligible structural heterogeneities are found at the solid–liquid interfaces for either fcc or bcc clusters in this study. However, when examining the hcp clusters (Fig. 3C), we find that the hcp-(0001) facets interface the liquid through a shell of fcc-(111) layers. This amounts to a shift in the stacking of the close-packed layers of hcp clusters, in favor of the fcc symmetry, in the vicinity of the melt. It is shown in SI Appendix that the energy cost associated with the change in stacking of the solid layers is an order of magnitude smaller than the fcc–liquid interfacial free energy. Later in this paper, we discuss a kinetic mechanism that drives the growth of metastable–fcc phase at (0001)-hcp–liquid interfaces at all temperatures below the fcc melting point. This behavior is quite at odds with the other hcp–liquid interface orientations. It lends credence to the unusual hcp–liquid interface structure, seen in Fig. 3C, as well as the existence of distinct MSCL basins of multiphase solid clusters composed of mixed fcc/hcp stackings (see Fig. 3D).

We conclude this section by pointing out that the existence of a multitude of MSCL basins demands that we quantify their relative thermodynamic stabilities or in other words, their phase space volumes. In the next section, we derive a rigorous and surprisingly general and convenient framework for extracting free energies of solid–liquid metastable equilibria from macroscopic observables extracted from equilibrium two-phase simulations, utilizing the one to one correspondence between the sizes and the coexistence temperatures of equilibrium solid clusters in the melt. It makes no assumption regarding the equilibrium shapes of the crystalline nuclei or the complexity of their crystal structures. Later in this paper, the power of this methodology is demonstrated by applying it to computer simulations of the fcc and the bcc-phase MSCLs of the EAM–Cu model. Consequently, the relative prevalence of these phases during the different stages of the nucleation and the growth processes is quantified. Furthermore, the validity of several standard phenomenological rules that are commonly used to analyze kinetics of solidification in experiments is discussed.

Validation by Computer Experiments.

We now show that while Turnbull coefficients $$ \alpha \left(T_m,P\right)$$ have distinct values for different phases, such as bcc and fcc, they are nearly constant over a wide range of pressures in a typical metal. Furthermore, by explicit calculations of the interfacial energy $$ {\alpha }_{ϵ}$$, we calculate the entropic contribution to the interfacial free energy and demonstrate the accuracy of the negentropic model.

In order to calculate the Turnbull coefficients of the fcc and the bcc phases, we have conducted careful calculations of fcc–liquid and bcc–liquid equilibria in the pressure range from 60 to 100 GPa. For these studies, cubic simulation boxes, containing in excess of 13 million atoms, have been utilized. At each pressure, three two-phase systems at different $$ H_{\text{sys}}$$ are equilibrated for 8 ns each using NPH MD simulations with a time step of 2 fs, and the coexistence temperatures $$ T_c$$ and the systems’ volumes $$ V_{\text{sys}}$$ are recorded. The coexistence temperatures $$ T_c$$ range from 30 to 60 K undercooling below $$ T_m$$. The measured values of $$ T_c$$ and $$ V_{\text{sys}}$$ are inserted into Eqs. 2 and 3 to obtain $$ N_{\text{S}}$$, which in turn, is used in Eq. 6 to obtain the interfacial free energy $$ {\gamma }_I\left(T_c,P\right)$$. The latter is further used in Eq. 9 for evaluation of $$ \alpha \left(T_c,P\right)$$, with $$ A_I^{\text{sphere}}={\left(36\pi \right)}^{1/3}{\left(N_{\text{S}}{v}_{\text{S}}\right)}^{2/3}$$. At each pressure, the calculated $$ \alpha \left(T_c,P\right)$$ are fitted to a linear temperature dependence, which is used to estimate $$ \alpha \left(T_m,P\right)$$.

Fig. 4A shows $$ \alpha \left(T_m\right)$$ in the pressure range from 60 to 100 GPa for both the bcc and the fcc phases. To within the statistical accuracy of our calculations, $$ \alpha \left(T_m,P\right)$$ for both phases can be considered nearly pressure independent, with $$ {\alpha }^{\text{fcc}}\left(T_m,P\right)\approx 0.54$$ and $$ {\alpha }^{\text{fcc}}\left(T_m,P\right)\approx 0.465$$. While this result confirms that the Turnbull coefficient of the fcc phase exceeds that of the bcc phase, their ratio of about 1.16 is significantly smaller than was found in ref. 59; Sun et al. (59) calculated the Turnbull coefficients of fcc and bcc iron using different empirical potential models and found $$ {\alpha }^{\text{fcc}}$$ to be in the range from 0.5 to 0.55 in agreement with this work, but they found $$ {\alpha }^{\text{bcc}}$$ to be in the range from 0.29 to 0.36, which is significantly smaller, leading to a ratio in excess of 1.5. Our result, however, is in closer agreement with the calculated relative Turnbull coefficients of the fcc and the bcc phases according to the polytetrahedral structural model of the liquid (5556–57), as well as particle systems interacting via inverse power law potentials (60, 61). We thus conclude that the value of the Turnbull coefficient for a crystal phase can be quite sensitive to particle interactions. Nevertheless, we find the pressure dependence of the Turnbull coefficient to be weak, in agreement with a recent study of flat fcc–liquid interfaces (62).

Fig. 4.

(A) Calculated Turnbull coefficients for the fcc and the bcc phases as a function of pressure. (B) Ratio of the interfacial energy $$ {\alpha }_{ϵ}$$ to the interfacial free energy $$ \alpha \left(T_m\right)$$ of the fcc and the bcc phases as a function of pressure. (C) Ratio of the latent heats of the fcc and the bcc phases as a function of pressure. (D) Comparison of two different algorithms (gauges) for extracting cluster radii from MD simulations. The filled circles are bcc-cluster radii extracted from two-phase simulations at different temperatures and pressures ranging from 60 to 100 GPa. The open squares are fcc-cluster radii extracted from two-phase simulations at different temperatures and pressures ranging from 60 to 100 GPa.

It should be noted that the preferential nucleation of the bcc phase at large undercoolings is driven by the interfacial free energy, which depends both on the relative Turnbull coefficients as well as the relative latent heats of the bcc and the fcc phases. The latter are shown in Fig. 4C, where it is observed that $$ \mathrm{\Delta }{H}_m^{\text{f}cc}>\mathrm{\Delta }{H}_m^{\text{b}cc}$$. Their ratio is, however, strongly pressure dependent. It is reduced from 1.37 at ambient pressure to 1.12 at 100 GPa. Nevertheless, at all pressures the larger latent heat and the larger Turnbull coefficient of the fcc phase contribute to $$ {\gamma }^{\text{f}cc}>{\gamma }^{bcc}$$.

In order to discuss the temperature dependence of the Turnbull coefficient, we calculate the energetic contribution $$ {\alpha }_{ϵ}$$ as follows:

${\alpha }_{ϵ}\frac{\mathrm{\Delta }{H}_m}{v_{\text{S}}^{2/3}}{A}_I^{\text{sphere}}=H_{\text{sys}}-N_{\text{S}}^c{H}_{\text{S}}^{\text{b}}-\left(N-N_{\text{S}}^c\right)H_{\text{L}}^{\text{b}},$ [10]

where $$ H_{\text{S}}^{\text{b}}$$ and $$ H_{\text{L}}^{\text{b}}$$ are the bulk enthalpies of the solid and the liquid phases. Fig. 4B shows the ratio $$ \frac{{\alpha }_{ϵ}}{\alpha \left(T_m\right)}$$ for both the fcc and the bcc phases in the pressure range from 60 to 100 GPa. Clearly, the energetic part of $$ \alpha $$ is very small, which validates the negentropic model for a typical metal. Also note that the bcc–liquid interfacial energy is negative and smaller than the fcc–liquid energy. The negative interfacial energy poses no conceptual problems, as it is much smaller in magnitude than the interfacial entropy. As a result, the interfacial free energy stays positive. It should further be noted that a prior study of flat solid–liquid interfaces of the Lennard–Jones system found that the entropic content of the interfacial free energy can vary significantly between interface orientations as well as densities (63). Hence, the validity of the negentropic model may be system dependent and hence, should be considered on a case by case basis.

We conclude this section by exploring different gauges for determining the cluster size $$ N_{\text{S}}$$. We have until now only used the zero-volume gauge, which assumes the interface is infinitely thin. This gauge is convenient to use and applicable to any multiphase system with the coexisting phases having distinct specific volumes. In order to reduce the statistical error, elaborate order parameters have been devised in the literature that allow for accurate characterization the solid region within the liquid in two-phase systems. We have used one such method, described previously in ref. 41, as an alternative to the zero-volume gauge for determining $$ N_{\text{S}}$$. In this method, the phase character of each particle is first determined by the a-CNA (5051–52). Next, a coordination filter based on a radial distribution function is applied to remove atoms not belonging to the solid cluster. Finally, an effective radius for the cluster is determined by identifying its surface through Voronoi analysis. Fig. 4D shows the comparison between the radii $$ R_{\text{S}}$$ calculated in this way and those calculated using the zero-volume gauge and the relation $$ R_{\text{S}}={\left(\frac{3}{4\pi }{N}_{\text{S}}{v}_{\text{S}}\right)}^{1/3}$$. The figure shows effective cluster radii extracted from a number of two-phase simulations of bcc- and fcc-phase clusters coexisting with liquid at different pressures and temperatures. Overall agreement is quite impressive. It is, however, interesting to note that the agreement between the two methods for the bcc cluster radii is clearly better than for the fcc ones. Furthermore, the common neighbor analysis seems to slightly underestimate the effective radii of the fcc-phase clusters compared with the zero-volume gauge. It is important to reiterate that the absolute value of the interfacial free energy is subject to the chosen gauge. The only gauge-invariant physical observables are the total system’s free energy, enthalpy, temperature, pressure, and volume. Therefore, it is important to use a consistent gauge when, for example, comparing the free energies of fcc- and bcc-phase nuclei coexisting with the liquid phase at a particular temperature and pressure.

Metastable–Solid Phase Diagrams

We are now ready to develop a general strategy for predicting kinetic stabilization of metastable crystal phases during solidification. For this purpose, we need to 1) quantify the relative phase stabilities of critical nuclei as a function of undercooling and 2) determine for each phase the extent of the MSCL basin or in other words, the range of solid-phase cluster sizes that can be considered metastable. In this way, we can predict the rate of transition from homogeneous undercooled liquid into any particular MSCL basin as a function of undercooling. As a result, we can identify the conditions that lead to metastable-phase growth by recognizing that whenever solidification occurs on a timescale shorter than the rate of transition out of an MSCL basin, the critical nuclei that belong to this basin are likely to grow to large sizes. To demonstrate the efficacy of this procedure, we apply it in the following two sections to the problem of metastable bcc-phase growth from the melt at pressures below 85 GPa, where one of the close-packed phases (fcc or hcp) is thermodynamically stable below the melting temperature (Fig. 1). It should be noted that in the vicinity of the hcp–fcc thermodynamic phase boundary (Fig. 1), stacking-fault free energies become so small that solidification into pure fcc or hcp phase becomes improbable. We conclude this paper with a discussion of intricacies arising from unusual structural features of hcp–liquid interfaces that promote mixed hcp/fcc-phase growth.

Solid–Nucleation-Phase Boundaries.

In this section, we describe the procedure for determining the undercooling temperature $$ T^{*}\left(P\right)$$, below which the nucleation rate of the metastable bcc phase exceeds those of the fcc and the hcp phases of EAM–Cu. Above this temperature, nucleation is dominated by the close-packed phases (fcc, hcp, or a combination thereof) and below it, by the bcc phase. Within the CNT, the rate of nucleation of critical solid clusters at temperature $$ T_c$$ is

$J_{\text{S}}\left(T_c\right)=\sqrt{\frac{\mathrm{\Delta }{G}_{\text{S}}^{{}^{″}}\left(T_c\right)}{2\pi k_B{T}_c}}{\left(N_{\text{S}}^c\left(T_c\right)\right)}^{2/3}\frac{\tau }{v_{\text{L}}}\exp \left(-\frac{\mathrm{\Delta }{G}_{\text{S}}\left(T_c\right)}{k_B{T}_c}\right).$ [11]

Above, $$ \mathrm{\Delta }{G}_{\text{S}}$$ is the excess free energy of a critical solid nucleus coexisting with the melt at temperature $$ T_c$$, and $$ \mathrm{\Delta }{G}_{\text{S}}^{{}^{″}}$$ is its curvature with respect to cluster size fluctuations. $$ \mathrm{\Delta }{G}_{\text{S}}$$ can be defined in terms of the quantities in Eq. 1: $$ \mathrm{\Delta }{G}_{\text{S}}=G_{\text{SL}}-N_{\text{sys}}{G}_L^{\text{b}}$$, which were numerically determined in previous sections. $$ \tau $$ denotes the rate of attachment from the liquid to a unit area of the solid cluster, and $$ v_{\text{L}}$$ is the specific volume of the liquid.

Eq. 11 allows the nucleation-phase boundary $$ T^{*}\left(P\right)$$ to be determined by solving two equations: $$ J_{\text{S}}^{\text{bcc}}\left(T_{\text{fcc}}^{*}\right)=J_{\text{S}}^{\text{fcc}}\left(T_{\text{fcc}}^{*}\right)$$ and $$ J_{\text{S}}^{\text{bcc}}\left(T_{\text{hcp}}^{*}\right)=J_{\text{S}}^{\text{hcp}}\left(T_{\text{hcp}}^{*}\right)$$. Hence, the phase boundary $$ T^{*}\left(P\right)$$ at any given pressure is the smaller of the two solutions $$ T_{\text{fcc}}^{*}$$ and $$ T_{\text{hcp}}^{*}$$ of the equations

$G_{\text{SL}}^{\text{cp}}\left(T_{\text{cp}}^{*}\right)-G_{\text{SL}}^{\text{bcc}}\left(T_{\text{cp}}^{*}\right)=\frac{k_B{T}_{\text{cp}}^{*}}{2}\log \left[\frac{\mathrm{\Delta }{H}_m^{\text{cp}}{\alpha }^{\text{cp}}}{\mathrm{\Delta }{H}_m^{\text{bcc}}{\alpha }^{\text{bcc}}}{\left(\frac{{\tau }^{\text{cp}}}{{\tau }^{\text{bcc}}}\right)}^2\right],$ [12]

where $$ T_{\text{cp}}^{*}$$ represents either of the two temperatures $$ T_{\text{fcc}}^{*}$$ or $$ T_{\text{hcp}}^{*}$$. We have solved the above equations by making two assumptions: 1) $$ {\tau }^{\text{cp}}\approx {\tau }^{\text{bcc}}$$ and 2) $$ {\alpha }^{\text{hcp}}\approx {\alpha }^{\text{fcc}}$$. These are both reasonable assumptions. Assumption 1 is good since the rate of attachment is mainly a property of the liquid. If desired, it can be calculated from size fluctuations of the critical clusters as described in refs. 64 and 65. We have found this unnecessary in the present context because we have found that neglecting the right-hand side of Eq. 12 altogether only leads to a slight rise in $$ T^{*}\left(P\right)$$. This increase is largest at low pressures. At 20 GPa, it amounts to 0.7 K, while at 40 GPa, it is reduced to only 0.25 K. Assumption 2 is reasonable due to the similarity of the two close-packed phases both structurally and energetically. In the pressure range of interest (around 70 GPa), we find that the latent heats of melting of the hcp and the fcc phases differ by only $$ 2\%$$, while that of the bcc phase is $$ 20\%$$ smaller. Furthermore, their respective equilibrium shapes shown in Fig. 3 A and C suggest that $$ {\alpha }_{\text{hcp}}>{\alpha }_{\text{fcc}}$$. This assertion is based on two observations: 1) the equilibrium shapes of fcc clusters in the melt are nearly spherical, from which it can be concluded that the fcc–liquid interfacial free energy only depends weakly on fcc crystal plane orientations, and 2) the equilibrium shapes of hcp clusters in the melt exhibit hcp-(0001) facets that interface the liquid through a shell of fcc-(111) layers. This implies that the liquid–fcc interfaces are favored relative to liquid–hcp ones and thus, have lower free energies.

Eq. 12 and its solutions are discussed in detail in SI Appendix. The resulting phase boundary $$ T^{*}\left(P\right)$$, separating the temperature–pressure regions where bcc nucleation dominates [$$ T<T^{*}\left(P\right)$$] from those where fcc or hcp nucleation prevails [$$ T>T^{*}\left(P\right)$$], is shown in Fig. 5. At pressures $$ P<71.9$$ GPa, $$ T^{*}\left(P\right)$$ (green line) is the nucleation phase boundary between the bcc and the fcc phases, and at $$ P>71.9$$ GPa (red line), it constitutes the boundary between the bcc- and the hcp-phase critical nuclei. Whenever a critical nucleus belongs to an MSCL basin, it is expected to be long lived. It can thus grow to large sizes and be observed during solidification and likely even after the phase transformation has completed. In the next section, we determine the extent of the bcc–MSCL basin upon very strong undercooling by studying the kinetic stability of postcritical bcc nuclei.

Fig. 5.

A solidification–kinetic phase map. The top and bottom solid lines bound the region where solidification can take place. The four lines in order from top to bottom are 1) the thermodynamic melt line $$ T_m$$, with the green part depicting liquid/fcc-phase boundary and the red part depicting liquid/hcp; 2) the boundary $$ T^{*}$$, where the CNT rate of bcc nucleation equals the fcc rate (green part) or the hcp rate (red part) with a critical nucleation bcc–hcp–fcc triple point at $$ P=71.9$$ GPa and $$ T=\mathrm{3,263}$$ K; 3) the boundary $$ T_u$$, below which the postcritical bcc nuclei become unstable and transform to fcc (green part) and mixed fcc/hcp (brown part); and 4) the boundary $$ T_{\text{HS}}$$, below which the liquid phase becomes unstable and crystallization occurs with vanishing nucleation barrier and thus, can be observed in MD simulations. The colors depict the observed crystal phase in the simulations: fcc: green; mixed fcc/hcp: brown; and bcc: blue. The solid parts of the $$ T^{*}$$ and $$ T_u$$ lines bound the region, where the bcc phase is nucleated inside the bcc–MSCL subspace.

Before ending this section, we remark on the fcc–hcp–bcc critical nucleation triple point at $$ P=71.9$$ GPa and $$ T^{*}=\mathrm{3,263}$$ K. The line connecting this point in Fig. 5 to the bulk thermodynamic fcc–hcp–liquid triple point at $$ P=71.6$$ GPa and $$ T_m=\mathrm{3,320}$$ K is the boundary that separates the temperature–pressure regions where nucleation (according to the CNT) is dominated by the fcc phase from those where the hcp phase is most favored. However, we expect kinetic factors beyond those considered in this section to promote mixed fcc/hcp solidification at $$ T>T^{*}\left(P\right)$$ within at least a 10-GPa-wide pressure window about the thermodynamic triple point at $$ P=71.6$$ GPa, due to diminishing energy cost of stacking faults. At first sight, this may seem only marginally relevant to solidification kinetics and have little effect on the liquid–solid interfacial kinetics. However, we will show in a later section that in regions where the CNT predicts hcp-dominated nucleation and growth, the liquid interface with special hcp plane orientations can in fact promote fcc-phase growth by a kinetic mechanism.

Extraneous-Phase Interface Layers: Interfacial Phase Transformations Coupled to Growth-Mode Transitions

Hitherto in this paper, we have discussed kinetic stabilization of metastable single-phase bcc solid at significant undercooling. This is the expected regime based on Ostwald’s step rule, which suggests that small solid clusters in liquid can prefer phases other than the bulk equilibrium phase, due to their interfaces with the melt contributing substantially to their stability. As an example, the line representing $$ T^{*}\left(P\right)$$ in Fig. 5 depicts the minimal undercooling required for nucleation of metastable–bcc phase from the melt for the model Cu system of this study. In this analysis, the interface assumes a passive role. As described earlier in this paper, it suffices to consider the interface as a dividing surface occupying no space and containing no particles, only contributing to the overall free energy of solid nuclei immersed in the melt. This is reasonable as long as the internal structure of the solid–liquid interface can be neglected.

However, the observation of hcp nuclei with heterogeneous interface structures (Fig. 3C) is suggestive of the existence of such phenomenology that cannot be explained via the above approach. In particular, Fig. 3C shows that the (0001)-hcp planes interface the liquid via an fcc shell. This can be rationalized by the observation that the stacking fault free energy is an order of magnitude smaller than the solid–liquid interfacial free energy (SI Appendix).

The observation of a shell structure is not new. In fact, in the earliest computational studies (29) of the Lennard–Jones system, the fcc nuclei were found to be coated by bcc layers at their interfaces with the melt. To our knowledge, the ramifications of this unusual interfacial structure for the kinetics of solidification have not been explored in the literature. As mentioned earlier, we have not been able to observe such bcc shell structures in any of our systems, which suggests that the energy cost of fcc–bcc interfaces must be relatively high in our systems. However, the fcc shell structure at the interface of hcp-(0,001) planes with the melt offers an opportunity to study its unusual features.

Unusual internal structures of solid–liquid interfaces in finite-sized nuclei (e.g., Fig. 3C) inevitably bring forth the question of whether such atomic arrangements are dependent on the interfacial curvature. If so, structural transformations are to be expected at the interfaces of the growing solid nuclei. This is undoubtably a noteworthy phenomenon with possible deep ramifications for crystal growth kinetics. If not the case, the shell structure should appear even in the flat interface limit. In the following, we present an in-depth study of this limit by two-phase simulations of periodic slab geometries for two interface orientations, hcp-(0,001) and hcp-$$ \left(\mathrm{1,1}\overline{2}0\right)$$, in the neighborhood of the melting point $$ T_m=\mathrm{3,497}$$ K at 80 GPa, where hcp is the bulk equilibrium phase (Fig. 1). Based on our previous discussions of Fig. 3C, the hcp-(0001) interface with the liquid should be in the fcc phase, contrary to the hcp-$$ \left(11\overline{2}0\right)$$ interface that stays isostructural with the bulk phase (SI Appendix, Fig. S11). A systematic study of the dynamic properties of these solid–liquid configurations at small undercooling and superheating reveals the remarkable result shown in Fig. 6. It depicts the velocity of each interface as a function of temperature in the vicinity of $$ T_m$$; positive sign of the velocity indicates growing liquid phase, and negative sign implies solid-phase growth. In the narrow range of temperatures studied in Fig. 6, about $$ 1\%$$ deviation from $$ T_m$$, linear response is expected. This is clearly observed for the hcp-$$ \left(11\overline{2}0\right)$$ interface mobility in Fig. 6. Note that due to the relatively small latent heat of the EAM–Cu system (i.e., $$ \mathrm{\Delta }H/k_{\text{B}}{T}_m\approx 1$$), which is typical of simple metals, the interface is atomically rough, and crystal growth is dominated by rates of attachment and detachment of atoms to the crystal surfaces that are controlled by liquid diffusivity, often only thermal velocity, and the latent heat (66). The available sites for attachment on the rough interfaces are large, commonly assumed to be about 1/4 of all of the sites, so-called repeatable step sites (67).

Fig. 6.

Temperature dependence of the mobility of two different hcp–liquid planar interfaces with orientations (0001) and $$ \left(11\overline{2}0\right)$$ at pressure $$ P=80$$ GPa and $$ T_m=\mathrm{3,497}$$ K. The calculations are performed in periodic slab geometries, with dimensions specified in SI Appendix. The (0001) interface kinetics is partitioned into three temperature regimes. 1) $$ T<\mathrm{3,489}$$ K (green solid line with filled circles): The liquid phase transforms to the fcc phase, which grows on the hcp substrate. 2) $$ \mathrm{3,489}\le T\le \mathrm{3,503}$$ K (green dashed line with open circles): No steady-state interface motion is observed. The hcp substrate interfaces the liquid through a film of fluctuating fcc clusters. 3) $$ T\ge \mathrm{3,503}$$ K (red solid line with filled circles): The hcp slab melts at steady rate.

In contrast, the (0001) interface velocity has strongly nonlinear temperature dependence. In particular, in the interval $$ {\mathrm{\Delta }}_a=\left(\mathrm{3,489,3,503}\right)$$ K, interface mobility vanishes on the simulation timescale exceeding 10 ns (dashed green line in Fig. 6). Outside of this interval, interface motion is observed to be normal, diffusive, and linearly changing with temperature variations. Above 3,503 K, the hcp phase shrinks smoothly at a speed proportional to the superheating, as shown in Fig. 6 as well as in SI Appendix, Fig. S13. Below 3,489 K, the interface region is solid fcc and grows steadily at a speed proportional to the undercooling from the fcc melting point at 3,489 K (Fig. 6 and SI Appendix, Fig. S14).

Fig. 7A explains the sudden dynamical arrest of the (0001) interface, observed at 3,503 K in Fig. 6, by invoking a first-order structural-phase transformation at the interface. The order parameter for this transformation is the fcc content of the interface. Above 3,503 K, the concentration of the fcc particles at the interface is dilute, which we designate as phase I. It undergoes a discontinuous jump at 3,503 K, transitioning to phase II, in which larger fcc clusters form and dissolve (SI Appendix, Fig. S12). While the rate of attachment to the fcc sites on the underlying hcp substrate is high, the clusters cannot become critical as long as the temperature is above the fcc melting point of 3,489 K. Further cooling toward the latter temperature causes the fcc content at the interface to steadily rise. As shown in SI Appendix, Fig. S11, the thicknesses of fcc-containing films at the interface are smallest near the transition temperature ($$ \approx \mathrm{3,503}$$ K) and increase upon cooling. Below 3,489 K, metastable fcc phase grows on the (0001) interface.

Fig. 7.

(A) The blue curve depicts the number of fcc particles in the system that mainly reside near the (0001)-hcp–liquid interfaces, and the black curve shows the interface widths $$ W_{\text{int}}$$ as a function of temperature. (B) Time evolution of an NPT MD simulation at $$ T=\mathrm{3,505}$$ K, initiated from an (0001)-hcp–liquid equilibrium slab configuration at $$ T=\mathrm{3,500}$$ K. The blue line shows the evolution of the number of fcc particles, while the black line represents for the same simulation the displacements of the (0001)-hcp–liquid interfaces as a function of time, derived from the size fluctuations of the simulation box. The symmetry character of the particles have been identified via the common-neighbor analysis using a cut-off radius of 2.875 Å. (C) The distribution density function $$ D\left(\overline{\dot{\eta }}\right)$$. The blue curves depict the distribution functions of the phase I atomically rough interfaces with relatively small OPF gradients, and the black curves show the sharp phase II interfaces with much higher OPF gradients. The dashed lines have been obtained by averaging over the first and the last 2,000 time steps of the simulation shown in B.

An important feature of first-order phase transformations is the presence of hysteresis, which implies that metastable phases exist in the vicinity of the transformation temperature. This is confirmed by NPT MD simulations at $$ T=\mathrm{3,503}$$ K, which can equilibrate in either phase I or II depending on initiation. In contrast, Fig. 7B clearly demonstrates that phase II is unstable at 3,505 K. It shows the time evolution of the fcc population at the interface in an NPT MD simulation starting from a configuration in phase II (large fcc content) generated at 3,500 K. During the course of the simulation, the fcc clusters dissolve, and the interfacial structure transforms to phase I (SI Appendix, Fig. S15). Similarly, we find that phase I is unstable at 3,501.5 K. Hence, the hysteresis window is at most 2 GPa.

Fig. 7B also records the motion of the interface with time. The intimate coupling of the interfacial structure to its dynamics is clearly evidenced in this figure. The fcc shell at the interface seems to constitute a protective shield for the hcp substrate, which starts dissolving at around 7,000 ps, when the fcc population has dwindled.

In order to explain the coupling mechanism between interfacial structure and its dynamics, we study the atomic roughness of the hcp–liquid interfaces at different temperatures. It correlates strongly with the width of an interface, which in turn, can be quantified by the extent of the region, where an order parameter field (OPF) exhibits large gradients. In the following, we outline a rigorous method for measuring interfacial thickness from atomistic simulations. As a result, we prove that hcp–liquid interfaces undergo an atomically diffuse to sharp transition concurrent with the structural transformation from phase I to phase II. Consequently, the growth mode changes abruptly from continuous (fast) to layer by layer with high nucleation barrier.

We choose as the order parameter the fraction of the particles in a volume element that are designated via the common-neighbor analysis as hcp. Let $$ \eta \left(z\right)$$ be the OPF along the hcp-(0001) plane normal, and define the OPF gradient $$ \dot{\eta }\left(z\right)$$ as the derivative of this function. It is now straightforward to define and calculate the distribution of OPF gradients from the simulations. For this purpose, it is convenient to work with normalized OPF gradients: $$ \dot{\overline{\eta }}\left(z\right)=\dot{\eta }\left(z\right)/{\dot{\eta }}_{\text{max}}$$, where $$ {\dot{\eta }}_{\text{max}}={\eta }_{\text{hcp}}/\mathrm{\Delta }{z}_{\text{hcp}}$$, with $$ {\eta }_{\text{hcp}}$$ being the value of the order parameter in bulk hcp solid and $$ \mathrm{\Delta }{z}_{\text{hcp}}$$ being the interlayer spacing of the hcp-(0001) planes. The physical representation of $$ {\dot{\eta }}_{\text{max}}$$ is a maximally sharp interface, with a bulk solid layer next to a liquid layer with negligible order. We can now define the distribution density function $$ D\left(\dot{\overline{\eta }}\right)$$ in the interval $$ 0<\dot{\overline{\eta }}<1$$, normalized by the relation

The result of the above analysis is shown in Fig. 7A, where a discontinuous drop in hcp–liquid interface thickness $$ W_{\text{int}}$$ is observed upon phase transformation. Fig. 7C shows the density functions $$ D\left(\dot{\overline{\eta }}\right)$$ for a few interfaces belonging to the two different phases. It is clear that upon transformation from phase I to phase II, the interface becomes smooth and planar and grows layer by layer through two-dimensional nucleation with large activation barrier so that the interfacial motion halts within the computational time scale of this study.

We conclude this section by noting that the structural anisotropy of the hcp phase gives rise to heterogenous growth patterns. In particular, in the system studied here, a growing hcp cluster always contains smooth (0001) facets, on which only metastable fcc phase can grow (Fig. 3C). At the same time, other hcp–liquid interface orientations behave normally and promote the growth of the hcp phase. Hence, interfacial kinetics can drive the growth of multiphase nuclei. This explains the prevalence of mixed fcc–hcp clusters in the pressure range from 60 to 80 GPa (Fig. 3D).

Discussion

In this paper, we have conducted an in-depth study of the kinetic processes during solidification that stabilize crystal phases that are not thermodynamically stable in the bulk. At the nucleation stage, these phases are stabilized by their relatively low interfacial free energies. Additionally, at the growth stage, the kinetic barrier to structural transformation is high and therefore, makes any transition to the thermodynamic ground state improbable. This happens most prominently near the triple points of the phase diagram. By constructing a polymorphic nucleation theory from atomistic first principles and by conducting a rigorous study of the kinetic stability of postcritical nuclei during their growth stage, we have successfully managed to devise a rational framework for characterizing and quantifying the conditions for nucleation and growth of metastable crystal phases during solidification. As a result, in addition to the regular thermodynamic phase diagram, we also construct a kinetic phase diagram that delineates the temperature–pressure regions where metastable phases are favored to grow from the melt.

As a proof of concept, we have applied the above methodology to an atomistic system described by a model EAM–Cu interatomic potential. We thus construct a metastable phase map complementing the thermodynamic phase diagram of this system (Fig. 5). It is quite sobering to observe the significant extent of the metastable region. It is bounded from above by the solid-nucleation phase boundary $$ T^{*}\left(P\right)$$ and from below by $$ T_u\left(P\right)$$, under which the growing bcc nuclei become kinetically unstable. The procedure adopted in this paper assumes near-equilibrium critical nuclei, whose thermodynamic as well as kinetic stabilities are investigated. However, no simplifying assumptions regarding their shapes or structures are made. No doubt, more sophisticated studies of the nonequilibrium kinetics during nucleation are necessary to determine the domain of validity of the assumption of near-equilibrium nuclei (68).

It should be noted that kinetic phase maps are not universal. They depend on the kinetic mechanism that is being considered. Conventionally, kinetic stabilization of metastable phases has been discussed as a result of phase selection during the nucleation stage described within the CNT. In this picture, thermodynamic driving forces govern the phase of a critical nucleus, which according to the CNT, grows to large sizes no fail. In this way, kinetic phase diagrams can be drawn and be interpreted just like thermodynamic phase diagrams. In this paper, we have generalized this idea by including kinetics in a restrictive way as to preserve the notion of a phase diagram. For this purpose, we have introduced the concept of SCL basins and their MSCL subspaces, defined as subdomains of the former that are weakly connected to other SCL basins. Hence, nucleation from the undercooled liquid into a particular MSCL subspace with near certainty is guaranteed to evolve without further phase transformations. The real power of this idea lies in the realization that MSCL subspaces contain all of the large clusters since the barrier to nucleation of solid–solid transformations only increases with cluster size. Hence, if critical nucleation occurs in the transient region of the SCL basin, unless kinetic instabilities occur within a short time after, the nucleus inevitably grows large enough to transition to the MSCL subspace. Hence, thermodynamic conditions that lead to nucleation into the MSCL subspace can be easily identified by exclusion of all those that lead to a finite number of clusters undergoing kinetic instabilities amid their early growths.

While the above is quite a reasonable guiding principle for construction of solidification–kinetic phase diagrams, it is too restrictive and sweeps all real complexities aside by excluding them from the MSCL subspaces. This shortcoming becomes apparent after the growth of metastable fcc in competition with an anisotropic crystal phase such as hcp is considered. In this case, the solid–liquid interface cannot be treated as uniform. Rather, there are special interfaces, whose behaviors are radically different from others. These interfaces can have complex structures distinct from the bulk of the solid clusters and can undergo phase transformations, which in turn, induce sharp transitions in the growth kinetics of competing crystal phases in their neighborhoods. Due to this explicit role of interfacial structure and dynamics, different crystal orientations may promote growth of disparate crystal phases. As a result, the fundamental notion of the kinetic phase boundary as an infinitely thin dividing surface in analogy with the thermodynamic phase boundary must be abandoned and generalized.

Materials and Methods

In this paper, we have presented results of a series of large-scale MD simulations using the LAMMPS code (44) and the EAM interatomic potential for Cu metal, as constructed by Mishin et al. (45). This model represents the bulk physical properties of Cu quite well at ambient condition. While it is not fitted to high-pressure ab initio data, it constitutes a reasonable model for atomic interactions in a generic close-packed metal under pressure. When not otherwise specified, simulations were performed with a 1-fs time step. Visualization and analysis of simulation snapshots were carried out using the OVITO program package (69). Solid-phase designation of each particle has been done via a-CNA (5051–52). Particles are colored according to their phase designation: fcc phase (green), bcc phase (blue), and the hcp phase (red).