INTRODUCTION
============

Atomic rearrangements are an essential property of glassy materials ([@R1], [@R2]); they are critical in controlling resistance to flow and are central to the evolution of many properties of glasses, such as their stability ([@R3], [@R4]) and mechanical performance ([@R5], [@R6]). Moreover, because of their disordered nature, glassy materials inherently feature diverse relaxation dynamics ([@R7]), ranging from picoseconds for fast processes to thousands of years for aging and densification. Understanding how the atomic rearrangements govern these processes constitutes an outstanding issue in glassy physics.

Particularly noteworthy is the so-called Johari-Goldstein (secondary or β) relaxation ([@R8]--[@R11]). Before the 1970s, the β relaxations were known in polymer glasses and assumed to originate from rotations of the side chains or functional groups of polymers. However, Johari and Goldstein ([@R8]) discovered β relaxations in several glasses of rigid molecules, suggesting that some β relaxations must involve the motion of the entire molecule instead of intramolecular degrees of freedom. They belong to a special class of relaxations having a strong connection to the structural α relaxation. To distinguish them from other secondary relaxations (mostly of intramolecular origin), Ngai and Paluch ([@R9], [@R10], [@R12]) introduced the name "Johari-Goldstein β relaxation."

The Johari-Goldstein β relaxations are intrinsic to almost all glassy materials ([@R10], [@R11]), and they are of significance to the understanding of many crucial unresolved issues in glassy physics and material sciences ([@R11], [@R13]). These include liquid dynamics and anomalies ([@R10], [@R14]), diffusion ([@R15]), physical aging ([@R16]), mechanical properties of polymeric and metallic glasses (MGs) ([@R17]--[@R21]), as well as the conductivity of ionic liquids and the stability of glassy pharmaceuticals and biomaterials ([@R22]--[@R24]). Consequently, understanding the structural origin of β relaxation is of practical importance for a wide variety of technologies. For example, polymers and MGs with pronounced β relaxations might exhibit good ductility and are favored for mechanical applications ([@R13], [@R18], [@R25]). On the other hand, for the pharmaceutical industry, the β relaxations are less desirable ([@R23]) and should be suppressed because they cause the recrystallization and destabilization of glassy medicines, which have better solubility and bioavailability compared with their crystalline counterparts.

Despite its importance, the structural origin of Johari-Goldstein β relaxation is still not established ([@R13], [@R26]). Some theoretical reflections on experimental or simulation results ([@R15], [@R27]--[@R33]) and a theoretical model ([@R34]) suggest a scenario of cooperative motions of atoms but without fully clarifying its nature, and direct evidence or convincing support is yet to be demonstrated. Other interpretations also exist ([@R35]--[@R38]).

For clarity, we note that in glass physics, there is another process called short-time (or fast) β relaxation occurring at picoseconds with a weak temperature dependence ([@R39]). We do not refer to this well-established process ([@R40], [@R41]). Here, all β relaxations refer to the Johari-Goldstein β relaxations.

In principle, molecular dynamics (MD) simulation is a powerful method to investigate atomic rearrangements of glasses at the microscopic level ([@R42], [@R43]). It would be helpful to reproduce the behavior of β relaxations in simulations and to scrutinize the details of structural rearrangements governing it ([@R13]). As demonstrated by experiments, the β relaxations can be discerned only when the observation time approaches (or exceeds) microseconds ([@R44], [@R45]), whereas most of the MD simulations conducted so far are within nanosecond time scales ([@R42], [@R43]), thus missing about three orders of magnitude. Moreover, in contrast to the case of structural α relaxations, we do not have a convenient time-correlation function to characterize the β relaxations for MD simulations \[except for a specially designed model of Fragiadakis and Roland ([@R46])\]. As a result, understanding the mechanism of β relaxation constitutes a challenging issue in glassy physics.

This work overcomes these difficulties by conducting MD simulations in a model MG at the time scales of microseconds (which is about 10^2^ to 10^3^ times longer than most of the conventional MD simulations in MGs) and reproducing the experimental dynamical mechanical spectroscopy in simulations. We capture the salient features of β relaxations in MGs and demonstrate that the structural excitations in MGs are string-like and that the most probable time of these string-like motions coincides with the β-relaxation time in a wide range of temperatures. Our combined findings suggest that string-like atomic rearrangements could be the structural origin of β relaxation in MGs, which might be the most convincing evidence to date.

RESULTS
=======

