Atomistic study of liquid fragility and spatial heterogeneity of glassy solids in model binary alloys

Abstract

Fragility is a fundamental property of glass-forming liquids. Here, we evaluated the liquid fragility and structural and dynamic heterogeneity of glassy solids for four model binary alloys. The most fragile alloy exhibited the maximum dynamic heterogeneity in the mechanical unfreezing process. The local atomic order contributed to structural and dynamic heterogeneities in the glassy solid. We observed that atomic displacement significantly correlated with degrees of clustering of local atomic orders. The clustering produced during the glass-forming quenching process enhanced structural and dynamic heterogeneities, especially in fragile glass alloys. Therefore, this alloy system exhibited correlations among liquid fragility, dynamic heterogeneity in liquid alloys, and dynamic and structural heterogeneities in glassy solids. We discussed the underlying physics of the correlation based on a theoretical model for fragility. These structural and dynamic analyses also provided deeper insights into the features of structural heterogeneity in glassy solids. The alloy with the most fragility exhibited the largest difference in atomic mobility between the densely and loosely packed local atomic orders, implying the greatest heterogeneity in the degree of packing density.

Introduction

Liquid alloys undergo a glass transition upon cooling from the high-temperature liquid when the cooling rate is sufficiently fast to avoid crystallization. This indicates the formation of solid metals with amorphous atomic configurations known as metallic glasses (or amorphous metals). The viscosity of the liquid alloy drastically increases during cooling with decreasing temperature. This sdependency of viscosity on temperature has been discussed using the concept of fragility^1,2. Fragility is correlated with glass-forming ability (GFA)³, elastic properties⁴, and plasticity⁵ of an alloy. In terms of log₁₀(viscosity) and ${T}_{{\rm{g}}}$/T (T and ${T}_{{\rm{g}}}$ are the temperature and the glass transition temperature, respectively), a liquid whose log₁₀(viscosity) value increases linearly with increasing ${T}_{{\rm{g}}}$/T is called a strong liquid. This linear increment indicates a constant activation energy (Arrhenius temperature dependency). In contrast, a liquid showing a nonlinear increment of log₁₀(viscosity) with increasing ${T}_{{\rm{g}}}$/T is called a fragile liquid, and this trend implies an increase in activation energy at lower temperatures (and thus non-Arrhenius temperature dependency). Fragility is correlated with certain structural features^6,7,8 and dynamic heterogeneity^9,10,11, both present in supercooled liquids.

Various studies have been conducted on the glass transition and fragility for decades. Event-driven Brownian dynamics simulations of self-propelled hard spheres showed that with increased activity, the relaxation dynamics could be significantly accelerated¹². Theoretical studies revealed the effects of repulsion and anharmonicity of the interatomic/interparticle interaction on viscosity and fragility^13,14,15. A theoretical model based on the steepness of repulsion suggests that interatomic interactions are one of the factors affecting the fragility of glass systems^13,14. The theoretical model also proposed a correlation between outer shell electronic structures and fragility. The application of these theoretical models to vast real glass systems is under discussion¹⁶. A study of Cu-Zr glasses study also investigated the relationship between thermal expansivity, fragility, and GFA¹⁷. Moreover, the correlations between the average linear thermal expansion coefficient and T_g (and between shear modulus and T_g) indicated that the glass transition of metallic glasses follows a Lindemann-like rule¹⁸. Understanding the glass transition phenomenon of metallic glasses based on Lindemann-type melting was also suggested based on experimental heat capacity measurements¹⁹. An experimental study also indicated that the fragility in the supercooled liquid state could be determined from the frozen behavior of interstitial-type defects during the glass transition²⁰. Recently, the static (or structural) heterogeneity of glassy solids has attracted much attention in experimental^{2,21,22,23,24,25} and computational studies^26,27,28 because heterogeneity affects the material properties of metallic glasses. Although structural heterogeneity in glassy solids could be correlated with liquid fragility, the relationship between heterogeneity and fragility remains unclear. Because local orders developed in the glass-forming quenching process should be a key to the fragility of liquids and structural heterogeneities in glassy solids, a unified understanding of them would provide deeper insights into the GFA and material properties of metallic glasses.

This study aimed to reveal the correlations among liquid fragility, dynamic heterogeneity in liquid and glassy solids, and structural (i.e., static) heterogeneity in glassy solids. Four binary Cu-Zr alloy models with different Cu contents were investigated using molecular dynamics (MD) simulations based on simple pairwise interatomic potentials. These four systems were chosen because they exhibited different fragilities and internal structures. In this study, we focused on the relationship among liquid fragility, dynamic heterogeneity in liquid and glassy solids, and structural heterogeneities in glassy solids for the model alloy system instead of precisely simulating the properties of real Cu-Zr alloys. We evaluated these properties during two processes, freezing a liquid (i.e., glass transition) by quenching and unfreezing a glassy solid by mechanical processing. The cooling process from the liquid state provided information on fragility and dynamic heterogeneity in the liquid. Meanwhile, mechanical unfreezing from a glassy solid provided information on the dynamic heterogeneity correlated with heterogeneous local atomic orders. Dynamic mechanical analysis (DMA)^29,30,31,32 at various temperatures yielded dynamic features in the glassy solids and during their phase transition to supercooled liquids (i.e., the unfreezing process or the inverse of the glass transition). In this study, we focused on the evolution of structural features beyond short-range order (SRO), herein called clustering, and revealed a significant correlation between clustering and the dynamic heterogeneity in liquids and glassy solids. Based on atomistic information from the two simulated processes, we explained the correlations between liquid fragility and structural and dynamic heterogeneities in glassy solids. Furthermore, we discussed features of structural heterogeneity in glassy solids from two viewpoints, the volume fractions of densely and loosely packed regions and the degree of densely and loosely packed states in these regions.

