Active microrheology of a bulk metallic glass

See allHide authors and affiliations

Science Advances  17 Jul 2020:
Vol. 6, no. 29, eaba8766
DOI: 10.1126/sciadv.aba8766


The glass transition remains unclarified in condensed matter physics. Investigating the mechanical properties of glass is challenging because any global deformation that might result in shear rejuvenation would require a prohibitively long relaxation time. Moreover, glass is well known to be heterogeneous, and a global perturbation would prevent exploration of local mechanical/transport properties. However, investigation based on a local probe, i.e., microrheology, may overcome these problems. Here, we establish active microrheology of a bulk metallic glass, via a probe particle driven into host medium glass. This technique is amenable to experimental investigations via nanoindentation tests. We provide distinct evidence of a strong relationship between the microscopic dynamics of the probe particle and the macroscopic properties of the host medium glass. These findings establish active microrheology as a promising technique for investigating the local properties of bulk metallic glass.


When a glass-forming liquid is rapidly cooled below its glass transition temperature, it falls out of equilibrium, and its viscosity, η, increases by more than 15 orders of magnitude (1, 2). Notably, the structural relaxation time of a glass, τα ∝ η, exceeds any other time scale in the realm of condensed matter physics and is comparable to cosmological time scales (3). The origin of this slowing down remains the primary challenge in the physics of glass transition, i.e., whether the long time scale originates from collective behavior (a thermodynamic description) (4), similar to critical phenomena, or the slowing down occurs at the microscopic level (a kinetic origin). A widely accepted picture of the slow progress of a supercooled liquid toward equilibrium is based on a trap-escape behavior: A particle is transiently trapped in a cage created by its neighbors; over a short time, τβ, the particle rattles inside the cage and escapes only after a long time, τα. This time scale diverges at the glass transition temperature, at which the cages become permanent. This interpretation is well supported by the characteristic two-step relaxation of a supercooled liquid, as observed in neutron scattering experiments and simulations (5). Several competing theories have aimed at explaining the origin of the slowing down of a glass. Among them, the mode-coupling theory (MCT) has successfully described power-law β-relaxation, stretched exponentials of α-relaxation, and the time-temperature superposition principle (6). An intriguing discovery in experiments (7, 8) and simulations (9, 10) has motivated a thermodynamic description of the glass transition: Localized particles, unable to diffuse away, form mesoscopic clusters whose typical size increases upon cooling (11). These mesoscopic structures, known as dynamic heterogeneities (DHs), whose origin is yet another controversy (static versus dynamic origin), suggest that a supercooled liquid can be treated as a complex fluid: An enlightening work by Furukawa and Tanaka (12) established the complex fluid view for a supercooled liquid. This interpretation contrasts with the original MCT, which was an extension of the simple liquid theory, although recent developments of the MCT have been reconciled with DHs (13).

Despite intense research in recent decades on the underlying origins of the glass transition, a consensus has not yet been established. Furthermore, understanding of the mechanical properties of glasses, which are investigated by shear deformation, is also elusive. The reason for this lack of consensus is fourfold: (i) Any global perturbation requires an extremely long time scale for relaxation; (ii) shear forces may rejuvenate the glass, thus relocating the glass to shallower energy minima (14, 15) and rendering the glass more prone to fluctuations; (iii) a realistic deformation rate for an atomistic glass is a technical impossibility in computer simulations (16, 17); and (iv), because of DHs, a glass is heterogeneous, and a global deformation allows for exploration of only an average mechanical response. Here, we overcome these problems by investigating the mechanical properties of a glass at the microscopic level. With this technique, known as active microrheology, a singled-out probe particle is pulled through a dense host liquid (18). Investigation of the probe particle dynamics reveals nonequilibrium properties of the host medium, such as viscosity.