The details of our model system are given in Materials and Methods. Briefly, we study the relaxation dynamics of a Ni~80~P~20~ model MG by MD simulations up to a time scale of 10 μs. (see fig. S1). The samples are quenched from a high temperature with a slow cooling rate \[0.1 K/ns (10^8^ K/s)\] to avoid unnecessary structural relaxations that might obstruct the determination of β relaxations in the model MG.

[Figure 1](#F1){ref-type="fig"} shows the α-relaxation time τ~α~ as calculated from the decay of the intermediate scattering function (ISF) (fig. S2). For *T* \> *T*~g~ ≈ 540 K, τ~α~(*T*) increases markedly with decreasing *T* (filled circles) and can be fitted by a Vogel-Fulcher-Tammann equation, whereas for *T* \< *T*~g~, τ~α~(*T*) shows a mild dependence on *T* (open circle). Here, *T*~g~ is the glass transition temperature of our model MG, indicating the crossover from ergodic to non-ergodic behavior associated with the structural relaxation ([@R47]).

![Relaxation map of the Ni~80~P~20~ MG.\
The plot shows the α-relaxation time determined from ISF and MD-DMS (dynamical mechanical spectroscopy), the Johari-Goldstein secondary (β) relaxation time, and the most probable relaxation time of string-like motions as functions of temperature (top axis) and reciprocal temperature (bottom axis). VFT, Vogel-Fulcher-Tammann; CM, coupling model.](1701577-F1){#F1}

[Figure 2A](#F2){ref-type="fig"} plots the displacement probability density function *p*(*u*) (see Materials and Methods) at *T* = 440 K (in the glassy state) for different time intervals δ*t* ranging from 0.1 to 10,000 ns (10 μs). Here, *u* is the atomic displacement that is defined as *u*(δ*t*) = \|\|**r**(δ*t*) − **r**(0)\|\|, where **r** is the position vector of atoms, and *p*(*u*)*du* has the physical meaning of the fraction of atoms whose displacements fall in the range of (*u* − *du*/2, *u* + *du*/2) in a given δ*t*. As shown in [Fig. 2A](#F2){ref-type="fig"}, one of the notable features is that *p*(*u*) has multiple peaks and valleys for longer δ*t*. For example, there is only one single dominant peak for δ*t* = 0.1 ns, whereas a second peak emerges when δ*t* ≥ 1 ns and a third one for δ*t* ≥ 10 ns. For δ*t* = 1000 ns (1 μs), one can see four peaks in *p*(*u*).

![String-like atomic motion in the Ni~80~P~20~ MG.\
(**A**) Probability density function *p*(*u*, δ*t*) of atomic displacements at *T* = 440 K for different time intervals δ*t* as indicated above. (**B**) Pair distribution function *g*(*r*) for *T* = 440 K. (**C**) Two-dimensional vector field of atomic displacements at *T* = 440 K and δ*t* = 100 ns. Red (blue) vectors are for those displacement magnitudes larger (smaller) than 1.8 Å. (**D**) Typical 3D string-like configurations from simulations.](1701577-F2){#F2}

A closer examination of the data in [Fig. 2A](#F2){ref-type="fig"} reveals that the peak positions have a direct connection with the structure of the MG. [Figure 2B](#F2){ref-type="fig"} shows the pair distribution function *g*(*r*) of the MG, which relates to the probability of finding an atom at a distance *r* away from a reference atom. One can see that the *n*th peak of *p*(*u*) quantitatively matches the (*n* − 1)st peak of *g*(*r*) (where *n* ≥ 2), as indicated by the dashed lines in [Fig. 2](#F2){ref-type="fig"} (A and B).

We note that *g*(*r*) characterizes the static structure of the system, whereas *p*(*u*) is a dynamical property. The match between them implies a structural origin of the atomic rearrangements in the model MG, suggesting 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. Atomic motions like this could rationalize the match between *p*(*u*) and *g*(*r*).

[Figure 2C](#F2){ref-type="fig"} shows a typical two-dimensional (2D) slice of the vector field of the atomic displacements at *T* = 440 K and δ*t* = 100 ns. We find that atoms with larger *u* values are heterogeneously distributed, and interestingly, many of them form string-like configurations (see fig. S3 for more string-like configurations in 3D). Coordination analysis (fig. S4) indicates that atoms with larger *u* have about 2.4 nearest neighbors on average (compared to 12.0 for all the atoms; see also fig. S4), confirming their string-like configurations. [Figure 2D](#F2){ref-type="fig"} displays some typical 3D structures of isolated strings. In addition, the linear-like geometry, semi- or complete loops, and aggregated ones can be found. We note that several previous studies have identified string-like motions in other model systems ([@R48]--[@R54]). Initially, string-like motions were first observed in a Lennard-Jones liquid model ([@R48], [@R49], [@R53], [@R54]). Subsequently, Zhang *et al.* ([@R50]), Starr *et al.* ([@R51]), and Hanakata *et al.* ([@R52]) reported similar motions in other complex systems.

[Figure 3](#F3){ref-type="fig"} quantifies the fraction (*N*~str~/*N*~fast~) of atoms in string-like configurations (*N*~str~) to the total number of fast-moving atoms (*N*~fast~) as functions of temperature and time *t*. Here, the fast-moving atoms are defined as those with displacements *u*(δ*t*) ≥ 1.8 Å, which corresponds to the position of the first minimum in *p*(*u*, δ*t*) in [Fig. 2A](#F2){ref-type="fig"}. One can see that *N*~str~*/N*~fast~ exhibits a finite maximum at each temperature studied, the height of the peak grows, and the associated time scale increases with lowering temperature. This indicates the increasing degree of cooperative motions and dynamical heterogeneity. For instance, at *T* = 440 K and δ*t* ≈ 100 ns, about 80% of the fast-moving atoms belong to strings, which implies that string-like motions dominate the low-temperature structural rearrangements of MGs.

![Dynamic evolution of string-like motions.\
The plot shows the fraction of atoms in string-like configurations to the total number of fast-moving atoms as a function of time for different temperatures as indicated.](1701577-F3){#F3}

Evidently, the peak position of *N*~str~/*N*~fast~ measures the most probable time τ~str~ for the observation of these strings. We summarize the temperature dependence τ~str~(*T*) in [Fig. 1](#F1){ref-type="fig"} (squares). These time scales separate from the α relaxation when τ~α~ \> 10 ps, and the difference between them grows with decreasing temperature. A crossover of both τ~str~(*T*) and τ~α~(*T*) can be observed near *T*~g~ as a result of the departure from ergodicity. Moreover, these features can be captured by the coupling model as proposed by Ngai and Paluch ([@R12]): τ~str~ = (τ~0~)^*n*^(τ~α~)^1−*n*^, where τ~0~ = 1 ps and *n* = 0.25. We notice that this equation was originally introduced for the relation between β and α relaxations. It is remarkable that string-like dynamics also follow. This implies that the observed strings are relevant for the origin of β relaxation in MGs.

We next establish this connection by studying the β relaxation from MD simulation of dynamical mechanical spectroscopy (MD-DMS), using the same approach applied in previous experiments ([@R13]). [Figure 4](#F4){ref-type="fig"} (A to D) shows the loss modulus *E*″ against temperatures derived from MD-DMS at four different periodic deformation times *t*~ω~: 100, 30, 10, and 3 ns (see Materials and Methods). For *t*~ω~ = 100 ns in [Fig. 4A](#F4){ref-type="fig"}, the prominent feature is that, besides the dominant α relaxation, there is an additional pronounced peak at lower temperature, suggesting bimodal relaxation dynamics.

![The β relaxation as revealed by MD-DMS.\
(**A** to **D**) Loss modulus *E*″ of the Ni~80~P~20~ MG at testing periods *t*~ω~ = 100 ns (A), 30 ns (B), 10 ns (C), and 3 ns (D). (**E**) A schematic illustration of the atomic structure and rearrangement of an MG-based cooperative shear model. (**F**) The corresponding potential energy as a function of atomic displacements.](1701577-F4){#F4}

As usually invoked in the literature ([@R26], [@R55]), we interpret the lower temperature process as the signature of the β relaxation, which is partly eclipsed by the α relaxation. Quantitatively, the curve *E*″(*T*) can be fitted by a superposition of functions containing two distinct peaks, which gives the characteristic temperature of β relaxation *T*~β~ = 450 ± 10 K at *t*~ω~ = 100 ns and that of α relaxation *T*~α~ = 540 ± 4 K.

[Figure 4B](#F4){ref-type="fig"} reports the same analysis for *t*~ω~ = 30 ns, which yields *T*~β~ = 485 ± 10 K and *T*~α~ = 545 ± 4 K. However, the fitted β relaxation is less pronounced than that in [Fig. 4A](#F4){ref-type="fig"}. With further decreasing *t*~ω~, the α and β relaxations become merely separable when *t*~ω~ ≤ 10 ns, as shown in [Fig. 4C](#F4){ref-type="fig"}. Eventually, for *t*~ω~ ≤ 3 ns, these two relaxations are already largely overlapped and a separation is no longer practical ([Fig. 4D](#F4){ref-type="fig"}). We note that previous MD-DMS simulations did not reveal the signature of β relaxations due to the insufficient *t*~ω~, similar to the cases in [Fig. 4](#F4){ref-type="fig"} (C and D).

The temperature-dependent relaxation time for both β and α relaxations as extracted from MD-DMS of [Fig. 4](#F4){ref-type="fig"} is displayed in [Fig. 1](#F1){ref-type="fig"}. Significantly, the data for the β relaxation (stars in [Fig. 1](#F1){ref-type="fig"}) and string-like excitations (squares in [Fig. 1](#F1){ref-type="fig"}) collapse onto a single master curve that follows the relation of the coupling model (solid curve). These results demonstrate that the two processes are inherently related and that string-like atomic movements could be the structural origin of β relaxation in MGs. Moreover, the values of τ~α~ (*T*) determined from MD-DMS agree with those from the ISF, as shown in [Fig. 1](#F1){ref-type="fig"}, thereby validating our MD-DMS.

DISCUSSION
==========

At this point, one fundamental question is why structural rearrangements of β relaxations are string-like. In a cooperative shear model ([@R56]), the potential energy of a cooperative atomic movement is described by a periodic sinusoidal elastic energy as a function of the displacements, as schematically shown in [Fig. 4](#F4){ref-type="fig"} (E and F). There is a barrier for the two nearest low-energy configurations ([Fig. 4F](#F4){ref-type="fig"}). The peak positions, both on *g*(*r*) and *p*(*u*), are associated with low-energy configurations, whereas the valleys represent higher ones. For this reason, it is understandable that *p*(*u*) has multiple peaks (valleys) ([Fig. 2](#F2){ref-type="fig"}) and that their positions are correlated with those of *g*(*r*) ([Fig. 1C](#F1){ref-type="fig"}), consistent with the formation of strings.

As schematically shown in [Fig. 4E](#F4){ref-type="fig"}, string-like motion of "*a*-*b*-*c*" can be considered as consecutive jump-like motions. We note that the dynamics associated with [Fig. 4E](#F4){ref-type="fig"} are such that atoms reside in the low-energy configuration for longer times (caging effect) and then jump to another low-energy configuration within a very short time (fig. S5).

More generally, this scenario indicates a correlation between confinement by neighbors and the time scales of larger-scale displacements of particles. Specifically, at high temperatures (for example, *T* \> 540 K in this Ni~80~P~20~ MG where the α relaxation is active), the thermal fluctuations weaken the confinement by neighboring atoms, and the string-like excitations become rare or absent, causing the disappearance of β relaxations (that is, the merging of β- and α-relaxation modes). On the other hand, at low temperatures (where the α relaxation is frozen), the neighboring confinements are strong and persistent, leading to low mobility of atoms. Thus, in this glassy state, one needs to observe the system for a sufficiently long time to identify the formation of string-like excitations and to probe β relaxations.

CONCLUSIONS
===========

In summary, by conducting MD simulations at the microsecond time scales, we have captured the structural signatures of string-like motions in MGs ([Fig. 2](#F2){ref-type="fig"}) and reproduced the salient features of β relaxation at low temperatures ([Fig. 4](#F4){ref-type="fig"}). Moreover, the string-like motions identified in MD simulations and the β relaxation observed in the mechanical response follow nearly the same dependence of relaxation time on temperature ([Fig. 1](#F1){ref-type="fig"}). Our findings establish that string-like motions could be the origin of β relaxation in MGs. This work also highlights the importance of studying glassy dynamics over longer time scales. Referring to the famous saying of Anderson ("more is different" for condensed matter physics) ([@R57]), we might expect that longer is different for glassy physics. We note that, very recently, Ninarello and co-workers ([@R58]) introduced the swap Monte Carlo simulation, which is quite efficient for equilibrating supercooled liquids and glasses. It could be useful in this regard.

MATERIALS AND METHODS
=====================

Simulation details
------------------

An open-source MD code, LAMMPS ([@R59]), was used for all the simulations. Our model system contains *N* = 32,000 atoms with the composition of Ni~80~P~20~ that interact with an embedded atom--method potential, as developed by Sheng and co-workers ([@R60], [@R61]). The glassy samples were prepared by quenching a liquid from 1000 to 200 K at a rate of 0.1 K/ns (10^8^ K/s). Such a slow cooling is about two to three orders of magnitude slower than most of the MD simulations of MGs conducted thus far. We used NPT ensembles (that is, constant number of atoms, pressure, and temperature) during the quenching. The external pressure was adjusted to around zero. Periodic boundary conditions were applied for all the simulations.

Dynamical mechanical spectroscopy
---------------------------------

The dynamics of β relaxations were studied by an approach of MD-DMS simulations ([@R62], [@R63]) that numerically reproduces the protocol of real DMS experiments ([@R11], [@R13]). Specifically, at a temperature *T*, we applied a sinusoidal strain ε(*t*) = ε~A~ sin(2π*t*/*t*~ω~), with a period *t*~ω~ (related to frequency *f* = 1/*t*~ω~) and a strain amplitude ε~A~, along the *x* direction of the model MG, and the resulting stress σ(*t*) and the phase difference between stress and strain δ were measured and fitted with σ(*t*) = σ~0~ + σ~A~ sin(2π*t*/*t*~ω~ + δ). The strains were applied by an appropriate change in length of the box size, whereas stresses were measured from each component of the pressure of the system, as implemented in LAMMPs. From these values, storage (*E*′) and loss (*E*″) moduli are calculated according to *E*′ = σ~A~/ε~A~ sin(δ) and *E*″ = σ~A~/ε~A~ cos(δ), respectively. We used a strain amplitude ε~A~ = 0.71% for all the MD-DMS, which ensures that the deformation is in the linear response regime (that is, the deformations do not change the structure of the MG). We used NVT ensembles (that is, constant number of atoms, volume, and temperature) for the MD-DMS. We used 10 cycles for each MD-DMS.

Displacement probability density function *p*(*u*)
--------------------------------------------------

The probability density function *p*(*u*) is defined as *p*(*u*) = \[*P*(*u* + Δ*u*) − *P*(*u*)\]/Δ*u*, where *P*(*u*) is the cumulative distribution that quantifies the probability of finding *X* ≤ *u*. We used Δ*u* = 0.01 Å for all the calculations. By definition, the probability density is normalized according to $\int_{0}^{\infty}\mathit{p}(\mathit{u})\mathit{d}\mathit{u} = \mathit{P}(\infty) = 1$.

String analysis
---------------

To identify string-like motions in glasses, we analyzed all the atomic trajectories at two times, *t*~0~ and *t*~0~ + δ*t*. One pair of atoms (for example, *i* and *j*) is considered within a string over the time scale \[*t*~0~, *t*~0~ + δ*t*\] if the position of atom *i* at time *t*~0~ + δ*t* is same as the position of atom *j* at *t*~0~ within a numerical accuracy (δ*r* = 0.1 Å in this work), which can be expressed as $\left\| {\overset{\rightarrow}{\mathit{r}}}_{\mathit{i}}(\mathit{t}_{0} + \delta\mathit{t}) - {\overset{\rightarrow}{\mathit{r}}}_{\mathit{j}}(\mathit{t}_{0}) \right\| \leq \delta\mathit{r}$. Different atom pairs are considered to be in a same string if they share at least one same atom.

Supplementary Material
======================

###### http://advances.sciencemag.org/cgi/content/full/3/11/e1701577/DC1

The computational work was carried out using the Tianhe-1A of the National Supercomputer Center in Tianjin, China, the Tianhe-2 of the National Supercomputer Center in Guangzhou, China, and the Gesellschaft für wissenschaftliche Datenverarbeitung mbH Göttingen (GWDG), Germany. **Funding:** H.-B.Y. acknowledges the support from the National Science Foundation of China (NSFC 51601064) and the National Thousand Young Talents Program of China. K.S. acknowledges the support from the German Science Foundation within the FOR 1394, P1. **Author contributions:** H.-B.Y. conducted simulations and analyzed the data. R.R. and K.S. proposed the research and contributed to the mechanisms. All the authors co-wrote the manuscript. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Supplementary material for this article is available at <http://advances.sciencemag.org/cgi/content/full/3/11/e1701577/DC1>

fig. S1. Properties change during quenching of the Ni~80~P~20~ MG.

fig. S2. Intermediate scattering function.

fig. S3. Typical string-like configurations.

fig. S4. Coordination analysis.

fig. S5. Typical single atomic trajectories at *T* = 440 K and the duration of the observation time *t* = 1000 ns for three single atoms.
