Universal origin of glassy relaxation as recognized by configuration pattern matching Open Access

ABSTRACT

Relaxation processes are crucial for understanding the structural rearrangements of liquids and amorphous materials. However, the overarching principle that governs these processes across vastly different materials remains an open question. Substantial analysis has been carried out based on the motions of individual particles. Here, as an alternative, we propose viewing the global configuration as a single entity. We introduce a global order parameter, namely the inherent structure minimal displacement (IS D_min), to quantify the variability of configurations by a pattern-matching technique. Through atomic simulations of seven model glass-forming liquids, we unify the influences of temperature, pressure and perturbation time on the relaxation dissipation, via a scaling law between the mechanical damping factor and IS D_min. Fundamentally, this scaling reflects the curvature of the local potential energy landscape. Our findings uncover a universal origin of glassy relaxation and offer an alternative approach to studying disordered systems.

INTRODUCTION

Glasses are disordered materials that lack the structural long-range order of crystals but behave mechanically like solids. They are typically produced by rapidly cooling a viscous liquid to prevent crystallization. However, the dynamic processes through which liquids acquire amorphous rigidity during cooling remain incompletely understood [1–5].

One hallmark of glasses and glass-forming liquids is their complex relaxation dynamics, which span ∼12 orders of magnitude in timescales accessible with current technologies. These relaxation processes significantly affect the mechanical and functional properties of glassy materials [6–10]. Understanding how structural rearrangements govern these relaxation processes is crucial for uncovering the nature of glass and designing amorphous materials with improved properties [6–8,11,12].

Microscopic explanations of these dynamics are the fundamental components of glass transition research. For instance, we have demonstrated that the β relaxation in metallic glass could arise from string-like cooperative motions [13]. Berthier and co-workers have recently proposed that excess wings result from rearrangements of rare and localized regions over broadly distributed timescales [14]. Guan and colleagues have shown that low-temperature relaxation might be attributed to revisable atomic extrusions with hybrid features of vibration and diffusion motions [15]. Chang and colleagues found evidence that low-temperature relaxation is correlated with the liquid's light temperature properties [16]. More recently, Lukenheimer et al. made a correlation between the glass transition temperature and the thermal expansion coefficient with the necessary inclusion of the cooperativity parameter [17], the fragility.

These studies have enhanced our understanding of the nature of glass and glassy motions, providing valuable insights into the design of amorphous materials. However, a general theoretical framework for relaxation dissipation in disordered systems is still lacking. It remains unclear whether there is a unified mechanism for the various relaxation processes observed in diverse glasses [2,5].

Central to these relaxation measurements using various dynamic spectroscopic techniques is the loss angle δ (also known as the phase-lag angle or damping angle) between the input and output signals, which characterizes the energy dissipation during cyclic perturbation. In particular, δ is widely used as the damping factor of internal friction and has practical engineering applications. It is directly related to density and shear stress fluctuations through the fluctuation-dissipation theorem [18]. Thus, in this work, we focus on a unified understanding of δ.

Developing a general framework that unifies the relaxation dissipation of different glassy systems would be a significant step forward in understanding the nature of glasses and the glass transition. Such a framework could provide a basis for designing new amorphous materials with tailored mechanical and functional properties [19–21].

One approach that may be useful in addressing this issue is the potential energy landscape (PEL) [22–24]. The PEL refers to a potential energy function of an N-body system Ф(r₁, …, r_N), where the vectors r_i comprise position, orientation and intermolecular coordinates. It is a multidimensional landscape. For example, in the case of N identical atoms, the landscape is a (3N + 1)-dimensional object.

An essential aspect of the PEL is the presence of local minima, also known as inherent structures (ISs). It has been suggested that how a system explores its landscape as a function of temperature can provide insights into its dynamic behavior. However, due to the high dimensionality of the PEL, it is a challenging task to establish quantitative connections between the features of the PEL and material properties [25–30].

In this study, we propose a novel concept called ‘inherent structure minimal displacement’ that utilizes a pattern-matching technique to quantify the variability of IS configurations. We demonstrate that this global order parameter can unify relaxation dissipation into a scaling law, providing a general principle for understanding relaxation in glass-forming liquids.

RESULTS

Definition of IS D_min