We perform active microrheology of a probe particle with a constant velocity U. Density map snapshots of the system are given in fig. S1. The size of the probe particle mimics that of a typical length scale of DHs. We obtain the response of the host medium by measuring the force acting on the particle, i.e., the friction, F=F·n̂, where n̂=U/U. The host medium is modeled by a Kob-Andersen mixture in three dimensions (3D), whose thermodynamic glass transition temperature, the Kauzmann temperature, is TK = 0.3 (more information regarding the model is given in Methods). Before performing microrheology simulations, we use a state-of-the-art method to deeply anneal the host medium with cyclic shearing (19). This step minimizes the occurrence of plastic events, thus rendering the system less prone to fluctuations. In Fig. 1A, the magnitude of the frictional force acting on the probe particle, F, is depicted versus the probe velocity U. Each color corresponds to a different temperature T. A complicated nonlinear velocity-force relationship is found. The velocity-force curves can be separated into two distinct classes. In the first class, at high T, the friction F vanishes as U → 0. However, in the second regime, at low T, the force F is finite at U → 0. To our knowledge, this is the first report of such behavior for a bulk metallic glass, although a similar threshold has been found for soft colloids (20, 21). This threshold force is akin to the yield stress of a paste (22). The emergence of a nonzero threshold force F is reminiscent of a phase transition and can be interpreted as an arrest transition: The probe particle will not delocalize if the probe is driven by a force smaller than the threshold force. This crossover from delocalization to an arrest regime occurs at low T, where the glass transition occurs. Therefore, a close relationship may exist between the microscopic dynamics of the probe particle and the macroscopic properties of the host medium. We will demonstrate the existence of such relationship in the remainder of this paper. The crossover region from the arrest regime to delocalization, which is marked by the shaded region in Fig. 1A, is wide. We will demonstrate that this wide region is an advantage of microrheology.

Fig. 1 Velocity-force curves.

(A) Friction, F, versus the probe velocity, U, for various temperatures across the glass transition temperature. Each point is an ensemble average over 500 independent realizations. (B) Successive slope of the velocity-force curves m=lnF/lnU versus the probe velocity. A negative-positive crossover of the curvature occurs in the range of T = 0.375 to 0.300, marking the critical transition region. Equation 2 is fitted to each curve to compute the corresponding asymptotic exponent q in the critical region (solid lines). (C) The asymptotic exponent q varies monotonically with temperature T across the transition region. The data are interpolated by splines, thus enabling us to obtain the asymptotic exponent for any temperature within the transition region. The shaded region shows the range over which the data can be collapsed. (D) An excellent collapse is obtained for Tmc = 0.36875, with an asymptotic exponent of q = Γ/Δ = 0.17303, in which Γ = 0.475 and Δ = 2.75.

Scaling ansatz

To fully capture the nonlinear behavior of the velocity-force curves, we seek a scaling ansatz reminiscent of that for magnetization in an Ising model. In this approach, the friction, F, depends only on a scaling variable x = U/∣TTmcΔ viaF(U,T)=TTmcΓF(UTTmcΔ)(1)in which δT = ∣TTmc∣ is the distance from a critical delocalization temperature Tmc, F() is a scaling function, Γ is a parameter-dependent exponent, and Δ is a universal gap exponent. At T = Tmc, scale invariance requires FUq, in which q is a critical exponent, and Eq. 1 confers q = Γ/Δ (a formal derivation is given in note S1). We previously developed a systematic framework to extract the critical density and critical exponent for the macrorheology of jamming (23). Here, we generalize this framework for microrheology to extract Tmc and the critical exponent q. This technique is based on the successive slope of the velocity-force curves, defined as m = ∂ ln F/∂ ln U: At T = Tmc, the curvature of the m − U curves crosses over from negative to positive. This aspect can be used as a criterion to find the critical region. As shown in Fig. 1B, we computed the successive slope of the velocity-force curves, m, for various temperatures. The slope m can be analytically calculated from Eq. 1, and at T = Tmc, the slope becomes m = q. From Fig. 1B, the curvature of mγ̇ curves crosses over from negative to positive in the range of T = 0.300 to 0.375. Clearly, at T = 0.375 (highest bold symbols), the data curve upward for the asymptotic limit of U → 0; similarly, at T = 0.300 (lowest bold symbols), the data curve downward for the same limit. Therefore, the critical temperature must be in the range 0.300 < T < 0.375. We also observe a strong correction to scaling for U > 10−2, indicating that at T = Tmc, a correction to scaling of the form F = Uq1 + κ2Uω/z) is required, in which ω is the leading correction-to-scaling exponent, and z is the dynamic exponent. Therefore, at T = Tmc, the successive slope becomesm=q+κUω/z(2)

