Temperature dependence of the solid-liquid interface free energy of Ni and Al from molecular dynamics simulation of nucleation

The solid-liquid interfacial (SLI) free energy, γ, plays a fundamental role in crystal nucleation and growth process.¹ It is also a key parameter required to model the formation of solidification microstructures.² Despite its importance, the measurement of the SLI free energy is extremely difficult in experiments. Therefore, computer simulation, which provides detailed atomistic information, remains heavily employed to quantitatively investigate γ.

A well-established method to compute γ is the capillary fluctuation method (CFM)³ which measures the SLI stiffness based on capillary wave theory.^4,5 While CFM makes an accurate determination of γ, it is only available at the melting point T_m and usually computationally expensive.⁶ To obtain γ at other temperatures, Laird and co-workers further extend the CFM results along the pressure-temperature coexistence curve using the “Gibbs-Cahn integration” method.⁷ However, the temperature dependence of γ at p = 0 remains unclear. Moreover, in the case when several crystal phases compete with each other, a large pressure can trigger a nucleation of the phase which was metastable at p = 0. On the other hand, one can make an indirect measurement of the SLI free energy from nucleation simulation with the classical nucleation theory (CNT).^8,9 This method utilizes the results of molecular dynamics (MD) simulations where the critical nucleus was actually observed. While the method is in principle reliable (see details below), the accuracy strongly depends on the measurement of the size and shape of the critical nucleus.¹⁰ In particular, this method faces the well-known difficulty associated with the fact that the nucleation is usually a very rare event. Recently we developed a persistent-embryo method (PEM)¹¹ to overcome this problem in moderately undercooled liquids. With the PEM, one can observe the actual fluctuations of the large critical nucleus without any biasing. In this work, using the PEM, we determined the average nucleus shape for two fcc crystals, Al and Ni, in the moderately undercooled regime. Then the temperature dependence of the SLI free energy was obtained in the framework of the CNT. These data were used in turn to predict the free energy barrier in a wide temperature range for both systems.

The rest of the paper is organized as follows: in Sec. II, we will introduce the persistent embryo method and provide the simulation details. In Sec. III, we will present the obtained temperature dependences of SLI free energy for Al and Ni. In Sec. IV, we will show that the obtained SLI free energy data lead to the nucleation barriers in agreement with the data determined using a very different technique. In Sec. V, we will discuss the obtained results, and we will provide the summary in Sec. VI.

According to the CNT,¹ a homogeneous nucleation involves a formation of the critical nucleus in the undercooled liquid. The formation of such a nucleus is governed by two factors. The first one is the thermodynamic driving force towards the lower-free-energy bulk crystal. This term is negative and proportional to the number of atoms in the nucleus. The other is the energy penalty for creating an interface between the nucleus and the liquid. This term is positive and proportional to the area of the interface. Therefore, the excess free energy to form a nucleus with N atoms is

where Δμ(<0) is the chemical potential difference between the bulk solid and liquid, γ is the solid-liquid interfacial free energy, and A is the interface area which can be evaluated as $A=sN/ρc2/3$⁠, where ρ_c is the crystal density and s is a shape factor. The competition between the bulk and interface terms leads to a nucleation barrier ΔG^* when the nucleus reaches the critical size N^*, i.e., $∂ΔG(N*)∂N=0$⁠, and

The CNT assumes the spherical shape (⁠$sCNT≡36π3$⁠) for the nucleus to relate ΔG^* with γ and Δμ. This assumption can be lifted by introducing the shape factor s, assuming that the averaged shape of the sub-critical nucleus does not change at the critical size. Mathematically, the interfacial free energy density γ and the shape factor s in Eq. (2), which are both difficult to compute, can be replaced by the critical nucleus size N^* at the critical point¹¹ based on the relation

resulting in $ΔG*=12ΔμN*$⁠. According to Eq. (3), four quantities (ρ_c, Δμ, N^*, and s) are needed to obtain from the MD to calculate the interfacial free energy γ at a given temperature. The determination of the crystal density, ρ_c, is trivial. The chemical potential difference, Δμ, can be calculated by integrating the Gibbs-Helmholtz equation from the undercooling temperature to the melting point.¹² The determination of the critical nucleus size N^* and the shape factor can be obtained from the PEM simulations which will be described in detail below.