We start by considering an instantaneous configuration at t = t₀, as shown schematically in Fig. 1a, where the atoms are labeled (1, 2, 3, 4, 5) (the blue disks). Later, at a time t = δt + t₀, these atoms move a certain distance, and a new configuration is arrived at, as indexed by (1', 2', 3', 4', 5') (the yellow disks) in Fig. 1. The displacement vectors are indicated by dashed lines for each atom.

Now, we consider the inherent structures of the two configurations and consider that all the atoms of the same chemistry are the same (identical particles), in other words, the atoms do not have labels. Then the minimal displacement between the two configurations (1, 2, 3, 4, 5) and (1', 2', 3', 4', 5') can be represented by the red solid arrows in Fig. 1 (referred to as IS D_min, inherent structure minimal displacement).

We note that the IS D_min reflects the minimal cooperative rearrangements, and it can be cast into a scalar value (see Equations (2) and (3) in Methods). It is smaller than the labeled displacement (each dashed vector in Fig. 1a). For instance, in cases where atomic motions involve precise position replacement and result in a loop-like motion, the IS D_min would be zero, while the labeled displacements would be a finite value.

Since experimental relaxation dynamics are probed by macroscopic tools, we expect that only the variations between the two configurations would have detectable effects on relaxation properties. Such a variation can be characterized by the IS D_min as outlined above. This is the major premise and the starting point of the present work. In the following, we show that such a principle indeed works reasonably well for a wide range of glass-forming liquids. In particular, a scaling law will be revealed.

We note that the root-mean-square displacement (RMSD), which is widely used in atomic dynamic analysis [1], is defined based on labeled displacements rather than IS D_min. On the other hand, the definition of IS D_min bears some similarity with the overlap parameter suggested by Parisi for spin systems [31], which has also been used in glass research [32]. Their differences are notable. The overlap parameter is a coarse-grained indicator for ‘similarity’ between two configurations. As will be discussed, IS D_min characterizes the shortest distance between two neighboring ISs at the same potential level.

The evaluation of IS D_min is not as straightforward as RMSD or the overlap parameter, because one has to find an optimum pattern matching between two configurations for the calculation of IS D_min. In this work, we use a linear sum assignment algorithm, which is also known as minimum weight matching in bipartite graphs and the ‘Hungarian Algorithm’ in computer science [33]. Computation details for IS D_min are presented in the Methods section. A simplified 1D model is given in Fig. 1e, which illustrates that the central idea is to find a transform matrix that operates on the configuration.

Figure 1b and d compare the labeled displacements (i.e. associated with RMSD) and the IS D_min based on the same configuration of an Al₈₅Sm₁₅ model glass-forming liquid over a waiting time of ∼2000 pico-seconds (ps) in a supercooled state [τ_α (T = 760 K) ∼2.5 × 10⁴ ps]. They are sliced from a 3D model, with a slice thickness of 10 Å. One can see that although the labeled displacements are large, the IS D_min is relatively small. This suggests that although the motions of individual atoms are significant, the global configuration is only slightly varied.

Figure 1c reinforces this observation by showing their distribution plot. It indicates that lots of atomic motions are position-replacing as suggested by the discrete peaks around 2.8 (the second peak) and 5 Å (the third) in the labeled displacements. As discussed previously (e.g. ref. [13]), the peaks and valleys in the displacement distributions p(u) suggest that the atomic movements proceed cooperatively, that is, one atom jumps to the position that was previously taken by another atom, such as its nearest or secondary neighbors. These can be seen from the fact that the peak positions of p(u) match those of g(r), the radial distribution function.

Scaling relation between δ and IS D_min in Al₈₅Sm₁₅

To substantiate this proposal, we use the molecular dynamics simulations of dynamic mechanical spectroscopy (MD-DMS, see Methods and Fig. S1) [13,34,35] to study the relaxation behaviors in a set of model glass-forming liquids. We start with a model system that has a composition of Al₈₅Sm₁₅ based on the force-field of ref. [36]. In a previous work, Sun et al. [37] illustrated that this force-field can yield relaxation processes that are consistent with the experimental DMS results and predict an additional process.

Meanwhile, the MD-DMS method uses the same protocol as in DMS experiments, and it has been applied in studies of several relaxation processes in glasses [13,15,34,35,37–39]. Thus, our simulations can be considered in-silico versions of DMS experiments.

As a typical example, Fig. 2 presents the relaxation properties as a function of time period t_w of the perturbation (which is inverse to the frequency f = 1/t_w) for temperature T = 740 K and pressure p = 0. It includes the storage shear modulus G' in Fig. 2a, loss shear modulus G'' in Fig. 2b, phase angle δ in Fig. 2c between the applied strain and response stress, as well as the IS D_min in Fig. 2d. The decay of G' (Fig. 2a) and the broad peak of G'' (Fig. 2b) are due to the so-called primary α relaxation. It is the dynamic process of the glass transition. For smaller t_w ≪ τ_α(T = 740 K) = 3 × 10⁵ ps (the peak time of G'') the system behaves like a solid, having a small value of δ→0, while the system is liquid-like (δ→π/2) for larger t_w ≫ τ_α.

Both δ and IS D_min have limiting value scopes: δ∈(0, π/2) where the low and high limited values correspond to the Hook elasticity (σ ∝ γ) and Newtonian fluid (σ ∝ dγ/dt) respectively; IS D_min ∈ (0,∼r₀/2) where r₀ is the average distance between the atoms for systems with isotropic interactions (e.g. without bonding).

Figure 3 shows 2D contour plots of the relaxation properties over the temperature range from T = 720 to 800 K at p = 0. All the data are from the equilibrium supercooled liquid states (Figs S2 and S3). We have examined a wide range of t_w, covering 5 orders of magnitude in timescales, from 100 ps to 2 × 10⁶ ps (i.e. 2 microseconds). We did not include the short-time data (t_w < 100 ps), as they might be influenced by the contributions from atomic vibrations.

We note that the wide time scale is essential for relating the simulation results to the experiments, as demonstrated in recent works [13,37].

One can see from Fig. 3a and b that the α relaxation moves to higher t_w when the temperature is lowered, suggesting the sluggish dynamics in the supercooled liquids. These features mimic the experimental observations. It is interesting to find that the phase angle δ in Fig. 3c and the IS D_min in Fig. 3d show a similar dependence on temperature T and period t_w.

Figure 4a plots the phase angle δ as a function of the IS D_min. It reveals a clear correlation between the two and all the data collapse on a single master line. The least square fitting to data yields a power law relation:

where the constants a, A and b are fitting parameters and the power b = 2.8 ± 0.1. In addition, r₀ = (V/N)^1/3 is a typical length scale (V and N are the volume and number of particles respectively). It is introduced such that the parameter A is dimensionless. This relation holds individually for each temperature as well (Figs S4 and S5). In this work, we focus on the scaling behavior and the b values, and do not delve into the details of A. As will be discussed later, this relationship can be explained based on PEL.

In addition, we note that some other functions, such as the labeled (traditional) RMSD and the IS RMSD (i.e. the RMSD based on the inherent structures), are not correlated with phase angle (Fig. S6). Computing the direct D_min (for instantaneous configurations, without evoking the minimization to obtain IS), we find that the direct D_min and phase angle are also correlated, while the curve is non-linear in the log-log plot (Fig. S6). In a previous work [35], we found that in the low temperature glassy state (non-equilibrium), δ is related to the number of atoms that jump with a distance larger than half of the mean interatomic distance. Such a relation is not found to collapse the data for the supercooled liquid studied in this work (Fig. S7). This might be due to the fact that the atomic motions in supercooled liquids are more abundant and cooperative.

These results suggest a fundamental relation between the relaxation property and the feature of PEL. The results in Fig. 4 indicate that the phase angle is predominantly controlled by the configurations overlapping. On the other hand, the atomic motions as probed by the labeled RMSD contribute to the thermal fluctuation and diffusion of the system. But only those atomic motions that influence the IS could contribute to the experiment-probed phase angle.

Generalization to different glass-forming liquids

Besides varying temperature at a constant external pressure p = 0, we have changed thermodynamic states of the system (in the equilibrium supercooled liquids as well) by tuning the external pressure p. As a typical example, Fig. 4b shows that p is in the range of (0, 4) GPa, while keeping the temperature a constant T = 800 K. As can be seen in Fig. 4b, all the data points with different pressures (as indicated by color squares) are overlapped with those if the temperature is changed (as indicated by gray circle symbols, which are the same in Fig. 4a). The relationship between IS D_min and δ, i.e. Equation (1), fits these data within the same degree of accuracy. Other combinations of temperature and pressure give the same results.

We note that our findings are not limited to the Al₈₅Sm₁₅ model liquids; we have validated our findings in six different model glass-forming liquids, including three metallic models, Ni₈₀P₂₀, Pd₈₀Si₂₀ and Y₆₅Cu₃₅, based on many-body embedded atom method potential [40,41], two models based on the Lennard-Jones (LJ) force-field, which are the Kob-Andersen (KA) model [42] and the Weeks-Chandler-Andersen (WCA) model [43], and a silica (SiO₂) model based on the three-body Tersoff force-field [44]. Results for these six systems are shown in Fig. 5. The detailed data are presented in Figs S8–S13 and the supplementary data. They confirm Equation 1 with reasonable accuracy. There are slight variations of the fitting parameters A and b, which might originate from the different interaction force-fields and the organizations of PEL.

Physical origin for the power-law scaling

Now we are in the position to elucidate what the fundamental cause is for the power-law scaling between δ and IS D_min. One usually evoked approach is to study how the exponent changes with the physical dimension. Among our studied examples, the KA-LJ system can form glasses for both two and three dimensions, and we conducted additional computations for this purpose (Fig. S13). Interestingly, we find that b ∼2.0 is nearly the same for the 2D and 3D glasses (Fig. S13). This implies that the presented scaling is invariant over different dimensions, calling for a fundamental understanding.

One possible explanation might be given based on the feature of PEL. As schematically shown in Fig. 6, IS D_min characterizes the shortest distance between two ISs that are at the same potential level for equilibrium state as investigated in this work. Since dissipation damping is a thermodynamically activated process with an activation energy ΔE_a, then we might have ΔE_a ∝ (IS D_min)^2+Δ where Δ is a number that describes the deviation of the PEL from the purely harmonic case. This is our theoretical premise based on PEL. Thus, we get δ ∝ (IS D_min)^2+Δ and b = 2 + Δ in Equation (1). Consequently, the value of b in Equation (1) characterizes the local curvature of the inherent state in the PEL. The observed power-law scaling is now explained.

We fit the G''(t_w) data with the Kohlrausch-Williams-Watts equation (KWW) f(t) ∝ exp[(−t/τ)^βKWW] and obtain the stretching parameter β_KWW (Fig. S15) [45]. In the coupling model of Ngai, n_CM = 1- β_KWW describes the coupling degree for the relaxation dynamics [46]. It is interesting to find that there is a modest correlation between β_KWW and b (Fig. 6b). Specifically, this shows that the Al₈₅Sm₁₅ and Ni₈₀P₂₀ models exhibit large b values and small β_KWW. This might reflect the curvature of the basin of PEL in these systems, as shown in the insets of Fig. 6b. The results in Fig. 6b are consistent with the coupling model [46].

DISCUSSION

We have shown that by introducing a parameter IS D_min, which characterizes the variability of inherent structure in the PEL, one can interpret the mechanical phase angle δ in a simple-yet-quantitative scaling law. Overall, the above results suggest a unified picture about the origin of mechanical phase angle δ: no matter what the thermodynamic state of the system is, the value of δ is uniquely determined by the variability of inherent structures over the probing time t_w, as quantified by IS D_min. For example, at low temperature and shorter t_w, the value of δ is small, as the variation between inherent structures is also small. The amorphous structures (packing, short- and medium-range orders) would influence δ through a change in the population of unique inherent structures.

Since δ is central to DMS and relaxation dynamics, our results provide a fundamental basis for the application of DMS in relaxation dynamics investigations in glass-forming liquids. The emerging physical picture is that only the variability of inherent structure in PEL provides relaxation properties that can be probed by mechanical spectroscopy. Moreover, the IS D_min, as a newcomer to the theoretical toolbox, might also be utilized to investigate other issues (e.g. aging, crystallization, deformation and glass-forming ability) in complex systems [47–50].

Finally, we anticipate that our discoveries could be investigated through experimental means, such as colloid experiments that can monitor particle positions, or by examining the contrast between relaxations and diffusions, which are primarily influenced by IS D_min and RMSD, respectively.

METHODS

Simulations of model glass-forming liquids

An open-source molecular dynamics (MD) simulation code, LAMMPS, was used for all the simulations. As listed in Table 1, we have studied seven different glass-forming systems. All these systems contain N = 9088 atoms. We used an NPT ensemble (that is, constant number of atoms, pressure and temperature) for the equilibration of the studied configurations at different temperatures (T) and pressures (p). Particular attention has been paid to making sure that the simulations are long enough for the configuration to reach the desired state. Periodic boundary conditions were applied for all the simulations.

The ISs are used in this work for the calculation of IS D_min. They are obtained by performing an energy minimization of the system, by iteratively adjusting atom coordinates. At that point, the configuration will be at a local minimum on the potential energy surface and the influence from atomic vibrations will be removed. Specifically, we used the conjugate gradient (CG) algorithm as implemented in LAMMPS for the calculations. We have checked that other algorithms, such as ‘fire’ and ‘quickmin’, give the same results.

Computation details for IS D_min

Briefly, we define a cost matrix for a pair of atoms i and j, via Equation (2):

Let X be a boolean matrix where X [i, j] = 1 if row i is assigned to column j and otherwise X [i, j] = 0. Then the IS D_min can be calculated by optimizing matrix X for the minimization of the sum of the product of C and X, as shown below in Equation (3).

In implementing the calculation of Equation (3), we used the linear_sum_assignment function from the scipy package for Python. It takes the cost matrix C_{i, j} as one of the input parameters and returns the optimum assignments that minimize the sum of the product of C and X in Equation (3) (i.e. one does not need to construct X by hand, more details can be found in the code). Obviously, if the optimum X is the unit matrix, then Equation (3) is reduced to the common RMSD. This condition occurs when all atoms have small displacements (Δr_i ≪ r₀/2). We note that the memory cost for the computation of IS D_minis at the order O(N²).

Dynamical mechanical spectroscopy and dissipation loss

The MD-DMS was performed on equilibrated configurations, covering a wide temperature range for the supercooled liquid. Specifically, at each configuration, we applied a sinusoidal shear strain

with an oscillation period t_w (related to frequency f = 1/t_w) and a strain amplitude γ_A, along the x-y direction of the metallic glass-forming liquid model. The resulting shear stress σ(t) was computed and fitted according to

wherein δ is the phase angle between stress and strain. The σ₀ is a linear term and is usually small. One typical example of the responding stress is shown in Fig. S1, together with the fitted line by Equation (5).

The strain was applied by a smooth change of the box shape, whereas stresses were directly measured from each component of the pressure. From these values, the storage (G′) and loss (G″) moduli are calculated as below:

A strain amplitude γ_A = 1.2% was applied in all the MD-DMS simulations. Such a value of strain amplitude ensures that the deformations do not change the structure and internal energy of the models. In other words, the deformation is in the linear response regime and the dissipation-fluctuation theorem is valid.

Technically, the NVT ensemble (constant number of atoms, volume and temperature) was used during the MD-DMS simulations. Ten cycles were applied for each MD-DMS to measure the storage and loss shear moduli as well as damping factor δ, fitted by Equation (6).

DATA AVAILABILITY

A Python script for the compunction of IS D_min is available in the supplementary data or at https://note.youdao.com/s/2RMcW3Oz.

ACKNOWLEDGEMENTS

The computational work was carried out on the public computing service platform provided by the Network and Computing Center of HUST. We are grateful to Prof. Ji-Chao Qiao (Northwestern Polytechnical University of China), Prof. Yang Sun (Xiamen University of China), Prof. Yun-Jiang Wang (Institute of Mechanics, Chinese Academy of Sciences) and Prof. Peng Tan (Fudan University of China) for advice, help and discussion.

FUNDING

This work was supported by the National Natural Science Foundation of China (52071147) and the National Thousand Young Talents Program of China. K.S. received support from Deutsche Forschungsgemeinschaft (DFG) via the Leibniz Program and SFB 1073.

AUTHOR CONTRIBUTIONS

H.-B.Y. and K.S. proposed the research; H.-B.Y. conducted the calculations and analyzed the data with input from L.G. and J.-Q.G. The manuscript was jointly written by H.-B.Y. and K.S.

Conflict of interest statement. None declared.

REFERENCES