Materials and methods

We used binary metallic glass models composed of Cu and Zr atoms. Employing an MD framework³³, we prepared four Cu_xZr_100-x alloys (x = 30, 40, 52, and 58). Instead of precisely simulating the properties of real Cu-Zr alloys, we focused on correlations among liquid fragility, dynamic heterogeneity in liquid and glassy solids, and structural heterogeneities in glassy solids for the model alloy system. The interatomic interactions were calculated using pairwise potentials³⁴ that had been used for Cu-Zr metallic glasses. We first calculated the α-relaxation time τ_α^8,9,11 over a wide temperature range from the high-temperature liquid to slightly above the glass transition temperature based on the self-intermediate scattering function and estimated the fragility. For each alloy composition, we conducted three independent simulations with different initial velocities and atomic configurations to evaluate the effects of these initial settings on the fragility and freezing behaviors. We also constructed metallic glass models by melt-quenching (MQ) from 2000 to 0 K at a constant cooling rate of 10¹¹ K/s. Specifically, the alloys were first equilibrated in the liquid state at 2000 K. Then, they were frozen into a glassy solid during cooling, indicating that a glass transition had occurred. The constructed glass models were relaxed at arbitrary temperatures for 0.1 ns. We then performed DMA simulations at each temperature by applying a time-dependent sinusoidal strain γ_Asin(2 πt/t_p) at γ_A = 0.02 and t_p = 1 ns (i.e., frequency = 1 GHz). DMA can provide atomic-scale dynamic information about glassy solids at various temperatures within the MD timescale^29,30,31,32. Ten cycles were employed for each DMA simulation. Moreover, atomic-scale dynamics in the DMA simulations were evaluated using atomic displacement Δr = |r(t)|. In all simulations mentioned above, the number of atoms was N_A = 30,000, and periodic boundary conditions were applied along all orthogonal directions. The probability density function p(Δr) was the histogram of Δr. The probability density function was normalized as ${\int }_{0}^{\infty }p\left(\Delta r\right)d\Delta r=1$. The non-Gaussian parameter α₂ was calculated by $3\left\langle \Delta {r}^{4}\right\rangle /5{\left\langle \Delta {r}^{2}\right\rangle }^{2}-1$ (Ref. ³⁵). We used OVITO³⁶ for the visualization of atomic models.

Results

Figure 1 shows the calculated α-relaxation time τ_α during cooling from 2000 to 1000 K. As the temperature decreased, τ_α drastically increased owing to the increase in liquid viscosity. For each alloy composition, the data were fitted by the Vogel-Fulcher-Tammann (VFT) equation:

$${{{\tau }}}_{{{\alpha }}}={{{\tau }}}_{0}{\rm{exp }}[{D}_{0}{T}_{0}/(T-{T}_{0})]$$

(1)

where τ₀, D₀, and T₀ are fitting parameters. The fragility parameter m, calculated by:

$${m={D}_{0}{T}_{0}{T}_{{\rm{g}}}/\{({T}_{{\rm{g}}}-{T}_{0})}^{2}{\rm{ln}}10\}$$

(2)

Fig. 1: Changes in the α-relaxation time of four liquid alloys.

The solid curves represent the VFT fitted results. Here, ${T}_{{\rm{g}}}^{* }$ is defined as the temperature at which the VFT equation satisfies τ_α = 10^–8 s (ref. ¹⁰). ${T}_{{\rm{g}}}^{* }$ = 920, 913, 949, and 966 K for Cu₃₀Zr₇₀, Cu₄₀Zr₆₀, Cu₅₂Zr₄₈, and Cu₅₈Zr₄₂, respectively.

characterizes the growth trend of relaxation time with decreasing temperature^37,38,39. A larger m means a more fragile (less strong) liquid, and vice versa^1,40,41. In general, the relaxation time at T_g is ~10² s in the laboratory⁴². To evaluate m, we used ${T}_{{\rm{g}}}^{* }$ instead of T_g. ${T}_{{\rm{g}}}^{* }$ and T_g are defined as the temperatures at which the VFT equation satisfies τ_α = 10^–8 s (Ref. ¹⁰) and 10² s, respectively. From Fig. 1, the estimated values of m were 42, 46, 70, and 84 at 30, 40, 52, and 58 at% Cu, respectively, indicating that fragility increased with increasing Cu content in the present model system. For comparison, we also calculated T_g and used it to evaluate m instead (see Fig. S1 in Supplementary Information). The m values obtained from T_g and ${T}_{{\rm{g}}}^{* }$ exhibited the same qualitative composition dependence. We used ${T}_{{\rm{g}}}^{* }$ in the later data for temperature normalization. Notably, the GFA of a realistic Cu-Zr alloy could not be targeted in this study. The GFA generally depends on the liquid properties, such as the fragility and thermodynamic stability of the corresponding crystalline phase. The experimental study also suggested that the density and atomic diffusivity of the amorphous phase have a significant correlation with the GFA in the Cu-Zr binary system, where Cu₅₀Zr₅₀, Cu₅₆Zr₄₄, and Cu₆₄Zr₃₆ have high experimental GFAs^43,44. For these reasons, the discussion of GFA was beyond the scope of the present study. Here, we revealed the correlations among liquid fragility, dynamic heterogeneity in liquid and glassy solids, and structural (i.e., static) heterogeneity in glassy solids.