The PEM utilizes the main CNT concept that homogeneous nucleation happens via the formation of the critical nucleus in the undercooled liquid. The PEM allows efficient sampling of the nucleation process by preventing a small crystal embryo (with N₀ atoms which are much smaller than the critical nucleus) from melting using external spring forces.¹¹ This removes long periods of ineffective simulation where the system is very far away from forming a critical nucleus. As the embryo grows, the harmonic potential is gradually weakened and is completely removed when the cluster size reaches a sub-critical threshold N_sc(<N^*). During the simulation, the harmonic potential only applies to the original N₀(<N_sc) embryo atoms. The spring constant of the harmonic potential decreases with increasing nucleus size as $k(N)=k0Nsc−NNsc$ if N < N_sc and k(N) = 0, otherwise. This strategy ensures that the system is unbiased at the critical point such that a reliable critical nucleus is obtained. If the nucleus melts below N_sc(<N^*), the harmonic potential is gradually enforced preventing the complete melting of the embryo. When the nucleus reaches the critical size, it has equal chance to melt or to further grow causing fluctuations about N^*. As a result, the N(t) curve tends to display a plateau during the critical fluctuations, giving a unique signal to detect the appearance of the critical nucleus. In addition, multiple plateaus can be collected before a critical nucleus eventually grows, allowing sufficient statistical analysis of the size and shape of nuclei.

All MD simulations in the present study were performed using the GPU (graphics processing unit)-accelerated LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) code.^13–15 The interatomic interaction was modelled using the Finnis-Sinclair potentials¹⁶ developed for the Ni¹⁷ and Al.¹² During the MD simulation, the NPT ensemble was applied with Nose-Hoover thermostats. The damping time in the Nose-Hoover thermostat is set as τ = 0.1 ps which is frequent enough for the heat dissipation during the crystallization (see the supplementary material). The time step of the simulation was 1.0 fs. The simulation cell contained up to 32 000 atoms which is at least 20 times larger than the critical nucleus size. This setting ensures that the effect of pressure change during the nucleation is minimal to the entire simulation box (see the supplementary material).

To identify the nucleus size during the MD simulation, we used the bond-orientational order (BOO) parameter.^18,19 In this approach, one first defines the correlation between the structures of two neighbor atoms i and j as

where

is the Steinhardt parameter, $Ylmr→ij$ are the spherical harmonics, N_b(i) is the number of nearest neighbors of atom i, and $r→ij$ is the vector connecting it with its neighbor j. Two neighboring atoms i and j are considered to be connected when S_ij exceeds a threshold S_c. To choose a reasonable value of S_c, Espinosa et al. suggested an “equal mislabeling” method²⁰ by plotting the population of mislabeled atoms in the bulk solid and liquid as a function of the threshold values. As shown in Fig. 1(a), the crossing point of the mislabeling curves of the bulk liquid and solid phases is chosen as the threshold, S_c, to provide that the probability of mislabeling atoms in the bulk liquid as solid-like atoms is the same as the probability of mislabeling atoms in the bulk solid as liquid-like atoms. This approach works very well when one needs to detect “solid” atoms within a bulk liquid. However, it tends to mislabel “solid” atoms at the cluster interface. To account for that, one can determine the number of solid-like neighbors an atom has. Figure 1(b) shows that this quantity, ξ, is quite different for a majority of atoms in the bulk solid and liquid phases and the number of mislabeled atoms is very small (see the inset in this figure). Intuitionally, it is natural to choose the threshold value, ξ_c, to be 6 for FCC-liquid interfaces. This approach is quite sufficient for the PEM which requires on-the-fly identification of solid-like atoms during the MD simulation. However, recent study shows that the choice of ξ_c considerably affects the value of N^* determined from the MD snapshots.²¹ We will return to this issue in Sec. V.

