Formation energetics/dynamics of icosahedral clusters in supercooled metallic liquids in the dynamic equilibrium regime: Gibbs free energy, entropy, enthalpy, and connection to coordination shells

Icosahedral (ICO) clusters are known to exist in many supercooled metallic liquids and believed to play an important role in stabilizing the liquid before it transitions into a glassy, crystalline or quasicrystalline solid. However, a detailed understanding of their formation energetics/dynamics is currently lacking and a set of key questions regarding these clusters remains to be answered. Here, we report our study on the formation energetics/dynamics of ICO clusters in liquid Cu64Zr36 and Ta by combining MD simulations with statistical and theoretical analysis. We present the formation Gibbs free energy, entropy, enthalpy of ICO clusters in the two liquids in the dynamic equilibrium regime (T > 0.75 Tm), determine the size of the spatial domain (number of coordination shells) surrounding the clusters from which the formation enthalpy is originated, and discuss the results in connection with liquid composition, degree of randomness, potential energy landscape, and glass transition.


Introduction
Icosahedral (ICO) clusters, each having one core and 12 shell atoms arranged with a fivefold symmetry (ideally), were first predicted by Frank in 1952 to be a local atomic structure potentially favored by many liquids [1]. To date, widespread presence of ICO clusters among supercooled metallic liquids has indeed been substantiated by both computation (mainly molecular dynamics, MD, simulations [2][3][4][5][6][7][8]) and experiments (including neutron scattering and synchrotron X-ray scattering [9][10][11][12][13][14]). The ICO clusters in supercooled liquids are believed to be closely related to the formation of metallic glasses and quasicrystals which are usually manufactured by liquid solidification. These non-conventional, amorphous, or quasicrystalline, metals/alloys have been found to contain disorderly packed (metallic glasses) [2][3][4][5][6][7][8][15][16][17] or elegantly assembled (quasicrystals, with rotational symmetry but no or reduced translational symmetry) [18][19][20][21] ICO clusters which are most likely inherited from the liquid precursor. ICO clusters are also known to exist in some crystalline alloys with topologically close pack (TCP) structures (e.g., [21,22]). Even for more conventional crystalline metals/ alloys with relatively simple structures (e.g., body-centeredcubic, face-centered-cubic) manufactured by liquid solidification, the ICO clusters in the liquid may also play an important role (delaying) in the incubation, nucleation, and early growth stages of the crystal phase(s) [4,10,11], and hence affect the resulting microstructure and properties. Therefore, understanding the ICO clusters in the liquid phase has broad relevance to physical metallurgy and materials science.
It is well established that the presence of ICO clusters stabilizes the liquid state and that the ICO cluster population in a liquid increases with decreasing temperature. However, a more detailed understanding of the formation energetics/dynamics of these clusters is currently lacking. A number of intriguing and practically important questions regarding these clusters are still awaiting answers, for example: Is the ICO cluster formation a thermally activated process? If so, why does the ICO population increase with decreasing temperature? If not, why is their population limited (and changing with temperature)? How is the formation energetics/dynamics of ICO clusters related to the potential energy landscape near the clusters? How is the continuous increase in the ICO population logically connected to glass transition?
In this paper, we report our study on the formation energetics/dynamics of ICO clusters in liquid Cu 64 Zr 36 and Ta by combining MD simulations with statistical and theoretical analysis (technical details of the approach are provided in Methods after the main text). These two materials are particularly interesting as they have excellent glass-forming ability compared with other binary and unary metals/alloys, respectively [23,24]. We limit ourselves to the dynamic equilibrium (> 0.75 T m here, including thermodynamic equilibrium) regime where the stable state of the liquid can be reached at cooling rates amenable to MD simulations. As follows, we will present the formation Gibbs free energy, entropy, enthalpy of ICO clusters in the two liquids, determine the size of the spatial domain (number of coordination shells) surrounding the clusters from which the formation enthalpy is originated, and discuss the results in connection with liquid composition, degree of randomness, potential energy landscape, and glass transition.