Because the critical temperature Tmc is unknown, we fit Eq. 2 to the data in the critical region, 0.300 < T < 0.375. The resulting lines are depicted by solid curves in Fig. 1B, and we display the values of the fit parameters in Table 1.

Table 1 Fitting parameters for Eq. 2 displayed by the solid lines in Fig. 1B.

View this table:

The magnitude of the prefactor in Eq. 2 is crucial in this analysis because it determines the reliability of the fit. For high-quality data, κO(1) must be obtained. From Table 1, κ is on the order of unity and does not change markedly with temperature. This behavior demonstrates the quality of our data, owing to the careful preparation of a well-annealed medium, as well as the rarity of shear rejuvenation in microrheology. The asymptotic exponent q changes rapidly with T, as illustrated in Fig. 1C. We interpolate the q − T data with splines, thus enabling us to obtain the corresponding q for any T in the critical region. If F scales according to Eq. 1, then rescaling the data with FTΓ versus UTΔ must lead to a collapse of the data into a master curve. This rescaling of the flow curves is of interest for both fundamental and technological applications for predicting the flow properties of materials (17, 24). The ratio of these exponents q = Γ/Δ is obtained through interpolation, and we vary the exponent Γ at a given T to achieve a collapse. We obtain a collapse of the data within only a narrow range of Tmc = 0.36875 ± 0.0062 (shaded area in Fig. 1C), with the optimal collapse at the center of the interval. The error bars account for a mere 1.5% uncertainty. The narrowness of the interval over which a collapse can be obtained is highly similar to the results for conventional critical phenomena. In contrast, the flow curves of athermal systems undergoing a jamming transition exhibit an uncertainty of approximately 5% (23), thus leading to the long-standing controversy regarding the exponents involved in jamming (25). We conclude that the critical delocalization temperature is Tmc = 0.36875 ± 0.0062, which is approximately 23% higher than the thermodynamic glass transition temperature, at TK = 0.3. The delocalization temperature Tmc is greater than the thermodynamic glass transition temperature TK, a result consistent with the findings of Habdas et al. (26). The critical exponent associated with this temperature is q = 0.17303 ± 0.0248, and the crossover exponents are Γ = 0.475 ± 0.068 and Δ = 2.75. With these exponents, an excellent collapse is obtained, as shown in Fig. 1D, thus confirming a phase transition for the probe particle. Although we obtain an equilibrium regime at high T, i.e., FU1, the velocity-force dependence is otherwise highly nonlinear. Therefore, this collapse is primarily over a nonequilibrium region. This behavior is remarkable because critical scaling collapses generally occur at equilibrium. We also observe that the velocities over which a linear response is seen decrease as Tmc is approached, in agreement with a report by Williams and Evans (27).

Equation 1 suggests that the inherent friction, Finh = F(U → 0, T), in the system scales as Finh ∝∣TTmcΓ. The numerical value of the exponent Γ within the measurement error bars is equal to the exponent β = 0.41 − 0.42 of the percolation theory in 3D. The latter describes the scaling of the order parameter P, the probability that a site belongs to the infinite cluster; it is zero below the critical point, and it scales as a function of the distance from the critical point P ∝ (ppc)β. This may indicate that the inherent friction is generated by percolation of immobile particles (DHs) near the glass transition point. To examine this possibility, a systematic percolation analysis would be required. Such an investigation is beyond the scope of this work, but we will address it in the near future.


