Critical scaling of icosahedral medium-range order in CuZr metallic glass-forming liquids

The temperature evolution of icosahedral medium-range order formed by interpenetrating icosahedra in CuZr metallic glassforming liquids was investigated via molecular dynamics simulations. Scaling analysis based on percolation theory was employed, and it is found that the size distribution of clusters formed by the central atoms of icosahedra at various temperatures follows a very good scaling law with the cluster number density scaled by S−τ and the cluster size S scaled by |1 − Tc/T|−1/σ, respectively. Here Tc is scaling crossover-temperature. τ and σ are scaling exponents. The critical scaling behaviour suggests that there would be a structural phase transition manifested by percolation of locally favoured structures underlying the glass transition, if the liquid could be cooled slowly enough but without crystallization intervening. Furthermore, it is revealed that when icosahedral short-range order (ISRO) extends to medium-range length scale by connection, the atomic configurations of ISROs will be optimized from distorted ones towards more regular ones gradually, which significantly lowers the energies of ISROs and introduces geometric frustration simultaneously. Both factors make key impacts on the drastic dynamic slow-down of supercooled liquids. Our findings provide direct structure-property relationship for understanding the nature of glass transition.

The nature of glass and glass transition is the deepest and most interesting unsolved problem in the solid state science [1][2][3][4] . With a generic definition, numerous systems spanning a broad range of length scale such as atomic and colloidal systems, foams, and granular materials, can be considered as glass when certain conditions are satisfied 3 . Metallic glass, as a relatively "simple" glassy system for scientific research of glass transition and promising industrial material 5 , has attracted much attention and interest of scientists from broad research fields 6,7 .
The CuZr metallic glass-former has been extensively investigated [8][9][10][11][12][13][14][15][16] . It exhibits good glass-forming ability, and its binary chemical composition reduces the complexity of local atomic structures, which makes this system a good model for the study of structure-property relationships in liquids and glasses. Previous studies demonstrate that icosahedral short-range order (ISRO) is closely correlated with the slow dynamics and dynamical heterogeneity during glass formation [8][9][10][11][12][14][15][16] . It has also been demonstrated that the formation of icosahedral medium-range structures via the interpenetration of ISROs plays a key role in the dynamic slowdown 10,11 , The more the ISROs are connected, the slower the dynamics of the connected ISROs is. This implies that the icosahedral medium-range order formed by the percolation of ISROs may be related to the glass transition. So far, the concept of percolation of ISROs has been widely used for understanding the relationship between structural evolution and glass transition in metallic glass-forming liquids. However, the question is that no specific percolation theory or scaling analysis has been derived to establish a direct link between the percolation of ISROs and glass transition in the metallic glass-forming liquids. In addition, the population of ISROs in CuZr metallic glass-forming liquids sometimes is not very high so that ISROs do not even percolate as glass transition occurs 14 , For example, Cu 50 Zr 50 is a good glass-former in both experiments 13 and computer simulations 17,18 . However, the fraction of ISROs in Cu 50 Zr 50 at 300 K is found to be less than 4% in simulation 9 . Even in the inherent structure of Cu 50 Zr 50 , it is less than 15% 6 . To establish the link between percolation of ISROs and glass transition, the nearest-neighbor atoms of ISROs were also taken into account for the percolation. Meanwhile, the percolation of ISROs in CuZr metallic glass-forming liquids could be influenced by system size in simulation, leading to uncertainty for understanding the strcuture-dynamics relationships in these systems. Therefore, it is highly desirable to develop a scaling analysis based on percolation theory to establish a quantitative description of percolation of ISROs and link to glass transition. Moreover, although ISROs are found to correlate with the slow dynamics and the connection of ISROs makes it even slower, the physical origin is still not very clear.
In this work, molecular dynamics simulations were performed for Cu 50 Zr 50 with realistic interatomic potential (see Methods). We investigated the connection of ISROs via interpenetrating or volume-sharing. Graph theory was introduced to characterize the the clusters formed by the central atoms of interpenetrating ISROs at different temperatures as the CuZr metallic glass-forming liquids are cooled down, and the equilibrium cluster size distribution was analyzed. Scaling analysis based on percolation theory 19,20 was conducted for the cluster size distribution. It is found that the cluster size distributions at various temperatures collapse together and follow a good scaling law, as the cluster size is S scaled by |1 − T c /T| −1/σ and the cluster number density is scaled by S −τ , respectively. The scaling analysis suggests that there could exist a geometric phase transition of percolation of locally favoured structures, once the metallic glass-forming liquids are quenched slowly enough (in the hypothetical limit of infinitely long relaxation time) but without crystallization intervening, and glass transition may be related to the percolation of locally favoured structures. Furthermore, it is revealed that as ISROs are connected together, the atomic configurations of connected ISROs are optimized towards more regular icosahedra. The optimized ISROs enhance the geometric frustration, and the local energies are significantly lowered, which stabilizes the structure of liquids and slows down the dynamics in glass transition.