Results and discussion
Overall cooling behavior and the dynamic equilibrium regime  Invited Feature Paper and 1 K/ps for Cu 64 Zr 36 , 1 and 10 K/ps for Ta, are all sufficient to prevent crystallization during the cooling process. In each case, the liquid goes through the thermodynamic equilibrium liquid (T > T m ) regime, the supercooled liquid (T < T m ) regime, and then, at a deep undercooling (0.5-0.7 T m ) transitions into a solid glass. The visible divergence between the red (dashed) and the blue (solid) lines in Fig. 1(a) and (b) near the end of the cooling process reveals the effect of the cooling rate on the glass transition temperature (or, temperature range) and the resulting state variables of the solid glass. However, above 0.75 T m (shaded areas), the converging lines in the figures indicate the establishment of a stable liquid state independent of the used cooling rates. Further reduction in the cooling rate (unless so low as to incur crystallization) will no longer affect this higher temperature regime. Therefore, the liquids have reached dynamic equilibrium in this higher temperature regime which covers the relatively shallow supercooled liquid regime (0.75 T m < T < T m here) and the thermodynamic equilibrium liquid regime (T > T m ). It is noted that the lower border of the dynamic equilibrium regime may be extended somewhat as the cooling rate is further reduced (at the cost of more computational time). However, using the data above 0.75 T m here is appropriate and adequate to probe the main features of the formation energetics/dynamics of ICO clusters in the dynamic equilibrium regime. Figure 2(a) and (b) exhibits the dynamic evolution of the ICO cluster population f ico (defined by the ratio of the counts of ICO core atoms to the total number of atoms in the system) in the cooling processes. In each case, f ico increases continuously as temperature drops. The increase appears to accelerate around the glass transition (0.5-0.7 T m ) and then slowed down afterward. For each material, a lower cooling rate results in a noticeably higher ICO population in the solid glass, but the f ico trajectories in the dynamic equilibrium regime (shaded areas) are converged with respect to the cooling rates used, attesting again to the establishment of a stable liquid state in this regime.
As can be seen from Fig. 2(a) and (b), and Table 1, the ICO population in the monatomic Ta liquid exceeds that in the binary Cu 64 Zr 36 liquid at a common homologous temperature (T homo = T/T m ) in the early dynamic equilibrium regime (T homo > 0.82). However, the trend is reverted in the late dynamic equilibrium regime (T homo < 0.82) due to a faster increase of f ico in the binary liquid there. This alludes to different formation energetics/dynamics of ICO clusters in the two liquids.

Formation Gibbs free energy
In the dynamic equilibrium regime, the population (number fraction) of ICO clusters f ico is directly determined by the Gibbs free energy of formation (G f ) as f ico = exp − G f k B T , similar to the equilibrium concentration of point defects in crystalline metals. Hence, Gibbs free energy (per cluster) of ICO formation is obtained by G f = −k B Tln f ico , as shown in Fig. 3. For both liquids, the G f is positive and decreases with decreasing temperature. The positive sign of G f suggests that the formation of ICO clusters from the liquid matrix is a thermally activated process.
The decreasing G f during cooling shows that the activation barrier for ICO cluster formation is diminishing at lower temperatures. This brings out the following three important insights.
First, it explains why the ICO cluster population increases even though the temperature is dropping. Normally, a thermally activated process would be facilitated by higher temperatures and suppressed by lower temperatures. The diminishing activation barrier found here flips this normal trend and makes the ICO cluster formation more favored at lower temperatures.
Second, it implies that the apparent acceleration of the increase in the ICO population near the glass transition (see Fig. 2) is actually not caused by the glass transition, but rather an inherent trend following the dynamic equilibrium behavior. The glass transition is in fact a departure from the dynamic equilibrium behavior, interrupting the growth of the ICO population.
Third, it provides a new alternative explanation for ideal glass transition, an event that would theoretically occur when a liquid was cooled infinitely slowly and yet no crystallization or quasi-crystallization was going to occur. Traditionally, the ideal glass transition is explained by the "Kauzmann paradox" [25] as a necessity in order to prevent an allegedly strange situation where the entropy of a supercooled liquid falls below that of the crystal counterpart. The Kauzmann criterion may become ambiguous, however, when there are multiple possible crystal counterparts (i.e., different crystal phases) or even quasicrystal counterparts for the same material that possess different entropies. Here, just focusing on the liquid itself, if the dynamic equilibrium regime were allowed to continue indefinitely upon cooling, the activation barrier for ICO cluster formation would drop to zero and then become negative, which would lead to f ico > 1, a structural catastrophe that is clearly not allowed. Thus, an ideal glass transition must occur to end the dynamic equilibrium behavior and prevent the structural catastrophe in this theoretical infinitely slow-cooling scenario. Real glass transitions (such as those shown in Figs. 1 and 2) occur sooner (at higher temperatures) upon cooling due to earlier departure/ deviation from the dynamic equilibrium behavior. Figure 3 also shows that at a common homologous temperature, the monatomic Ta liquid has a higher G f , i.e., higher thermal activation barrier for ICO cluster formation, than the binary Cu 64 Zr 36 liquid. However, it must be noted that the T m is vastly different for the two materials, 1233 K for Cu 64 Zr 36 , 3290 K for Ta. The comparison of f ico in the two liquids at common homologous temperatures as discussed in Sect. 2.1 is co-determined by both the G f and the T m as