We now focus on the properties of the time series of the force acting on the probe particle F = {F(t)}, and we investigate whether the structural relaxation process can be traced via a friction time series. In fig. S2 (A and B), we display the friction time series at T = 0.150 and 0.466, respectively, for a probe velocity of U = 0.0379269. As shown in fig. S2B, in the supercooled state, the friction changes rapidly with time; in contrast, in fig. S2A, corresponding to the glass phase, the friction exhibits a strong correlation. The time interval over which the friction is correlated is marked by an arrow as a guide to the eye (the shaded area). This apparent correlation may be related to the formation of cages. To investigate this possibility, in Fig. 2 (A and C), we plot the autocorrelation function for friction, CF(τ) = 〈F(τ)F(0)〉, for T = 0.200 and 0.466, respectively. 〈⋯〉 presents an ensemble average over 500 independent realizations (in fig. S3, we show a similar plot for T = 0.350). Each curve begins with an initial relaxation, which corresponds to the ballistic motion of particles at very small time scales. All the different curves with different U values are superimposed for a short time scale, thus indicating that the short-term dynamics is unaffected by the probe motion. This behavior is analogous to that of sheared liquids at very small time scales, which is similarly unaffected by shear (28). The curves level off with a dip followed by a damped oscillation and eventually a plateau. The dip and the damped oscillation correspond to the oscillatory nature of the initial system preparation, owing to shear annealing (29). The damped oscillation of the correlation function is known as quenched echoes (30). Kob and Barrat (31) have reported such oscillations for Lennard-Jones (LJ) glasses annealed with oscillatory protocols. However, Shiba et al. (32) have shown similar oscillations originating from the propagation of acoustic phonons. Kob and Barrat ruled out the possibility of this dip being related to the boson peak. After the damped oscillations, a plateau with a length of several orders of magnitude appears. This trend is akin to the β-relaxation process (5). After this rattling period, a long relaxation process terminates the curves. This final long relaxation is akin to an α-relaxation process, corresponding to the escape of a particle from a cage. The behavior of the system both below and above the glass transition temperature is consistent with that in shear-driven structural glasses (28, 33). Therefore, the relaxation of friction is described by a two-step relaxation process, thus revealing a strong relationship between dynamics of the probe particle and the two-step structural relaxation of the host medium glass. In the remainder of the paper, we provide more evidence of this intimate relationship.

Fig. 2 Two-step relaxation of friction.

The autocorrelation function CF(τ) = 〈F(0)F(τ)〉 of the time series of the force acting on the probe particle F versus the time window τ for T = 0.200 and 0.466 is depicted in (A) and (C), respectively. Each curve represents an average of 500 realizations. Each curve has an initial relaxation, which corresponds to ballistic transport at very small time scales. This fast relaxation is followed by a plateau that is akin to a β-process (rattling in cages). Each plateau starts with a damped oscillation, known as a quenched echo (30, 31), which corresponds to the oscillatory nature of the initial annealing of the system or phonon propagation (32). The plateau terminates with a long relaxation process akin to an α-process (escape). We plot CF(x) versus x = τ × Uν, in which ν = −(q − 1) = 0.95 ± 0.1 is the force thinning exponent in (B) and (D). The data collapse onto a compressed exponential function of the form Eq. 3 [dashed lines in (B) and (D)]. The collapse is akin to the time-shear superposition principle in spin glass and MCT (1, 33). The compressed (faster) exponential is a manifestation of the super-diffusive motion of the probe particle, corresponding to super-diffusive relaxation. In (E), the compressed exponent b is depicted versus the probe velocity for various values of T.

Time-force superposition

In Fig. 2 (B and D), CF(τ) is plotted versus x = τ × Uν, in which U is the velocity of the probe particle and ν = −(q − 1) = 0.95 ± 0.1 is the force thinning exponent (q is the critical exponent of the scaling ansatz in Eq. 1); here, an excellent collapse of the data into a master curve is obtained. The collapse of the correlation function is a time-force superposition, which is reminiscent of the “time-shear superposition principle” predicted by both the MCT (1) and spin glass theory (33). The force thinning exponent ν is consistent, within the error bars, with the shear thinning exponent reported by Furukawa et al. (12), ν = 0.8, and is very close to the values reported by Berthier and Barrat (16), ν = 2/3, and Varnik (28), ν = 1. Accordingly, our force thinning exponent, ν, and the critical exponent, q = 1 − ν, within the error bars are consistent with findings reported in the literature for macrorheology.

Compressed exponential

