Universal origin of glassy relaxation as recognized by configuration pattern matching

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 Dmin), 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 Dmin. 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][2][3][4][5].
One hallmark of glasses and glass-forming liquids is their complex relaxation dynamics, which span approximately 12 orders of magnitude in timescales accessible with current technologies.These relaxation processes significantly affect the mechanical and functional properties of glassy materials [6][7][8][9][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 of 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 spectroscopy 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][20][21].
One approach that may be useful in addressing this issue is the potential energy landscape (PEL) [22][23][24].The PEL refers to a potential energy function of an N-body system Ф(r1, ..., rN), where the vectors ri comprise position, orientation, and intermolecular coordinates.It is a multidimensional landscape, for example, in case of N identical atoms, the landscape is a (3N+ 1)-dimensional object.
Essential aspects of the PEL are the presence of local minima, also known as inherent structures (IS).It has been suggested that the way in which 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][26][27][28][29][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 Dmin.We start by considering an instantaneous configuration at t=t0, as shown schematically in Fig. 1(a), where the atoms are labeled as [1,2,3,4,5] (the blue disks).Later at a time t = δt+t0, these atoms move some distances, and a new configuration is arrived as index 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 the labels.Then the minimal displacement between the two configurations [1, 2, 3 ,4, 5] and [1', 2', 3', 4', 5'] can be represent by the red solid arrows in Fig. 1 (referred as IS Dmin, inherent structure minimal displacement).
We note that the IS Dmin reflects the minimal cooperative rearrangements, and it can be cast in to a scalar value (see Eqs.2 and 3 in method).It is smaller than the labeled displacement [each dashed vector in Fig. 1(a)].For example, when atomic motions are exactly position replacing and forming a loop-like motion, the IS Dmin 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 two configurations would have detectable effects on relaxation properties.Such a variation can be characterized by the IS Dmin 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, especially a scaling law will be revealed.
We note that the root-mean-square displacement (RMSD) which is widely used in the atomic dynamic analysis [1], is defined based on the labeled displacements rather than IS Dmin.On the other hand, the definition of IS Dmin 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 Dmin characterizes the shortest distance between two neighboring ISs at the same potential level.
The evaluation of IS Dmin 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 Dmin.In this work, we use a linear sum assignment algorithm, which is also known as minimum weight matching in bipartite graphs and the "The Hungarian Algorithm" in computer science [33].Computation details about IS Dmin are presented at the method section.A simplified one-dimensional (1D) model is given in Fig. 1(e), which illustrates that the central idea is to find a transform-matrix that operates on the configuration.As discussed previously (e..g, ref. [13]) that the peak and valleys in the displacement distributions p(u) suggest that the atomic movements proceed in a cooperative manner, that is, one atom jumps to the position that is 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 these of g(r), the radial distribution function.

Scaling relation between δ and IS Dmin in Al85Sm15.
To substantiate this proposal, we use the molecular dynamics simulations of dynamic mechanical spectroscopy (MD-DMS, see Method 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 Al85Sm15 based on the force-field of ref. [36].In a previous work, Sun et al [37] have 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 studying several relaxation processes in glasses [13,15,34,35,[37][38][39]. Thus, our simulations can be considered as an in-silico version of DMS experiments.Figure 3 shows 2-dimensional 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 (Fig. S2 and Fig. S3).We have examined a wide range of tw, covering 5 orders of magnitude in timescales, from 100 ps to 2×10 6 ps (i.e.,2 microseconds).We did not include the short-time data (tw < 100ps), as they might be influenced by the contributions from atomic vibrations.We note that the wide time-scale is essential for relating the simulations results to the experiments, as demonstrated in recent works [13,37].
One can see from Fig. 3(a) and (b) that the α relaxation moves to higher tw 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. 3(c) and the IS Dmin in Fig. 3 (d where the constants a, A and b are fitting parameters and the power b = 2.8 ± 0.1.In addition, r0 = (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 (Fig. 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 could be explained based on PEL.
In additional, 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 Dmin (for instantaneous configurations, without evoking the minimization to obtain IS), we find that the direct Dmin and phase angle are also correlated, while the curve is nonlinear 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 the 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 due to that the atomic motions in supercooled liquids are more abundant and cooperative.

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. 4 (b) shows the results that p in the range of [0, 4] GPa, while keeping the temperature a constant T = 800 K.As can be seen in Fig. 4(b), all the data points with different pressures (as indicated by color squares) are overlapped with those by changing the temperature [as indicated by gray circle symbols, which are the same in Fig. 4(a)].The relationship between IS Dmin and δ, i.e., Eq.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 Al85Sm15 model liquids; we have validated our findings in other six different model glass forming liquids, including three metallic models Ni80P20, Pd80Si20 and Y65Cu35 based on many-body embody atom method potential [40,41], two models based on Lennard-Jones force-field, which are the Kob-Andersen (KA) model [42] and the Weeks-Chandler-Andersen (WCA) model [43], and a silica (SiO2) model based on the three-body Tersoff force-field [44].Results for these six systems are shown in what the fundamental cause is for the power-law scaling between δ and IS Dmin.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 2 and 3 dimensions, and we conducted additional computations for this purpose (Fig. S13).Interestingly, we find that b ~2.0 is nearly the same for the 2 and 3 dimensional 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 Dmin characterizes the shortest distance between two ISs which are at the same potential level for equilibrium state as investigated in this work.Since dissipation damping is a thermodynamic activated process with an activation energy ΔEa.We might have ΔEa ∝ (IS Dmin) 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 Dmin) 2+Δ and b = 2+Δ in Eq.1.
Consequently, the value of b in Eq.1 characterizes the local curvature of the inherent state in the PEL.The observed power-law scaling is now explained.

Discussion
We have shown that by introducing a parameter IS Dmin, 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, these 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 tw, as quantified by IS Dmin.For example, at low temperature and shorter tw, the value of δ is small as the variation between inherent structures is also small.The amorphous structures (packing, shortand medium-range orders) would influence δ through the change 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 Dmin, as a newcomer to the theoretical toolbox, might also be utilized to the investigation of other issues (e.g., aging, crystallization, deformation and glass forming ability) in complex systems [47][48][49][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 Dmin and RMSD, respectively.

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 I, we have studied 7 different glass-forming systems.All these systems contain N = 9,088 atoms.We used a NPT ensemble (that is, constant number of atoms, pressure, and temperature) for the for the equilibration of the studied configurations at different temperatures (T) and pressures (p).Particular attention has been paid that the simulations are long enough to ensure the configuration reach the desired states.Periodic boundary conditions were applied for all the simulations.The inherent structures (ISs) are used in this work for the calculation of IS Dmin.
They are obtained by performing an energy minimization of the system, by iteratively adjusting atom coordinates.At that point the configuration will be in local potential energy minimum and the influence from atomic vibrations are removed.Specifically, we used the conjugate gradient (CG) algorithm as implemented in LAMMPS for the calculations.We have checked that other algorithm, such as 'fire' and 'quickmin' give the same results.Briefly, we define a cost matrix for a pair of atoms i and j, via Eq.2: Ci, j (t) = | ri (t+t0) -rj (t0) | 2 Eq.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 Dmin can be calculated by optimizing the matrix X for the minimization of the sum of the product of C and X, as shown below in Eq. 3.
Eq. 3 In implementing the calculation of Eq. 3, we used the linear_sum_assignment function from the scipy package for python.It takes the cost matrix Ci,j as one of the input parameters and returns the optimum assignments that minimize the sum of the product of C and X in Eq.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 Eq.3 is reduced to the common RMSD.This condition occurs when all atoms have small displacements (Δri << r0/2).We note that the memory cost for the computation of IS Dmin is at the order O(N 2 ).

Dynamical mechanical spectroscopy and dissipation loss.
The MD simulation of DMS (MD-DMS) was performed d on equilibrated configurations, covering a wide temperature range for the supercooled liquid.
Specifically, at each configuration, we applied a sinusoidal shear strain () = ( 2/ ) Eq. 4 with an oscillation period tw (related to frequency f = 1/ tw) 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

FIG. 1 .
FIG.1.Definition and quantification of IS Dmin.This schematic plot in (a) shows the labeled displacements (dashed vector) and the IS Dmin (solid vector) for each atom.Blue disks with label i (1, 2,.. ) indicate the initial position at t = t0, while yellow disks with label i' (1', 2',.. ) indicate the position at t = δt + t0.(b) Typical atomic displacements (black arrows) for an Al85Sm15 glass-forming liquid at T = 760 K and δt

Figure 1 (
Figure 1(b) and (d) compares the labeled displacements (i.e., associated with RMSD) and the IS Dmin for a same configuration of an Al85Sm15 model glass forming-liquids over a waiting-time about 2000 pico-seconds (ps) in a supercooled state [τα (T = 760 K) ~2.5×10 4 ps].They are sliced from a 3D model, with slice thickness 10 Å.One can see that although the labeled displacements are large, the IS Dmin is relatively small.This suggests that although the motions of individual atoms are significant, the global configuration is only slightly varied.

Figure 1 (
Figure 1(c) 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]) that the peak and valleys in the
) show a similar dependence on temperature T and period tw.

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

FIG4.
FIG4.Scaling relation between δ and IS Dmin.(a) data for configurations at different temperatures and the pressure is fixed at p = 0; (b) data for configurations at different p at a constant temperature T = 800 K as indicated by the color bar; the gray disk symbols are the same data from (a) p = 0.

FIG. 6 .
FIG. 6. PEL interpretation for the scaling relation.(a) Schematic PEL for dissipation relaxation.Relaxation damping is considered as a thermal activated process from two local PEL basins, which are separated by an energy barrier with the activation energy ΔEa.The green dashed line is a first approximation of the local PEL with a squared formula (i.e., the harmonic approximation).(b) Relation between βKWW and b for the studied systems; the insets are schematic illustration for the local PEL for systems with lower and higher b respectively.
(b) indicates a modest correlation between βKWW and b.Specifically, it shows that Al85Sm15 and Ni80P20 models exhibit large b values and small βKWW.This might reflect the local feature PEL in these systems, as shown in the insets of Fig.6(b).The results in Fig.6(b) are

TABLE 1 .
Models studied in this work.These include metallic glass-forming models based on the embedded atom method (EAM) potential, Kob-Andersen model based on Lennard-Jones potential (KA-LJ), Weeks-Chandler-Andersen model based on Lennard-Jones potential (WCA-LJ) and the three-body tersoff potential for SiO2.