Formation entropy
The formation entropy (per cluster) S f of ICO clusters can be calculated using S f = − dG f dT since the MD simulations here are performed under constant pressure (1 bar). As plotted in Fig. 4, the S f is negative for both liquids, reflecting the fact that the local ordering is enhanced when an ICO cluster is formed from the liquid matrix and hence the system loses certain degree of randomness. The S f is almost constant in liquid Cu 64 Zr 36 , while it increases (absolute value decreases) with decreasing temperature in liquid Ta. The latter trend is indicative of a diminishing effect of ICO cluster formation on the degree of randomness in liquid Ta upon cooling, which is likely because the monatomic liquid matrix (the bulk of the liquid excluding ICO clusters and related atoms) relaxes more efficiently (in terms of entropy here) than the binary liquid matrix. It is helpful to briefly relate the formation entropy of ICO clusters here to Richard's rule for entropy of fusion. According to Richard's rule [26], common metals have an entropy of fusion approximately 2.2 cal/(K mol), i.e., 9.5 × 10 -5 eV/(K atom). The latter number is essentially the (absolute value of) formation entropy (per atom) of a crystal (common structures in Richard's rule) from the liquid. A simple comparison with the ICO formation entropy is made in Fig. 4 by scaling this number (negative sign added) with the ideal number of atoms in one ICO cluster, 13. The ICO formation entropy in both liquids lies clearly above this crystal formation entropy, with lower absolute values. This is consistent with the expectation that creating an ICO cluster causes a smaller reduction in the degree of randomness than creating a crystal (with the same number of atoms) as the latter has higher symmetry and order. It is noted that this simple comparison is meant to put the ICO formation entropy obtained here into a general context, not to exactly compare it with the crystal structure formation entropy in the two specific materials.

Formation enthalpy
With the formation Gibbs free energy G f and formation entropy S f calculated, the formation enthalpy H f of ICO clusters is obtained by H f = G f + TS f , and the result is shown in Fig. 5. The H f is negative for both liquids, meaning that the formation of ICO clusters causes a decrease in the enthalpy of the system. The H f is nearly independent of temperature in liquid Cu 64 Zr 36 , while it increases (absolute value decreases) with decreasing temperature in liquid Ta. These main features of H f are similar to those of S f discussed above. The stronger temperature dependence of H f in liquid Ta is again indicative of the more efficient relaxation (here in terms of enthalpy) of the monatomic liquid matrix than the binary liquid matrix upon cooling.
Fundamentally, the enthalpy H of an atomic system is related to the internal energy U, pressure P, and volume V, as H = U + PV . Under a normal pressure (here 1 bar), the PV term for metallic liquids or solids has a much smaller (by 3-4 orders of magnitude) absolute value than the internal energy, and hence H ≈ U = E p + E k , where E p is the potential energy and E k the kinetic energy of the atoms. Under NPT (constant/ controlled number of particles, pressure and temperature) conditions as used in this study (in accordance with typical experiments), formation of ICO clusters has barely any impact on the atomic kinetic energy (the kinetic energy overall is regulated by the thermostat in MD, analogous to the thermal environment in experiments; our statistics reveals that the kinetic energy difference between ICO-related atoms and the matrix is one to two orders of magnitude lower than the difference in the potential energy and it has a random positive or negative sign at different temperatures, similar to thermal fluctuation). Therefore, the formation enthalpy H f is essentially the change in the local potential energy at/around the ICO clusters, with respect to the rest (matrix) of the liquid. However, it is not clear how large this "local" spatial domain really is that determines the formation enthalpy of the ICO clusters, whether it encompasses the ICO core and shell atoms (1st coordination shell) only, or also includes some non-nearest neighbors in the 2nd, 3rd … coordination shells.
In Fig. 6(a), we compare for liquid Cu 64 Zr 36 the ICO formation enthalpy H f with the statistically averaged local potential energy change per cluster using an effective cluster domain (territory) size of 3.8 and 8.5 Å out from the ICO cores. These two cutoff distances represent the borders of the 1st and the 3rd coordination shells, respectively, as determined from the (valleys of) radial distribution function (RDF). It is evident from Fig. 6(a) that, when the cluster domain is considered only up to the 1st coordination shell, the potential energy change per cluster is significantly more negative than the H f . When the cluster domain is considered up to the 3rd coordination shell, the potential energy change per cluster matches the H f quite well. This means that the energetics/dynamics of ICO cluster formation in the liquid is not solely determined by the ICO core and shell atoms, but also by more neighboring atoms, up to the 3rd coordination shell for liquid Cu 64 Zr 36 . Similar conclusions can be drawn for liquid Ta, as shown in Fig. 6(b), although a smaller cluster domain up to the 2nd coordination shell (cutoff distance: 6.5 Å) turns out to be the best match for the formation enthalpy in the monatomic liquid. The smaller effective ICO cluster domain is suggestive of the greater ability of the liquid matrix to dampen the impact of a local ICO cluster, which aligns with the expectation for the monatomic liquid based on its chemical simplicity. In a binary or multicomponent liquid, the formation of an ICO cluster involves not only topological but also chemical rearrangements of atoms, and hence its impact is expected to be more extensive in space.

Distinctions from previous studies
Several significant differences between this work and previous studies in the field are worth pointing out. First, the previous treatment of a liquid structure considers every atom as a core of some cluster (and determines a cluster type for each atom), but neglects the fact that each cluster has shell atoms simultaneously belonging to other clusters. The different types of clusters defined this way are inherently dependent on each other and should not be treated/analyzed all at once as if they were independent. Here, we focus on the most important type of cluster (ICO in the liquid Cu 64 Zr 36 and liquid Ta) and treat everything else as a liquid matrix, hence eliminating the ambiguity in the definition/treatment of clusters. Second, some previous studies [31] defined "formation energy" of a cluster to be the difference between the cluster potential energy (per atom) and the chemical potential in reference crystal structures. In this work, all the different forms of formation energies (e.g., free energy, enthalpy) are defined relative to the liquid matrix, i.e., using the liquid matrix as the reference state. This choice is more suitable for the purpose of understanding ICO cluster formation mechanisms/dynamics in metallic liquids, as the clusters are forming directly from the liquid matrix, not from crystal phases. Finally, to the best of our knowledge, this is the first time that Gibbs free energy, enthalpy, and entropy of ICO formation are determined and interpreted in relation to atomic potential energy and coordination shells.

Conclusions
We have studied the formation energetics/dynamics of ICO clusters in liquid Cu 64 Zr 36 and liquid Ta in the dynamic equilibrium regime by combining MD simulations with statistical and theoretical analysis. In both liquids, we find that the ICO formation Gibbs free energy is positive and diminishes with decreasing temperature. This shows that ICO cluster formation is a thermally activated process but the activation barrier diminishes upon cooling, which explains why the ICO cluster population increases with decreasing temperature. A structural catastrophe with ICO fraction exceeding 1 would occur if the dynamic equilibrium behavior were allowed to continue indefinitely at deep undercoolings. This provides a new alternative explanation for ideal glass transition at the theoretical limit of infinitely slow cooling. Both the formation entropy and the formation enthalpy of ICO clusters in the liquids are negative, demonstrating that the ICO formation reduces both the degree of randomness and the configurational energy of the liquids. The ICO formation entropy and formation enthalpy display a stronger temperature dependence (with absolute values decreasing with declining temperature) in liquid Ta than in Cu 64 Zr 36 , suggestive of more efficient relaxation of the monatomic liquid matrix (reducing the impact of ICO formation on system entropy and enthalpy). In both liquids, the ICO formation enthalpy is less negative than the potential energy difference between the clusters themselves (core plus the first coordination shell) and the liquid matrix, implying that the ICO formation enthalpy is originated from a larger effective cluster domain/territory that goes beyond the first coordination shell. By matching the local potential energy change with the formation enthalpy, the size of the effective cluster domain is determined to be the 3rd and the 2nd coordination shells for liquid Cu 64 Zr 36 and liquid Ta, respectively. The larger effective cluster domain in the binary liquid is indicative of a more spatially extensive impact of ICO cluster formation as can be expected based on the involvement of both topological and chemical rearrangements therein. The study sheds new light on the ICO cluster formation/dynamics, liquid structural dynamics, glass transition, and demonstrates a way to analyze and interpret the energetics of atomic-clusters in a dynamic equilibrium system.

Methods
Molecular dynamics (MD) simulations of thermal processes, including heating, isothermal relaxation, and cooling, are performed for the binary Cu 64 Zr 36 alloy and pure Ta (with 108,000 and 432,000 atoms, respectively) using the LAMMPS code [27] and EAM potentials [28,29]. The NPT (controlled particle number, pressure, 1 bar, and temperature) ensemble and periodic boundary conditions are applied in all the processes. The thermodynamic equilibrium liquid phase is first obtained and relaxed at 2000 K for Cu 64 Zr 36 and 4000 K for Ta, both above T m (Cu 64 Zr 36 : 1233 K; Ta: 3290 K). The preparation of the Cu 64 Zr 36 liquid, specifically, involves heating and melting a pure Cu crystal, randomly replacing 36% of the atoms in the molten Cu with Zr and subsequently relaxing the binary liquid at 2000 K. The thermodynamic equilibrium liquid for each material is then cooled down toward 300 K, during which the liquid goes through the supercooled liquid state and later solidifies into a metallic glass. The potential energy, kinetic energy, and spatial coordinates of each and all the atoms are exported during the cooling process. The cooling simulations are performed at two different cooling rates for each material. As shown in the main text, the system state as indicated by the statistical quantities (e.g., internal energy, population of ICO clusters) in the dynamic equilibrium regime (> 0.75 T m ) is independent of the cooling rates used.
Subsequent to the MD simulations, the Ovito (Open Visualization Tool) program [30] is used for visualization, ICO cluster identification, structural analysis, and atomic energy statistics. The ICO cluster cores are identified using the Voronoi tessellation algorithm built in Ovito and specifically a Voronoi index of < 0 0 12 0>. The theoretical formulations used for calculating the ICO formation free energy, entropy, and enthalpy are provided in the main text. To find the change in the local potential energy caused by the formation of ICO clusters (so as to compare with the formation enthalpy), ICO cluster-related atoms are selected by applying a target cutoff radius to the neighbors of the ICO core atoms, and the liquid-matrix atoms are identified by removing the ICO cluster-related atoms from the system. The atomic potential energy (per atom) is averaged within each group, and the difference (per atom) between the two groups is then scaled by the averaged number of related atoms per cluster (i.e., the number of ICO cluster-related atoms divided by the number of ICO core atoms). The cutoff radius used in the selection of cluster-related atoms is determined from the (valleys of) RDF. For example, the shell atoms of icosahedral clusters can be identified by using the first nearest neighbor distance (the distance at the valley between the first and the second RDF peaks).