The master function (dashed lines in Fig. 2, B and D) is given by a compressed exponential function of the following form (dashed lines in Fig. 2, B and D)CF(x)e(x/σ)b(T,U)(3)in which b(T, U) ≥ 1 is a temperature- and velocity-dependent (compressed) exponent. In Fig. 2E, the compressed exponent is plotted versus the probe velocity, U, for various temperatures. The exponent is generally larger than unity and reaches a ballistic value of 2 in the supercooled regime for large U values. This result is in contrast with the stretched-exponential form b < 1 predicted by the MCT. The faster-than-exponential relaxation b > 1 stems from the super-diffusive motion of the probe particle δrprobe2t2. This result sheds light on an open debate on the origin of the compressed exponential relaxation in colloidal systems and glass formers (34, 35). Another important characteristic of a glass is aging (31): The longer an observer waits, the slower the glass becomes. This characteristic is simply due to particles trapped in the cages. However, shear rejuvenates a glass and stops aging (28, 36): Cages are continually crushed by shear rejuvenation. Although glass is at nonequilibrium, a sheared glass reaches a stationary state similar to equilibrium, and the correlation functions depend only on the time difference. We examined aging in active microrheology with a two-point correlation function, CF(τ, τw), in which τw is the waiting time. In agreement with findings for sheared glasses, we found that CF(τ, τw) is time invariant and does not exhibit any aging effects. However, we do not suggest that this behavior is related to shear rejuvenation because our probe is local. This behavior is expected to arise because the activated probe particle continually breaks the cages and is not trapped in any cage.

Power spectral density

In Fig. 3, we display the power spectral density (PSD) of the friction time series F = {F(t)}, defined as PSD=F{F(t)}2, with F being the Fourier transform to the frequency domain v. At low frequencies, which correspond to long-term relaxations, we obtain a power law distribution ν−θ, where θ is an exponent. The exponent θ has a temperature dependence that can be rationalized by the absolute value of the Fourier transform of the compressed exponential function of the α-relaxation time, CF(τ) ≃ e−(τ/τα)b(T). At low temperatures, the exponent θ converges to unity, which corresponds to a 1/f colored-noise process. However, the exponent θ decreases with increasing temperature. At very high temperatures, the friction is expected to become a random process described by white noise, i.e., θ = 0, thus justifying the decreasing trend in θ as a function of T. At high frequencies, corresponding to the ballistic motion of particles at very short time scales, all curves with different T superimpose. This finding is consistent with the short-term behavior of CF(τ). This time series analysis approach can be generalized to investigate many other interesting properties of glass transition with the time series analysis obtained by active microrheology. For instance, a comparison of the PSD of the friction time series with that of the velocity time series can be used to investigate the violation of the Stokes-Einstein relation. These methods may provide an interesting future avenue for research aiming to connect microrheology with the time series analysis and glass transition.

Fig. 3 Power spectral density (PSD) of friction time series.

Absolute value of the Fourier transform of the autocorrelation function, PSD(ν) = ∣∫ CF (τ)e−2πiντdτ∣, the PSD, versus frequency ν for various T. The probe velocity is U = 10−3. At low frequency, a temperature-dependent scaling behavior is observed, owing to the temperature dependence of the compressed exponent b = b(T) of the autocorrelation function CF (τ) = exp (−(τ/τα)b). At very low T, the dependence converges to the well-known 1/ν colored-noise behavior because of the correlation in the data. As might be expected, the exponent of PSD decreases with increasing temperature. At very high T, the friction must be completely random, and therefore, the corresponding PSD must be a horizontal line of white noise, i.e., θ = 0.