Figure 2(a) shows a typical PEM simulation. The plateau indicates the appearance of the critical nucleus. Therefore the critical size N^* can be directly measured by averaging the size at the plateau.¹¹ To make a statistically sound description of the nucleus shape, we first averaged the nucleus by superposing the configurations collected in a short time interval (Δt₀ = 10 ps) during the plateau. As shown in Fig. 2(b), the superposed configuration shows a clear non-spherical nucleus shape. Since the crystalline order fades at the interfacial region, it results in a less dense atomic distribution at the outer shell of the nucleus. In order to see the averaged nucleus shape more clearly, a Gaussian smearing scheme^22–24 was applied to convert the atomic distribution into the atomic density in the 3D space. By applying a fast-clustering algorithm²⁵ on the density profile, we were able to extract the high-density points, which are essentially the as-formed crystalline sites. Then the crystalline sites, which were occupied in at-least half of the snapshots collected during the time interval Δt₀, were used to construct the surface of the nucleus by the geometric surface reconstruction method²⁶ integrated in the OVITO (Open Visualization Tool) software package,²⁷ as shown in Fig. 2(b). Finally, the shape factor s was computed based on the surface area A and the volume V of the polyhedron computed from OVITO as $s=A/V2/3$⁠. Figure 2(c) shows the measured shape factor and the critical nucleus size as a function of the temperature for Ni. The shape factor clearly demonstrates a non-spherical shape. However, while the critical nucleus size dramatically increases with the increase of the temperature, the shape factor shows only a slight decrease.

With the measured shape factor and the critical size, the interfacial free energy γ can be calculated by Eq. (3). Figure 3 shows the obtained data for both Ni and Al. In both systems, the interfacial free energy γ shows a nearly linear dependence on the temperature. Therefore, we fit the data with a linear relation to the temperature and extrapolate to the melting point. Figure 3 shows that within the accuracy of the measurement, the extrapolated interfacial free energies agree very well with the data obtained by CFM³ for both Al and Ni.²⁸ 

A straightforward application of the temperature dependence of the interfacial free energy γ is to estimate the free energy barrier at very small and very large supercoolings where the PEM cannot be applied. The case of very small supercooling is interesting because it corresponds to the experimental conditions of solidification. The only way to judge about the reliability of the calculations here is to compare with the experimental data although both experimental and computational data will be affected by the factors not related to the CNT (e.g., the quality of the employed semi-empirical potential in the case of simulation or the presence of impurities in the case of experiment). The case of very large supercooling in the case of pure metals is interesting because the nucleation rate can be directly obtained from the MD simulation. In this case, the quality of the employed semi-empirical potential is not an issue. However, the extrapolation to this temperature range may not work because of several other issues. For example, the temperature dependence of the SLI free energy can be different than the one observed at higher temperatures. Another issue is associated with the fact that the critical nucleus at low temperatures becomes so small that the entire CNT concept may not be applicable.

In the extrapolation of the nucleation barriers [see Eq. (2)], we used a linear fitting for the temperature dependences of the SLI free energy and the shape factor (see Fig. 2). The obtained temperature dependences of the nucleation barriers are shown in Fig. 4. The obtained temperature dependences well describe our PEM data, which was expected because these dependences were obtained by fitting to the PEM data. The question is, if these dependences can be useful to predict the nucleation barrier in a temperature range where the PEM is not applicable? In the case of Ni, the nucleation barrier for the same semi-empirical potential was obtained at T = 1180 K²⁹ using the combination of the mean first-passenger time (MFPT) method^30–32 and the Fokker-Planck equation^30,31,33 directly from an unbiased MD simulation.^34,35 In the present work, we used exactly the same approach to obtain the nucleation barrier for Al at T = 580 K. Figure 4 shows that the obtained MFPT data are in excellent agreement with the data we obtained using the temperature dependences of the SLI free energies.

In the present study, we obtained the temperature dependence of the SLI free energy at the moderate undercooling range where other existing techniques are not applicable. Therefore, to validate the obtained results, we extrapolated the obtained temperature dependences to the temperatures where well-established methods can be applied. Figure 3 shows that the extrapolation to the melting temperature agrees very well with the CFM data. It should be noted that contrary to the CFM which provides the SLI free energy as a function of the interface orientation, in the present study, we obtained the SLI free energy averaged over all orientation using the CNT framework [see Eq. (3)]. Therefore, we compare the current results to the γ₀ value from the CFM [see Eq. (1) in Ref. 3]. This was reasonable for pure Ni and Al since the anisotropy of the SLI free energy is not very large for the pure fcc metals^36,37 at least at the melting temperature. Moreover, the PEM provides ample statistics to measure the shape of the nucleus in the temperature range where it is applicable, and in the present work, we did not observe a very large deviation from the spherical nucleus shape. However, one should be cautious in the interpretation of the SLI free energy value obtained from the PEM in the case crystal phase with very anisotropic SLI free energy (e.g., see Fig. 10 in Ref. 38).