Figure 2a shows the probability density function (i.e., histograms) of the atomic displacement Δr obtained by DMA simulations at various temperatures from 300 to 1300 K. Here, Δr was calculated based on the simulation snapshots, which were outputted at the end of each period (simulation frequency: 1 GHz). All models displayed a peak in the histograms at low temperatures. Upon increasing the temperature, both the whole histogram and the peak gradually shifted to the right, indicating increased atomic displacement by heating. At temperatures above 1100 K, all models entered the supercooled liquid state. Therefore, the figures showed changes in the histograms of atomic mobility in the unfreezing process (i.e., the inverse of the glass transition). In the medium-temperature range (around the unfreezing temperature), a shoulder appeared in the histograms on the right side of the peak. This shoulder was more apparent in the most fragile model (Cu₅₈Zr₄₂), resulting in the presence of two clear peaks at ${T}_{{\rm{g}}}^{* }$/T = 0.92 for this model. In the unfreezing process, the atoms could be classified as unfrozen atoms with large mobility and still-frozen ones with small mobility⁴⁵. Because the shape of the histograms reflected atomic mobility during unfreezing dynamics, Fig. 2a suggested distinct differences in the unfreezing behaviors of the four alloys. Cu₃₀Zr₇₀ (the strongest alloy) showed the most homogeneous unfreezing behavior, whereas Cu₅₈Zr₄₂ (the most fragile alloy) showed the most heterogeneous (or inhomogeneous) unfreezing behavior. Structural heterogeneity in metallic glasses has attracted much attention. Different types of local regions at the nanometer scale have been suggested to be present in the microstructure of metallic glasses including weakly bonded regions (WBRs; also called liquid-like regions or loosely packed regions) and strongly bonded regions (SBRs, also called densely packed regions)^{23,26,45,46,47}. Such structural heterogeneity is apparent in Fig. 2a. For instance, the fragile alloy exhibited two peaks in the medium-temperature region. The peaks on the left and right sides indicated characteristic atomic configurations with low and high atomic mobilities (correlated to densely and loosely packed regions), respectively.

Fig. 2: Dynamic heterogeneity obtained by DMA simulations.

a Histograms of atomic displacements during DMA simulations at T = 300–1300 K in increments of 50 K (${T}_{{\rm{g}}}^{* }$/T values are labeled on the plots). The frequency of the DMA simulation was 1 GHz. b Average atomic displacement Δr_ave and position of the first peak Δr_1st of the histograms obtained from Panel a. c Heterogeneity parameter χ (=Δr_ave/Δr_1st) calculated from b and non-Gaussian parameter α₂ (inset) calculated from a. The blue and red curves in a indicated the histograms corresponding to the peak temperatures of α₂ and χ in Panel c, respectively.

For the cases displaying one peak in Fig. 2a, such as in the low- or high-temperature regions, we defined Δr_1st as the peak position in the histogram. For cases displaying two peaks, such as those in the medium-temperature region, Δr_1st was defined as the left peak position, and it was an indicator of atomic displacement in regions with low atomic mobility. We also defined Δr_ave as the average atomic displacement of all atoms in the model by taking the statistical average of Δr. Figure 2b plots the Δr_1st and Δr_ave values from DMA simulations. In all the models, Δr_ave was larger than Δr_1st over the entire temperature range. Both Δr_1st and Δr_ave increased as the temperature increased. Compared to Δr_1st, Δr_ave displayed a drastic increase at a relatively lower temperature. Because this drastic increment in DMA simulations meant unfreezing from a glassy solid to a supercooled liquid, there was a significant time lag of unfreezing points within the same alloy model, indicating heterogeneous unfreezing dynamics. Using Δr_1st and Δr_ave, the strength parameter of heterogeneity in the unfreezing dynamics was defined as χ = Δr_ave/Δr_1st, and its values are summarized in Fig. 2c. A large χ indicated significant heterogeneity in the unfreezing process, and the largest value occurred at a temperature with the most heterogeneous unfreezing dynamics. In Fig. 2c, the more fragile models have larger χ values, suggesting that they exhibited significant heterogeneity in their unfreezing behaviors. The inset in Fig. 2c displays the non-Gaussian parameter α₂ calculated from the histograms in Fig. 2a⁴⁸. Because α₂ describes the deviation of a histogram from a Gaussian distribution, it provides information on the heterogeneity in atomic displacement⁴⁹. From the calculated α₂ values, the more fragile models also exhibited significant heterogeneity in their unfreezing behaviors. Notably, χ and α₂ had peaks at different temperatures, which are indicated by the red and blue curves in the histogram of Fig. 2a, respectively.