A conventional method of exploring the mechanical properties of supercooled liquids is macrorheology. We now compare the results obtained from active microrheology and macrorheology. We performed simulations of a simple shear applied to the host medium by using the Lees-Edwards boundary conditions. The resulting flow curves are shown in Fig. 4A. For most of the flow curves at low temperature, the shear stress saturates to a shear rate–independent regime at low γ̇. This saturation signals the emergence of yield stress and a transition from fluid to solid states. Moreover, the crossover from fluid to solid regimes occurs in an extremely narrow range between T = 0.400 (squares) and 0.375 (diamonds). This region is marked by a shaded area whose width can be compared to a similar region in Fig. 1A. The comparison notably reveals that macrorheology provides a much narrower critical range, if any, than microrheology. We will see that the narrowness of the transition region makes a systematic scaling analysis nearly impossible for macrorheology. In Fig. 4B, we plot the successive slope of the flow curves. In contrast to that in microrheology, the successive slope is very noisy, with large fluctuations for U → 0; these fluctuations are characteristic of shear rejuvenation and slow dynamics at low T. With macrorheology, similarly to microrheology, the crossover region is defined as the range at which the curvature of mγ̇ data cross over from positive to negative (black), which occurs in the range between T = 0.400 and 0.375. A systematic analysis similar to that in microrheology is not possible because of the large fluctuations in the data. Therefore, we cannot obtain the transition temperature or the critical exponent with a systematic scaling analysis. As many previous authors have noted, we cannot observe a critical transition. In Fig. 4C, we collapse the flow curves by using the exponents obtained from a systematic microrheology scaling analysis in Fig. 1. The collapse clearly is not as strong as that with microrheology. Moreover, to achieve a data collapse, we adopt a higher critical temperature Tc = 0.41, thus indicating that shear rejuvenation prevents the system from reaching equilibrium; as a result, the system falls out of equilibrium at higher temperatures.

Fig. 4 Flow curves.

(A) Shear stress, σxy, versus the shear rate, γ̇, for different temperatures. In the shaded region, a crossover from fluid to solid can be seen in a very narrow range. As a result, the critical region is very small. (B) Successive slope of the flow curves, m, versus the shear rate γ̇. In contrast to the successive slope observed for active microrheology, large fluctuations are observed for U → 0. These large fluctuations prevent a systematic scaling analysis to extract the critical exponents. (C) Flow curves collapsed with the exponents obtained from active microrheology at a critical temperature Tc = 0.41. The collapse has lower quality than that for microrheology.

Bridge between micro- and macrorheologies

Creating a link between macro- and microrheology is crucial. To this end, we seek a direct comparison of microviscosity (ηm) and viscosity (η), defined asηm=(16πa)FU,η=σxyγ·(4)where a is the radius of the probe particle. A scaling ansatz for microviscosity can be obtained from Eq. 1ηm(U,T)=TTmcΓΔFη(UTTmcΔ)(5)

A similar scaling ansatz can also be written for viscosity if ηm, Tmc, and U are replaced by η, Tc, and γ̇, respectively (a formal derivation is given in note S2). To compare the numerical values of ηm and η, we define an effective shear rate γ̇e=U/D, where D = 2a is the diameter of the probe particle. This effective shear rate enables us to compare microrheology with macrorheology. In Fig. 5, the scaling collapse suggested by Eq. 5 for microviscosity (ηm/δTΓΔ versus γ̇e/δTΔ) superimposes on the scaling collapse for viscosity (η/δTΓΔ versus γ̇/δTΔ). This finding is remarkable because the microviscosity and viscosity superimpose without any additional rescaling along the x or y axis. Furthermore, this finding is notable because (i) it reveals a robust consistency of active microrheology and macrorheology after intensive averaging over time and different realizations, and (ii) it provides a simple means for direct comparison of velocity-force curves of microrheology with flow curves of macrorheology.

Fig. 5 Superimposed η and ηm.

To bridge the gap between micro- and macrorheology, we examine the scaling collapse of microviscosity (ηm bold symbols) and viscosity (η) [defined in Eq. 4]. We define an effective shear rate as γ̇e=U/D, where D = 2a is the diameter of the probe particle, for microrheology. The scaling collapse of viscosity superimposes with that by microviscosity, without any additional rescaling along the x or y axis.

Distribution of microviscosity

We now demonstrate the superiority of microrheology to macrorheology, the latter of which yields only an average response. So far, we have demonstrated that the macrorheology can be attained by microrheology when the local properties, such as microviscosity, are spatially and temporally averaged. However, microrheology is superior to macrorheology in that it can probe the local properties of the system, such as local microviscosity and the local elastic modulus. To demonstrate this capability, in Fig. 6, using time series of microviscosity ηm, we plot the probability distribution function (PDF) of ηm. The horizontal axis is (ηmηm¯)/σηm, where the microviscosity is subtracted from the mean, ηm¯, and σηm is the width of each distribution. The vertical dashed line shows the limit of macrorheology given by a delta function, δ(ηηm¯). We expect that the heterogeneity of the medium at low temperature must manifest in PDFs of microviscosity. To quantitatively investigate this possibility, in the inset of Fig. 6, we depict the skewness of PDFs, defined as (ηmηm¯)3/σηm3. The skewness increases as the temperature decreases, a result that we attribute to the growth of DHs at low T. Therefore, microrheology can be used to investigate the heterogeneous distributions of the mechanical properties; a similar analysis can be easily performed for the elastic modulus.