Another possibility to validate our results was to extrapolate the obtained temperature dependences to low temperatures and compare the obtained nucleation barrier free energies with the data obtained from the brute force MD simulations. The obtained excellent agreement is rather surprising because it suggests that the CNT still works at these temperatures despite the fact that critical nucleus size (only tens of atoms²⁹) is so small that it is not really possible to distinguish between the bulk and the interface regions within a nucleus. In this case, even the concept of the SLI free energy is not clear. Yet, one can always describe the change in the free energy associated with the nucleus formation as the sum of two contributions: the product of the difference in the bulk free energy per atom and the number of atoms in the nucleus and a contribution, which accounts for the nucleus interface. The latter can be treated as the flat interface free energy corrected for the high interface curvature (e.g., see Refs. 35, 39, and 40). In fact, this is the quantity we obtained from the PEM. At high temperatures, where the nucleus is large and the correction for the high interface curvature is negligible, we obtained good agreement with the flat interface free energy data from the CFM. At low temperatures, we obtained good agreement with the brute-force MD simulation data, but the quantity we extracted includes not just the flat SLI free energy but also corrections associated with the SLI curvature. The authors of Ref. 41 argued that these corrections explain why the value of the SLI free energy obtained from the seeding simulations is always below that estimated from the Turnbull correlation⁴² which was proposed in Ref. 34 to use to estimate the temperature dependence of the SLI free energy. The temperature dependences of the SLI free energy obtained in the present study are also below the predictions based on the Turnbull correlation (see Fig. 3).

The main source of the uncertainty in the determined value of the SLI free energy comes from the uncertainty in determination of the number of atoms in the critical crystal cluster, N^*. This quantity can be rather sensitive to the choice of order parameters, as has been noted in Refs. 21 and 43 and can be seen in Fig. 5. In addition to the BOO parameter, we also employed the cluster-alignment (CA) method²³ in which minimal root-mean-square deviations (RMSD) between the atom cluster and the perfect packing templates such as FCC, HCP, and BCC polyhedral are calculated for crystal-structure recognition. Interestingly, the CA order parameter leads to almost identical results comparing to the use of the BOO parameter with ξ_c = 6 which was assumed to be the most reasonable value.

Figure 6 shows how the uncertainty in N^* caused by the choice of the order parameters propagates in the uncertainty of the SLI free energy determined within the present study. A vivid systematic difference can be seen. However, it is important that the temperature dependence remains qualitatively the same: no matter what order parameter we used, the obtained temperature dependence was linear. What is even more important is that all lines come to almost the same point which is in excellent agreement with the CFM value of the SLI free energy.

The extrapolation to low temperatures is shown in Fig. 7. The obtained results indeed depend on the choice of the order parameter, and extreme choices can lead to considerable overestimations or underestimations of the nucleation barrier. However, the reasonable choice of the threshold value in the BOO parameter (ξ_c = 6) or using the CA order parameter provide excellent agreement with the brute-force MD simulation. Moreover, using slightly different values of the threshold value in the BOO parameter (ξ_c = 5 or ξ_c = 7) lead to the variations in the nucleation barrier value obtained by extrapolation of the PEM data within uncertainty of the brute-force MD simulation data.

In summary, based on the atomic configurations of critical nucleus obtained by the PEM method, we estimate the shape factor of the nucleus and find that both critical nucleus of Ni and Al deviate from the spherical shape. With the framework of CNT, we obtain a nearly linear temperature dependence of the orientation-averaged interfacial free energy for both Ni and Al. Using this temperature dependence, we predict the free energy barrier in a wide temperature range, which shows good agreement with the value obtained from brute-force MD simulations.

See supplementary material for the latent heat dissipation, the analysis of the statistical uncertainty, and the pressure on the nucleus.

The work at Ames Laboratory was supported by the US Department of Energy, Basic Energy Sciences, Materials Science and Engineering Division, under Contract No. DE-AC02-07CH11358, including a grant of computer time at the National Energy Research Supercomputing Center (NERSC) in Berkeley, CA. K.M.H. acknowledges support from USTC Qian-Ren B (1000-Talents Program B) fund. The Laboratory Directed Research and Development (LDRD) program of Ames Laboratory supported the use of GPU computing.