Results and Discussions
Scaling analysis for the percolation of icosahedral network. First, we analyzed the size distributions of clusters formed by the connection of central atoms of the icosahedra at various temperatures. In our analysis, graph theory was applied for the construction of connection of icosahedra and the resulting icosahedral network (see Methods). As shown in Fig. 1, the size distributions decrease monotonically as cluster size increases and exhibit similar behaviour at different temperatures. The size distributions in small size part follow power-law behaviour, but deviate in larger size part, which signals there is a finite characteristic cluster size in the system and, for a given temperature T, the characteristic cluster size marks the crossover between a power-law behaviour and a rapid decay of n(S, T), qualitatively. As temperature decreases, more and more larger clusters are formed in the supercooled liquids. It would be expected that the cluster size distribution could asymptotically approach a pure power-law behaviour as temperature further decreases and approaches some critical point if possible. The similarity of the clusters size distributions at different temperatures is reminiscent of a general scaling behaviour. Therefore, to explore the scaling law behind the cluster size distributions, a general scaling ansatz 19 for the cluster number density n(S, T) based on percolation theory was employed. That is, c where the characteristic cluster size S c diverges as a power-law in terms of the distance of T from T c : and the function f is known as the scaling function for the cluster number density. Generally, the expression of f varies from system to system, and dimension to dimension. Analytic solution of f is non-trivial in most cases. However, the asymptotic behaviour of f can provide sufficient information of the scaling behaviour of the size distributions. According to percolation theory, for S/S c ≪ 1, the scaling function is approximately constant, and for S/S c ≫ 1 it decays rapidly. To our knowledge 19 , except for one-dimension percolation, this behaviour of f is quite universal in the scaling analysis based on above scaling ansatz. To reveal the underlying scaling behaviour of the cluster size distributions, following above scaling ansatz, we scaled the cluster number density n(S, T) with S −τ and the cluster size S with |1 − T c /T| −1/σ , respectively. Figure 2 shows that the scaled cluster size distributions at different temperatures collapse onto a master curve representing the graph of scaling function f(S/S c ). Here T c = 700 K is the scaling crossover-temperature. τ = 2.05 and σ = 0.4 are the scaling exponents. These two exponents are different from the values obtained in three-dimensional site percolation [21][22][23] , because the ISROs percolate in a continuous space, but not on a discrete periodic lattice. As shown in Fig. 2, the scaling function f(S/S c ) decays rapidly as the scaled cluster size S T ≫ 1, which indicates that there is a characteristic cluster size in the system. The evolution of the characteristic cluster size S c shows a divergence behaviour as temperature is approaching the critical point T c . At T = T c , the characteristic cluster becomes infinite, so that n(S, because it can be seen from Fig. 2 that the scaling function f(S/S c ) will approach a non-zero constant for S T ≪ 1.
The scale-free behaviour of the cluster size distribution at threshold shows some information about the geometric properties of the percolating cluster, indicating that the incipient infinite cluster has an internal fractal geometry.
The fractal dimension d f ~2.86 of the structures formed by connected ISROs at the critical point can be calculated from the scaling relation of f where d = 3 is the spatial dimension 19 . It can be seen that the value of d f is not that small, which indicates that the atomic packing of incipient infinite cluster is not so loose. On the other hand, just as the characteristic cluster size S c diverges as T is approaching T c , the correlation length associated with the connected ISROs also diverges. For a particular characteristic cluster size S c , the associated radius of gyration defines a characteristic length scale that is proportional to the correlation length ξ. According to percolation theory, one has so that the temperature dependence (as T → T c ) of the correlation length ξ can be characterized as c which shows a power-law behaviour with a divergence at T = T c . The corresponding critical exponent ν = 0.875 can be determined by the scaling relation 19 of f As shown above, in percolation, once the cluster number density is known, all other quantities can be derived. All results indicate that the percolation of ISROs forming the so-called icosahedral medium-range orders indeed correlates with the glass transition in CuZr metallic glass-forming liquids. We argue that d f of the incipient infinite cluster emerging at T c cannot be calculated directly by the box-counting method. The reason is as follows: Since the crossover temperature T c is blow the glass transition temperature (over 900 K) obtained in the previous study 24 , the atomic structures at 700 K obtained in simulations are in non-equilibrium states and sensitive to the cooling-rate. As a result, d f obtained from box-counting analysis for the atomic structures at 700 K will be cooling-rate dependent. This is not exactly the same d f derived in percolation theory for the equilibrium structural phase transition. Therefore, d f obtained from the universal scaling relation (τ = d/d f + 1) is generic. It is also worth noting that T c is not cooling-rate dependent, because it is a critical temperature derived from the equilibrium structural phase transition of the modelling system, and all data for determining T c are generated from the relatively equilibrium supercooled liquid states. This is similar to T 0 in Vogel-Fulcher-Tamman (VFT) 2,3 equation.
Physical origin of dynamic slowdown associated with percolation. The above scaling analysis demonstrates that the percolation of ISROs in CuZr metallic glass former during cooling is closely related to glass transition. It is still not clear why the percolation of ISROs could contribute to slowing down of dynamics. As mentioned above, the connectivity of ISROs significantly influences the dynamical property of local structures 10 . In order to emphasize the importance of the medium-range structures formed by the connection of ISROs to the dynamic slowdown, we calculated the self-intermediate scattering functions 25 (SISFs) of the central atoms of icosahedra with different node degree . Figure 3 shows the SISFs as a function of time for central atoms of the icosahedra with different k values (due to the very rare connection of k > 5, SISFs of k > 5 were not calculated 10 ) in supercooled metallic glass-forming liquid at 1000 K. It can be seen that all SISFs for different k values exhibit non-exponential decay behaviour, and the decay becomes much slower as k value increases, dramatically depending on k. The inset in Fig. 3 explicitly shows that the scaled relaxation time τ k /τ α increases exponentially with the value of k/12 (τ α is the average structural relaxation time, and 12 is the maximum value of node degree) 10 . It can be seen that the relaxation time of the atoms with large k is even more than 10 times of the τ α . This finding reveals that the medium-range structures formed by the connection of ISROs fundamentally influences the relaxation dynamics of local atomic structures. As the icosahedral medium-range order percolates during quenching, the resulting dynamic slowdown spreads in the whole system, leading to the sluggish dynamics, which finally contributes to glass transition.
It is also very interesting to investigate how the icosahedral medium-range order influences the atomic symmetry of icosahedra themselves. All icosahedra in the MD modeled samples are actually distorted from the ideal icosahedron, because from a geometrical viewpoint, icosahedral clusters cannot fill the entire three-dimensional space without partially breaking of the five-fold rotational symmetry 26 . As observed by Frank over half a century ago, the ideal icosahedral arrangement indeed has a significantly lower energy than fcc atomic arrangement 27 . The question is whether the connection of ISROs will promote dense packing and five-fold local symmetry of ISROs. If so, the self-aggregation effect of icosahedra 9 will naturally tend to minimize the local energy density, slow the atomic dynamics, and lead to great geometric frustration. In order to get deep insight into the above discussion, we analyzed the local atomic symmetry of ISROs and its dependence on the connectivity degree k. To analyze the local atomic symmetry of ISROs, the bond orientational order (BOO) parameter introduced by Steinhardt et al. 28 was adopted, in which the BOO of the -fold symmetry is defined as a +  2 1 vector:    where the term in brackets in the invariants (11) is the Wigner 3j symbol 29 . For fcc, bcc, and hcp symmetry, the value of Ŵ 6 is close to zero, while for perfect icosahedron, = − . W 0 169757 6 , so that Ŵ 6 is sensitive to reflect the five-fold atomic symmetry features of ISROs in metallic glass-forming liquids. Figure 4 shows that the average values of Ŵ 6 of ISROs are not − 0.169757, but in between − 0.06 and − 0.10, suggesting that the ISROs in metallic glass-forming liquids are not perfect, but distorted with partially fcc symmetry, in agreement with the previous studies 26 . Furthermore, as shown in Fig. 4, with increasing node degree k, Ŵ 6 decreases to more negative values. This indicates that as the node degree increases, on the range we studied, the atomic configurations of ISROs are optimized towards more regular icosahedra. The configuration optimization of ISROs produces significant geometric frustration in the supercooled liquids.
We also calculated the formation energy E p,c of ISROs in the supercooled metallic liquids and investigated the dependence of the formation energy on node degree k. The formation energy of ISROs can be calculated by the formula of where E p,j is the potential energy of the jth atom in an icosahedron and E ref,a is the reference energy of the element a 30 . In our calculations, the chemical potentials of Cu and Zr in crystal structures (fcc for Cu and hcp for Zr) were used as the reference energies for Cu an Zr atoms, respectively. Therefore, in equation (13), the effect of the chemical compositions in an icosahedron was eliminated, and the formation energies of ISROs with different chemical compositions can be comparable. Thus, we can analyze the formation energies of ISROs with different degree k. The inset histogram in Fig. 4 clearly shows that with increasing degree k, the formation energy of ISROs decreases monotonically, indicating that the formation of the icosahedral medium-range order leads to lower energies, and the structures become more stable. All of above results indicate that the connection between different icosahedra will indeed generate a positive feedback for optimization of the local ISROs. This effect will minimize the local energy density, slow down the atomic dynamics, cause great geometric frustration, and finally contributes to glass transition.
In summary, we carried out MD simulations for CuZr metallic glass-forming liquids to investigate the relationship between percolation of ISROs and glass transition. As the system is supercooled, the ISROs tend to connect with each other to form big clusters. The cluster size distributions at different temperatures follows an excellent scaling law, which indicates that the cluster size distribution evolves toward a power-law behaviour, and an infinite size cluster is formed as the system is approaching a critical point. Therefore, there would be a structural phase transition manifested by percolation of locally favoured structures underlying the glass transition if the liquid could be cooled slowly enough but without crystallization intervening. Furthermore, the percolation of ISROs and the formation of medium-range orders optimize the atomic configurations of ISROs towards more regular icosahedra, enhancing the geometric frustration and minimizing the local energy density, which leads to the dynamic slowdown of the metallic glass-forming liquids. Our findings suggest that the geometric phase transition manifested by percolation of locally favoured structures could be critical for the understanding of the nature of glass transitions.

Methods
Molecular dynamic simulations. In this work, molecular dynamic (MD) simulations were carried out for Cu 50 Zr 50 metallic alloy using the LAMMPS package 31 . The interatomic interaction was described by the well-developed embedded-atom method potential for CuZr alloys 32 . The sample contains 10000 atoms being randomly distributed in a cubic box with periodic boundary condition applied in three dimensions. The MD step is 2 fs. At first, the sample was equilibrated at 2000 K for 4 ns (2,000,000 MD steps) in NPT (P = 0) ensemble with Nose-Hoover thermostat and barostat. The liquid was then quenched at a rate of 1 K/ps down to its target temperature (1600 K, 1500 K, 1400 K, 1300 K, 1200 K, 1100 K, and 1000 K, respectively). During cooling, the box size was adjusted to maintain zero pressure. At each temperature, the atomic configuration was relaxed in NPT (P = 0) ensemble for another 2 ns (1,000,000 MD steps) for the analysis of physical properties (500 atomic configurations were collected for ensemble average). The structural relaxation time (τ α ) of our modelling system at 1000 K was order of magnitude of 10 ps, so the ensemble average window for analysis was much longer than τ α of the system at each temperature we studied. Note that at temperatures higher than the glass transition point, the structures corresponding to (metastable) equilibrium can be achieved quite rapidly. Therefore, in general, the effect of cooling rate on the analysis of structural, dynamic, or other physical properties of our modelling system is eliminated.
Icosahedral network and node degree. The local atomic structures in supercooled liquid samples at different temperatures were analyzed by the Voronoi tessellation method [33][34][35] and identified in terms of the Voronoi index 〈 n 3 , n 4 , n 5 , n 6 〉 , where n i (i = 3, 4, 5, 6) denotes the number of i-edged faces of a Voronoi polyhedron. In our analysis, a cutoff distance of 5 Å was chosen so that the Voronoi index distribution was converged. To characterize the connectivity of ISROs in the system, we introduced the graph theory [36][37][38] . In our scheme 10,39 , the central atom (with Voronoi index 〈 0, 0, 12, 0〉 ) of an icosahedron is treated as a node, and two nodes are considered to be connected if they are the nearest neighbors with each other, that is, the two icosahedra are interpenetrating or volume-sharing. The choice of connection criterion is reasonable based on recent experimental observation 40 , in which the extent of icosahedral short-range order to form medium-range order is consistent with a facing-sharing or interpenetrating configurations. Therefore, the case of interpenetrating configuration was considered in our analysis. Property variation depending on the connection criterion will be discussed in detail in future works. With the definitions of nodes and edges abstracted from the atomic modelling system, we established the icosahedral network. Based on the scenario of graph theory, the node degree k was defined as the number of other nodes directly connected to it. The maximum value of node degree for the central atom of an icosahedron is 12. Our results suggest that node degree introduced from graph theory is a good nonlocal (to some extent) order parameter for classifying atoms with distinct properties.