Fig. 6 Probability distribution function of microviscosity.

Main: Probability distribution function of local microviscosity ηm, defined by Eq. 4, versus (ηmηm¯)/σηm, where ηm¯ is the mean and σηm is the width of the distribution. Different colors correspond to different temperatures given in the legend. All distributions are bell shaped. Inset: The skewness of PDFs increases as the temperature decreases. These results provide direct evidence of the heterogeneous distribution of mechanical properties due to DHs.


Here, we investigated the dynamics of a probe particle driven into a host glass-forming liquid. The friction acting on the probe exhibits a second-order phase transition. The microscopic origin of friction is one of the fundamental questions in tribology (37). Our results establish a nontrivial correspondence between the emergence of friction as a result of the approach to glass transition. We discovered a two-step relaxation process, a characteristic of the structural relaxation of a supercooled liquid, in the time autocorrelation function of the friction. Moreover, the autocorrelation function exhibits a time-force superposition reminiscent of the time-shear superposition principle in MCT (1) and spin glass theory (33). This finding strengthens the relationship between the dynamics of the probe particle and the structural relaxations of the host medium, on the basis of existing theoretical lines of argument. We further demonstrated that the autocorrelation function of friction is a compressed exponential, thus shedding light on a recent controversy regarding the origin of compressed exponential relaxation in active matter (34, 35). We compared the results obtained by active microrheology and macrorheology, revealing that spatiotemporally averaged microviscosity becomes equivalent to viscosity, thus providing clear evidence of the consistency of micro-macro rheologies. However, microrheology becomes superior to macrorheology because the former can provide the spatial distribution of microviscosity. We found that the skewness of PDF of microviscosity becomes larger as the temperature is decreased, a direct sign to the growth of heterogeneities at low temperatures.

The delocalization transition of the probe particle investigated by active microrheology has many common features with the so-called depinning transition. The depinning transition can be used to investigate collective behavior in many nonequilibrium systems. Specifically, it examines interacting particles driven over heterogeneous substrates, and it has abundant applications in vortex lattices in type II superconductors and driven colloids over optical periodic substrates. This topic has parallels with nonequilibrium phase transitions and directed percolation. More details can be found in an excellent recent review by Reichhardt and Reichhardt (38).

The following applications can be expected in microrheology studies of glasses or supercooled liquids. The spatial heterogeneities of the elastic modulus and viscosity can be observed in the systems. A previous study (39) has shown that elastic moduli are expected to be spatially heterogeneous in glass because the response on the strain fields for weak macroscopic flow becomes heterogeneous. In another study on supercooled liquids (12), a characteristic length scale has been observed from the wave number dependence of Fourier transformed viscosity. However, no studies have directly obtained the spatial distribution of the elastic modulus and viscosity. In the present study, we statistically confirmed the heterogeneous distribution of these mechanical properties. In the future, we expect to obtain the spatial distribution of these quantities. The heterogeneous distribution of these static mechanical properties suggests a relationship with the DH characterized by the particle-displacement field. Our analysis should be essential to understanding the origin of the DHs.

The low-temperature properties of glasses have recently been the subject of intensive investigations. Originally known as the Gardner transition in spin glass theory, the free energy landscape of a glass has been proposed to become fractal-like at low temperatures: A local minimum divides into many metabasins in the free energy landscape. As a result of this proliferation of the metabasins, local cages divide into several smaller cages upon cooling (or compression) (40). Therefore, the Gardner transition is associated with a change in the local morphology. As a result, a global deformation would not only hinder the exploration of such delicate local morphological changes but also might result in a loss of local information, thus potentially explaining why the existence of the Gardner transition at finite spatial dimensions has been debated. We believe that active microrheology provides a strong experimental and computational tool for investigating the low-temperature properties of glasses and, eventually, shedding light on the controversy regarding the existence of the Gardner transition at finite dimensions. Our work provides a basis for such explorations.