To obtain information on the SRO, we evaluated the Voronoi polyhedron of each atom in the atomic configuration from the DMA simulations. Each Voronoi polyhedron can be classified into six types (i = 1–6)⁵⁰. Types 1 and 6 are called Z clusters and non-Kasper clusters, respectively. The Z clusters with small disclination density are densely packed and have stable local structures, whereas the non-Kasper clusters with large disclination density are loosely packed and have relatively unstable local structures⁵⁰. Icosahedra are a well-known local order with a high packing density, and they are included in type 1. For the detailed classification of types 1–6, see Fig. 3a and Table S1 in the Supplementary Information. Figure 3b–e show clustering in the atomic system consisting of only type i atoms, as well as the average and histograms of the atomic displacement for each local order. First, the type of local order of atom α is represented by type (α). Then, using the function,

$${\delta }_{{ij}}=\left\{\begin{array}{l}1\,\left(i=j\right)\\ 0\,\left(i\,\ne\, j\right)\end{array}\right.$$

(3)

Fig. 3: Clustering, fraction, and atomic displacement for each type of local order.

Here, the atoms were classified into six types of local order (1–6). a Typical atomic structures of types 1–6. The center atom and neighboring atoms are colored brown and green, respectively. b Two-dimensional schematic diagram of the local order distribution and the definition of θ (the degree of clustering between the same local orders). c Relationship between the average clustering parameter $\left\langle {n}_{i}\right\rangle$ and the fraction of each type-i atom f_i. d Average Δr for each type of local order. e Histograms of atomic displacements in the four models. Insets: the same data plotted on log-log scale. In e, DMA simulations were conducted at T_χ, corresponding to the peak temperature χ in Fig. 2c. We used OVITO³⁶ for the visualization of atomic models in Figs. 3, 5, and 7.

The number of atoms belonging to the type i local order (i = 1–6) in the model is calculated by:

$${N}_{i}=\mathop{\sum }\limits_{\alpha =1}^{{N}_{{\rm{A}}}}{\delta }_{{\rm{type}}(\alpha )i}$$

(4)

Here, N_A is the total number of atoms in the model. The fraction of type i atoms in the model is defined by f_i = N_i/N_A. For an atom α of type i, the number of neighboring atoms of the same type i is represented by θ and calculated by:

$$\theta =\mathop{\sum }\limits_{\begin{array}{c}\beta =1\,\left(\ne \alpha \right)\\ {r}_{\alpha \beta } < {r}_{{\rm{c}}}\end{array}}^{{N}_{{\rm{A}}}}{\delta }_{{\rm{type}}(\alpha ){\rm{type}}(\beta )}$$

(5)

where r_αβ and r_c are the distance of atoms α–β and the cutoff distance, respectively (see Fig. 3b). Here, r_c is set to 0.4 nm and almost equal to the distance at the minimum point between the first and second peaks of the radial distribution function (RDF) (see Fig. S2 in Supplementary Information). Then, for type i atoms, the average θ value of the model is represented by $\left\langle {n}_{i}\right\rangle$ and calculated by the following equation:

$$\left\langle {n}_{i}\right\rangle =\frac{1}{{N}_{i}}\mathop{\sum }\limits_{\alpha =1}^{{N}_{{\rm{A}}}}\left\{{\delta }_{{\rm{type}}(\alpha )i}\mathop{\sum }\limits_{\begin{array}{c}\beta =1\,\left(\ne \alpha \right)\\ {r}_{\alpha \beta } < {r}_{{\rm{c}}}\end{array}}^{{N}_{{\rm{A}}}}{\delta }_{{\rm{type}}(\beta )i}\right\}$$

(6)

Hereafter, we used θ and $\left\langle {n}_{i}\right\rangle$ (named the clustering parameter and average clustering parameter, respectively) to evaluate the clustering tendency between the same type of local orders: a large θ or $\left\langle {n}_{i}\right\rangle$ value means that the local order forms clusters. Figure S3 in the Supplementary Information shows the spatial distributions of atoms belonging to type i in the Cu₃₀Zr₇₀ and Cu₅₈Zr₄₂ models at 300 K. For these models at 300 K, Fig. 3c shows the relationships between $\left\langle {n}_{i}\right\rangle$ and f_i. The solid gray line indicates a linear fit to the Cu₃₀Zr₇₀ plots. For both the Cu₃₀Zr₇₀ and Cu₅₈Zr₄₂ models, $\left\langle {n}_{i}\right\rangle$ increased monotonically with increasing average fraction f_i; this trend was intuitive from a statistical viewpoint. Meanwhile, order types 1 and 6 in the Cu₅₈Zr₄₂ model exhibited larger $\left\langle {n}_{i}\right\rangle$ than the solid gray line of Cu₃₀Zr₇₀, indicating a relatively larger $\left\langle {n}_{i}\right\rangle$ and a stronger tendency for clustering. We evaluated the average atomic displacements from the DMA simulations (simulation frequency: 1 GHz). The average displacement of atoms belonging to type i is shown in Fig. 3d. The atomic displacement is correlated to the type of order, with the types at low i numbers mostly showing smaller atomic displacements. In the DMA simulations, the atomic displacement reflects atomic bonding or packing density: a large (small) atomic displacement indicates weak (strong) atomic bonding or a low (high) packing density. Therefore, Fig. 3d suggested that types 1 and 2 may have mainly contribute to SBRs or densely packed regions, while types 5 and 6 may have mainly contributed to WBRs or loosely packed regions. Furthermore, Fig. 3d shows that Cu₅₈Zr₄₂ had a stronger correlation between the atomic displacement and the local order type than Cu₃₀Zr₇₀. In our previous study⁵¹, the icosahedra preferred to interpenetrate one other and formed medium-range orders (MRO). In the same study, we also reported that interpenetrating icosahedral MRO exhibited high local elastic stiffness. Here, two conclusions could be drawn based on Fig. 3c, d: (i) local atomic order with small and large disclination densities tended to cluster with the same type and induce structural heterogeneity in a fragile alloy, and (ii) the degree of structural heterogeneity and atomic mobility significantly depended on the alloy composition.

As shown in Fig. 2, the unfreezing behaviors exhibited heterogeneity depending on the alloy composition. Based on the classified types 1–6 mentioned above, we determined the origin of dynamic heterogeneity in the mechanical unfreezing process. For the temperature at which χ had a single maximum value, Fig. 3e shows histograms of the atomic displacement of each atom type i. In the case of Cu₅₈Zr₄₂ with heterogeneous unfreezing behaviors, we observed a significant difference in histogram shape among the various local orders. For instance, type 1 showed a high probability in the region with small Δr and a small probability in the region with large Δr, while type 6 exhibited the opposite trend. In contrast to Cu₅₈Zr₄₂, histograms of different local orders in the Cu₃₀Zr₇₀ model showed only small differences in their shapes. In Fig. 3e, the fragile alloys (i.e., Cu₅₈Zr₄₂) exhibited a large change in histogram shape between different local orders, suggesting a structural origin of their heterogeneous unfreezing dynamics seen in Fig. 2.

Figure 4a depicts the relationship between the average atomic displacement (i.e., average Δr or Δr_ave) and the clustering parameter (θ) for each local order during mechanical unfreezing of the Cu₅₈Zr₄₂ model in DMA simulations at ${T}_{{\rm{g}}}^{* }$/T = 1.07, 1.38, and 3.22 (the corresponding data for the Cu₃₀Zr₇₀ model are plotted in Fig. S4a in Supplementary Information). At ${T}_{{\rm{g}}}^{* }$/T = 3.22, there are weak correlations between Δr_ave and θ. For atoms of type 6, their Δr_ave increases with increasing θ, while atoms of icosahedra, type 1, and type 2 exhibited the opposite trend. The correlation between θ and Δr_ave for icosahedra, type 1, and type 6 became more significant with increasing temperature, reaching a strong correlation at ${T}_{{\rm{g}}}^{* }$/T = 1.07. This result provided important insights into the structural heterogeneity of glassy solids. In addition to features of the local order itself, the clustering of regions with the same local order correlated with dynamic heterogeneity during the mechanical unfreezing process. Notably, the clustering was correlated with atomic displacement in the cases of local orders with both large and small disclination densities.

Fig. 4: Average atomic displacement vs. clustering parameter and changes in the fraction and atomic displacement of each local order type in the Cu₅₈Zr₄₂ model obtained from DMA simulations.

a Relationship between the average atomic displacement and the clustering parameter θ. The values of ${T}_{{\rm{g}}}^{* }$/T are denoted on each figure (for Cu₃₀Zr₇₀, see Fig. S4a in Supplementary Information). Δr is the atomic displacement within a 1 ns interval. Only plots obtained from enough samples are shown. Panels in the far right also include two-dimensional schematics of clustering with θ = 4, 9, and 13. b Fraction and c normalized atomic displacement of atoms belonging to type 1 (Z) and type 6 (non-Kasper) clusters. d Difference in normalized atomic displacement between type 6 and type 1 clusters.

Figure 4b, c show the fraction of type 1 (Z) cluster and type 6 (non-Kasper) cluster atoms, as well as their average atomic displacements in the DMA simulations from low to high temperatures beyond ${T}_{{\rm{g}}}^{* }$. In Fig. 4c, we show the normalized atomic displacement of the type i local order defined by:

$${\Delta r}_{{\rm{n}}}^{i}=\left({N}_{{\rm{A}}}/{N}_{i}\right)\left(\mathop{\sum }\limits_{\alpha =1}^{{N}_{{\rm{A}}}}{\Delta {r}_{\alpha }\delta }_{{\rm{type}}(\alpha )i}/\mathop{\sum }\limits_{\alpha =1}^{{N}_{{\rm{A}}}}\Delta {r}_{\alpha }\right)$$

(7)

where $\Delta {r}_{\alpha }$ is the displacement of atom α. Figure 4d plots the difference in normalized atomic displacement between the type 1 (Z) clusters and type 6 (non-Kasper) clusters, ${\Delta r}_{{\rm{n}}}^{6}-{\Delta r}_{{\rm{n}}}^{1}$. In low-temperature glassy solids, the more fragile alloy had a larger fraction of type 1 (Z) clusters and a smaller fraction of type 6 (non-Kasper) clusters than the stronger alloy (Fig. 4b). In Fig. 4c, d, the normalized atomic displacement of type 1 (Z) clusters decreased, and that of type 6 (non-Kasper) clusters increased with increasing temperature in the case of fragile alloy, which suggested that dynamical heterogeneity in atomic mobility increased as the temperature approached the unfreezing point. Meanwhile, the Cu₃₀Zr₇₀ model did not show a significant increase in the heterogeneity of atomic mobility with increasing temperature. In all the alloys, the difference in atomic displacement between type 1 (Z) clusters and type 6 (non-Kasper) clusters became negligibly small at ${T}_{{\rm{g}}}^{* }$/T < 1 because the entire system became a supercooled liquid at this temperature.

Discussion

Figures 1–4 suggest the correlation among liquid fragility and structural and dynamic heterogeneities in glassy solids. Here, we revisited the glass-forming quenching process and discussed the origin of this correlation based on the development of local orders in the glass-forming quenching process. We evaluated the atomic displacement and the clustering parameter θ of Cu₅₈Zr₄₂ in the freezing process from liquid and summarized them in Fig. 5a (for the Cu₃₀Zr₇₀ model, see Fig. S4b in Supplementary Information). In high-temperature liquids, there was no clear difference in atomic displacement between different types of local order or any significant correlation between θ and the average displacement. This suggested that high-temperature liquids displayed neither dynamic heterogeneity nor significant structural features. However, the correlation between θ and the average displacement became significant upon lowering the temperature, especially in the Cu₅₈Zr₄₂ model. Figure 5a implies that in the freezing process, the formation of types 1 and 2 local orders and the development of clustering within them were correlated with the local atomic mobility and dynamic heterogeneity. Notably, not only local orders with low atomic mobility (such as types 1 and 2) but also those with high atomic mobility (such as type 6) exhibited a clear correlation between the clustering parameter θ and average displacement. Therefore, clustering between local orders should be one of the causes for forming the heterogeneous glass matrix with SBRs and WBRs.

Fig. 5: Average atomic displacement vs. clustering parameter and changes in the fraction and atomic displacement of each local order type in the Cu₅₈Zr₄₂ model obtained during cooling from the liquid.

a Relationship between the average atomic displacement and the clustering parameter θ (for Cu₃₀Zr₇₀, see Fig. S4b in Supplementary Information). Δr is the atomic displacement within a 1 ps interval. The value of ${T}_{{\rm{g}}}^{* }$/T is denoted in the lower left corner of each figure. Only plots obtained from enough samples are shown. Insets: distribution of type 1 atoms colored by θ, omitting the atoms of types 2–6. b Fraction and c normalized atomic displacement of atoms belonging to type 1 (Z) and type 6 (non-Kasper) clusters. d Difference in normalized atomic displacement between type 6 and type 1 clusters.

The fractions of local atomic orders and their atomic displacements in the cooling process from the high-temperature liquid are shown in Fig. 5b, c. In Fig. 5b, upon cooling from the high-temperature liquid, the fraction of type 1 (i.e., Z) clusters monotonically increased, while that of type 6 (i.e., non-Kasper) clusters decreased. Notably, fragile models exhibited a significant increase in type 1 clusters and a decrease in type 6 clusters at ${T}_{{\rm{g}}}^{* }$/T > 0.8. Figure 5c shows the normalized atomic displacements for atoms belonging to type 1 and type 6 clusters. The atomic displacements were evaluated at intervals of 1 ps and normalized by the average atomic displacement of each model. Figure 5d shows the difference in the normalized atomic displacements of type 1 and type 6 clusters. In Fig. 5c, d, as the temperature decreased, the normalized atomic displacement of type 1 clusters decreased while that of type 6 clusters increased, indicating increased heterogeneity in the atomic displacement as the temperature approached ${T}_{{\rm{g}}}^{* }$ during cooling. It was interesting that the decrease in the normalized displacement of type 1 clusters at ${T}_{{\rm{g}}}^{* }$/T > 0.8 significantly depended on the alloy composition, with the fragile alloys exhibiting a more significant decrement. In Fig. 1, the relaxation time of the fragile alloy significantly increased (i.e., the atomic mobility significantly decreased) at ${T}_{{\rm{g}}}^{* }$/T > 0.8, implying that fragility was correlated with an increase in atomic orders with a densely packed state and a decrease in atomic mobility of the atomic orders. Figures 4d and 5d suggest two features: the increasing dynamic heterogeneity of fragile alloy as the temperature approached ${T}_{{\rm{g}}}^{* }$ from both the glassy solid on heating and the supercooled liquid on cooling and that dynamic heterogeneity in liquid and glassy solid was correlated with the local orders.

Recently, a theoretical model has been developed for the fragility of glass-forming melts^13,14,15,52. This model predicts that alloys with a larger thermal expansion coefficient or steeper repulsive interatomic interaction have larger fragility, as expressed in the following equation,

$$m={V}_{{\rm{C}}}{C}_{{\rm{G}}}\left[1+\left(2+\lambda \right){\alpha }_{{\rm{T}}}{T}_{{\rm{g}}}\right]/\left({k}_{{\rm{B}}}{T}_{{\rm{g}}}{\rm{ln}}10\right)$$

(8)

where V_C, C_G, and k_B are the characteristic atomic volume, shear modulus at T_g, and Boltzmann constant, respectively¹³. The parameter λ represents the steepness of the repulsive interatomic interaction, and α_T is the thermal expansion coefficient calculated by α_T = (∂V/∂T)/V. Equation (8) is derived from the Krausser-Samwer-Zaccone (KSZ) model, which can represent T-dependent viscosity as a function of V_C, C_G, k_B, λ, α_T, and T¹³. To evaluate the temperature-dependent shear modulus, the KSZ model uses the Born-Huang ansatz⁵³. Previous studies based on the KSZ equation mainly focused on the repulsive interatomic interaction and discussed the fragility of various systems^13,14,15. Although the present study focused instead on the correlation between liquid fragility and spatial heterogeneity in glass solids using four alloy models constructed by the same interatomic interaction, Eq. (8) theoretically supported the findings of this study. To discuss a connection between Fig. 1 (i.e., fragility) and Figs. 2–5 (i.e., local orders, clustering, and atomic mobility) based on Eq. (8), we calculated the RDF at ${T}_{{\rm{g}}}^{* }$ and the temperature-dependent thermal expansion coefficient α_T. The parameter λ was the steepness of the repulsive side (i.e., the left-hand side) of the first peak of the RDF. Figure 6a represents the RDF at ${T}_{{\rm{g}}}^{* }$, and it suggested that Cu₅₈Zr₄₂ had the largest steepness and Cu₃₀Zr₇₀ had the smallest steepness (i.e., λ) among the four models. Equation (8) predicted that alloys with large λ had large fragility, which agreed with the present results. In addition, Fig. 6b shows the change in the RDF of the Cu₃₀Zr₇₀ and Cu₅₈Zr₄₂ models between ${T}_{{\rm{g}}}^{* }$/T = 0.7 and 1.0 during the MQ processes. The heights of the first and second peaks changed significantly in the Cu₅₈Zr₄₂ model compared to the Cu₃₀Zr₇₀ model, suggesting that the atomic packing density significantly increased as ${T}_{{\rm{g}}}^{* }$/T increased from 0.7 to 1.0 during the cooling processes. This trend was confirmed in Fig. 6c: a significant increase in the thermal expansion coefficient of models with high fragility was observed when the temperature decreased from ${T}_{{\rm{g}}}^{* }$/T = 0.7. Equation (8) predicts that a large α_T induced a large fragility, which agrees with the present result. The KSZ equation was derived from the relationship between viscosity and high-frequency shear modulus. An early experimental study also indicated a correlation between the elastic constant and fragility of glass-forming liquids⁵⁴. The shear modulus is related to the atomic connectivity, Z, which is affected by the thermal expansion behavior. From a theoretical viewpoint, the temperature-dependent change in atomic connectivity increases with increasing λ or α_T, as denoted by δZ ~ (−ϕα_TδΤ)^1+λ when δΤ < 0 (ϕ is the atomic packing fraction)13. In the present system, alloys with large λ had large α_T as seen in Fig. 6a–c, and thus they had a large atomic connectivity. In addition, Fig. 6c predicts that the atomic connectivity in a fragile alloy significantly increases when T approaches ${T}_{{\rm{g}}}^{* }$ in the cooling process. The increment of the clustering trend or atomic connectivity reduced atomic mobility, as seen in Figs. 4 and 5, and increases viscosity. Figure 6d schematically depicts the factors affecting fragility as observed in this alloy system. During the cooling process, local orders (i.e., SRO and MRO) with high packing density appeared. The development of local order in supercooled liquids was more significant in fragile systems such as Cu₅₈Zr₄₂. The local order increased packing density in the local region and the whole model. As a result, the thermal expansion coefficient of the whole model increased, as shown in Fig. 6c. The development of local order was also attributed to clustering or atomic connectivity. The areas of local order with high packing density and high atomic connectivity in supercooled liquid reduced atomic diffusivity (or mobility), increased viscosity, and affected fragility. We also obtained the thermal expansion coefficient at ${T}_{{\rm{g}}}^{* }$ from the volume vs. temperature relation in the heating process from the glassy solid. Figure 6e represents volume changes in the heating process and the obtained relationship between m and α_T${T}_{{\rm{g}}}^{* }$. The fragility m increased with increasing α_T${T}_{{\rm{g}}}^{* }$, as predicted by Eq. (8). This study can add a new aspect to the discussions related to Eq. (8): The development of local order (and thus heterogeneity) in supercooled liquids affects the thermal expansion coefficient, atomic connectivity, atomic mobility, and fragility. The development of local order in a supercooled liquid also induces structural heterogeneity in glassy solids after vitrification. Therefore, there is a correlation between the liquid fragility and structural heterogeneity of the glassy state in the present glass-forming systems. This discussion should also be relevant to other glass-forming metallic systems when there are areas of local order with low atomic mobility in the supercooled liquid and related structural heterogeneity in the glassy solids.

Fig. 6: RDF and thermal expansion coefficient.

a RDF of each model at ${T}_{{\rm{g}}}^{* }$ obtained from snapshots during the MQ process. b Changes in the RDF of the Cu₃₀Zr₇₀ (upper) and Cu₅₈Zr₄₂ (lower) models between ${T}_{{\rm{g}}}^{* }$/T = 0.7 and 1.0 during the MQ processes. c Changes in model volume during the MQ processes. The volume V is normalized by the volume at 0 K (denoted by V₀). The inset figure shows changes in the thermal expansion coefficient obtained from the V vs. T curves. d Schematic diagrams of changes in local order, clustering (connectivity), and atomic mobility of the present alloy systems. e Changes in model volume during the heating process and the relationship between m and α_T${T}_{{\rm{g}}}^{* }$ (inset).

The spatial distributions of atoms with large and small atomic displacements are shown in Fig. 7a at ~${T}_{{\rm{g}}}^{* }$ in the DMA simulations. Here, we only show alloys with the smallest and largest fragilities (i.e., Cu₃₀Zr₇₀ and Cu₅₈Zr₄₂). In the figure, regions with large and small atomic displacements are distributed on the nanometer scale, revealing nanoscale heterogeneity in the unfreezing dynamics. The heterogeneity in atomic displacement was more obvious in Cu₅₈Zr₄₂ than in Cu₃₀Zr₇₀. This agreed well with the results shown in Figs. 2–4. In addition, the present structural and dynamic analyses provided interesting insights into structural heterogeneity. The structural heterogeneity of metallic glasses includes aspects such as (1) the volume fractions of loosely and densely packed regions and (2) the difference in the degree (or intensity) of the loosely and densely packed states. Because denser atomic packing restricts atomic displacement and vice versa, Figs. 2–4 provide information on structural heterogeneity in terms of both the volume fractions and the packing intensity (aspects (1) and (2) mentioned above) in each alloy.

Fig. 7: Features of structural heterogeneities.

a Snapshots showing the distributions of atomic displacements during DMA simulations of Cu₃₀Zr₇₀ (the strongest model) at ${T}_{{\rm{g}}}^{* }$/T = 0.98 and Cu₅₈Zr₄₂ (the most fragile model) at ${T}_{{\rm{g}}}^{* }$/T = 1.05. b Local orders, packing density, atomic mobility, and temperature (temp.) of freezing and unfreezing in the microstructure of metallic glasses. c Schematic illustrations of structural heterogeneity in glass alloys that are less heterogeneous (left) and more heterogeneous (right). The red and blue regions correspond to loosely packed regions with high atomic mobility and densely packed regions with low atomic mobility, respectively.

Figure 7b schematically compares the structure, packing density, atomic mobility, and freezing/unfreezing behaviors in the microstructure between the SBR and WBR suggested by Figs. 2–5. In the heterogeneous microstructure with SBRs and WBRs, types 1 and 2 orders with high clustering parameter θ are located around the inner region of the SBR, while types 5 and 6 orders with high θ are located around the inner region of the WBR. In other words, the fractions (or θ) of types 1 and 2 local orders were larger in the SBR than the WBR. Therefore, the SBR had a high packing density and low atomic mobility. The structure, packing density, and atomic mobility should all gradually change from around the inner region of the SBR to around the inner region of the WBR, indicating a spatial gradient in the microstructure and properties of metallic glasses at the nanometer scale. In the MQ process, freezing occurs around the inner region of the SBR to the WBR, while mechanical unfreezing occurs in the opposite direction. In other words, freezing during the MQ process (or unfreezing during DMA) in SBR occurs at a higher temperature than in WBR. Ke et al. reported a critically reduced free volume value for the glass transition in various metallic glasses, providing an understanding of the glass transition in metallic glasses similar to Lindemann-type melting behavior¹⁹. The results in Figs. 3 and 4 reveal the effect of free volume (or packing density) on the glass transition and provide additional information on heterogeneous melting depending on the local free volume content (or packing density, SRO).

Figure 7c schematically depicts structural heterogeneity in terms of the volume fraction and degree of densely/loosely packed regions. The image on the right has more significant structural heterogeneity than that on the left. Our previous study²⁶ suggested that among the present four alloys, both the average excess free volume and average excess potential energy decreased with increasing Cu content, reaching the lowest values in the Cu₅₈Zr₄₂ model. The present study indicated that the alloy with smaller excess free volume and potential energy had larger densely packed regions and smaller loosely packed regions (Fig. 7c), as well as a more significant difference in packing density (or atomic displacement) between the loosely and densely packed regions. In other words, the fragile alloy had a large fraction of densely packed and a small fraction of loosely packed regions, as well as a significant difference in the degree between the densely and loosely packed states in this model system. During the cooling process, the increased fraction of Z clusters with low atomic mobility induced an increase in the liquid viscosity, as shown in Fig. 1. Namely, in the present alloy system, the densely packed structure induced both high fragilities in the liquid and high structural heterogeneities in the glassy solid. Therefore, there was a correlation among liquid fragility, dynamic heterogeneity in the freezing and unfreezing behaviors, and structural heterogeneity in glassy solids. The glassy solids that had significant structural heterogeneity induced by large differences between densely and loosely packed states exhibited considerable dynamic heterogeneity in their freezing and unfreezing behaviors.

Conclusions

Using binary alloy systems, we theoretically investigated the correlations among fragility, dynamic heterogeneity in the freezing process from liquid, and dynamic and structural heterogeneities in glassy solids. The DMA simulations indicated heterogeneous unfreezing behaviors in the glassy solids, with the most fragile alloy exhibiting the largest dynamic heterogeneity and vice versa. Heterogeneity in atomic displacements increased as the temperature approached the unfreezing point. We found a correlation between the clustering of local orders and atomic displacement during freezing from the liquid and unfreezing from the solid, suggesting one of the origins of the dynamic and structural heterogeneities. The alloy with the largest fragility contained the largest fraction of regions with small atomic displacement and showed the largest difference in atomic displacement between the densely and loosely packed orders, implying that the degree of packing density displayed the largest difference among the four models. The high fragility was affected by the development of clustering between the local orders with lower atomic mobility in the quenching process. These results were supported by a theoretical model for fragility, which depended on the effects of thermal expansion and atomic connectivity. This study indicated that the development of local order, which appeared in supercooled liquid and induced structural heterogeneity in glass solids, can affect fragility. Because the atomic orders and the clustering with low and high atomic mobility are found in the glass-forming quenching process as crucial factors for developing structural heterogeneity in glassy solids, the present alloy system exhibited correlations among liquid fragility and dynamical and structural heterogeneities in glassy solids.