We performed large-scale 3D simulations of a 80 : 20 binary mixture (Kob-Andersen model) interacting via a smoothed LJ potentialUij(r)={4ϵij[(σijr)12(σijr)6]+4ϵij[σ0ij+c2ij(rσij)2],r<rcij0,r>rcij(6)where i, j = A, B. The parameters are represented in dimensionless LJ units in which the length, energy, mass, and time units are σ, ϵ, m, and τ=σ2m/ϵ, respectively. We fix the particle number density to ρ = (NA + NB)/V = 1.2/σ3. The time step of integration is Δt = 0.01τ, and ϵij and σij are set to ϵAA = 1.0ϵ, ϵAB = 1.5ϵ, ϵBB = 0.5ϵ, σAA = 1.0σ, σAB = 0.88σ, and σBB = 0.8σ. Regardless of the particle type, the mass is m = 1. The coefficients c0ij and c2ij are determined via Uij(rcij) = 0 and dUij/drr=rcij = 0. The cutoff, rcij, is set to rcij = 2.5σij, and the total number of particles is N = 4000. The thermodynamic glass transition temperature for this system, the Kauzmann temperature, is TK = 0.3. All simulations were performed with large-scale atomic/molecular massively parallel simulator (LAMMPS). Units are not shown explicitly throughout the text for readability unless otherwise specified.

Preparation of initial configurations

We prepared the initial configurations in three steps. (i) Particles are randomly inserted in the simulation box. (ii) A cosine soft pair potential of the form Uij(r) = A[1 + cos (πr/rc)] is applied to all pairs of interacting particles within the cutoff, rc, regardless of particle type. We then integrate the equations of motion of the particles by using the microcanonical ensemble (NVE) integration scheme followed by canonical ensemble (NVT) at T = 0.466ϵ/kB. This step is crucial to minimize the overlap between particles. After minimization, the soft cosine potential is replaced with the smoothed LJ potential given by Eq. 6. (iii) The system is annealed with cyclic-shearing deformation (15, 29) to obtain a stable glass less prone to fluctuations. The SLLOD algorithm was used to perform nonequilibrium molecular dynamics (MD) simulations with the Lees-Edwards boundary conditions.


We use a singled-out probe particle to perform microrheology. The probe particle is activated by a constant velocity U to mimic constant-velocity active microrheology (18). The probe particle interacts with the following expanded and shifted LJ pair potentialU(r)={,rΔ4ϵ[(σrΔ)12(σrΔ)6+14],Δ<rΔ+rc0,r>Δ+rc(7)

Unless otherwise stated, Δ = 2.0σ and rc = 21/6σ. Accordingly, the probe particle is fivefold larger than the A-type particle. This is a typical length scale for DHs in supercooled liquids (12). At t = 0, the probe particle is placed at the boundary, and after activation, the probe particle moves into the host at the prescribed velocity. We record the force acting on the probe particle in the direction of the velocity, F = −Fn where n = U/∣U∣. Here, F = ∑jfj, where fj is the force exerted by a neighboring particle, and the sum spans all neighbors of the probe particle. The force is averaged over time for 500 independent realizations. The total CPU time for the simulations is approximately 2 months for a cluster of 1656 cores. We gathered a dataset of 20 TB for all the simulations.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We especially thank T. Voigtmann for many fruitful conversations during the course of this work and acknowledge communications with I. Procaccia and L. Berthier. We thank the Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Linux cluster system) for this work and are grateful for consultations with H. Kim. Funding: KIAS authors acknowledge KIAS individual grant nos. PG054602 (S.H.E.R.) and PG013604 (H.P.). T.K. acknowledges Japan Society for the Promotion of Science (JSPS) KAKENHI (No. 19K03767 and 20H05157). This work was supported by the Samsung Research Funding Center for Samsung Electronics 566 under project number SRFC-MA1602-02. Author contributions: S.H.E.R. and J.W.Y. conceived the project. J.W.Y. performed the simulations. S.H.E.R. and J.W.Y. analyzed the data. S.H.E.R. wrote the manuscript. All authors discussed the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All the 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.

Stay Connected to Science Advances

Navigate This Article