## Abstract

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.

## INTRODUCTION

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.

## RESULTS

### Microrheology

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, *T*_{K} = 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.

### 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*/∣*T* − *T _{mc}*∣

^{Δ}via

*T*= ∣

*T*−

*T*∣ is the distance from a critical delocalization temperature

_{mc}*T*,

_{mc}*T*=

*T*, scale invariance requires

_{mc}*F*∝

*U*, in which

^{q}*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

*T*and the critical exponent

_{mc}*q*. This technique is based on the successive slope of the velocity-force curves, defined as

*m*= ∂ ln

*F*/∂ ln

*U*: At

*T*=

*T*, the curvature of the

_{mc}*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*=

*T*, the slope becomes

_{mc}*m*=

*q*. From Fig. 1B, the curvature 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*=

*T*, a correction to scaling of the form

_{mc}*F*=

*U*(κ

^{q}_{1}+ κ

_{2}

*U*

^{ω/z}) is required, in which ω is the leading correction-to-scaling exponent, and

*z*is the dynamic exponent. Therefore, at

*T*=

*T*, the successive slope becomes

_{mc}Because the critical temperature *T _{mc}* 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.

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, *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 *F*/δ*T*^{Γ} versus *U*/δ*T*^{Δ} 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 *T _{mc}* = 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

*T*= 0.36875 ± 0.0062, which is approximately 23% higher than the thermodynamic glass transition temperature, at

_{mc}*T*= 0.3. The delocalization temperature

_{K}*T*is greater than the thermodynamic glass transition temperature

_{mc}*T*, a result consistent with the findings of Habdas

_{K}*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.,

*F*∝

*U*

^{1}, 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

*T*is approached, in agreement with a report by Williams and Evans (

_{mc}*27*).

Equation 1 suggests that the inherent friction, *F*_{inh} = *F*(*U* → 0, *T*), in the system scales as *F*_{inh} ∝∣*T* − *T _{mc}*∣

^{Γ}. 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*

_{∞}∝ (

*p*−

*p*

_{c})

^{β}. 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.

### Relaxations

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, *C _{F}*(τ) = 〈

*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.

### Time-force superposition

In Fig. 2 (B and D), *C _{F}*(τ) 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)*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 *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, *C _{F}*(τ, τ

_{w}), in which τ

_{w}is the waiting time. In agreement with findings for sheared glasses, we found that

*C*(τ, τ

_{F}_{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 *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, *C _{F}*(τ) ≃

*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

*C*(τ). 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.

_{F}### Macrorheology

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 *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 *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 *T _{c}* = 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.

### 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*a* is the radius of the probe particle. A scaling ansatz for microviscosity can be obtained from Eq. 1

A similar scaling ansatz can also be written for viscosity if η_{m}, *T _{mc}*, and

*U*are replaced by η,

*T*, and

_{c}_{m}and η, we define an effective shear rate

*D*= 2

*a*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

*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.

### 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} is the width of each distribution. The vertical dashed line shows the limit of macrorheology given by a delta function, *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.

## CONCLUSION AND DISCUSSION

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.

## METHODS

### Simulations

We performed large-scale 3D simulations of a 80 : 20 binary mixture (Kob-Andersen model) interacting via a smoothed LJ potential*i*, *j* = *A*, *B.* The parameters are represented in dimensionless LJ units in which the length, energy, mass, and time units are σ, ϵ, *m*, and *N _{A}* +

*N*)/

_{B}*V*= 1.2/σ

^{3}. The time step of integration is Δ

*t*= 0.01τ, and ϵ

*and σ*

_{ij}*are set to ϵ*

_{ij}*= 1.0ϵ, ϵ*

_{AA}*= 1.5ϵ, ϵ*

_{AB}*= 0.5ϵ, σ*

_{BB}*= 1.0σ, σ*

_{AA}*= 0.88σ, and σ*

_{AB}*= 0.8σ. Regardless of the particle type, the mass is*

_{BB}*m*= 1. The coefficients

*c*

_{0ij}and

*c*

_{2ij}are determined via

*U*(

_{ij}*r*) = 0 and

_{cij}*dU*/

_{ij}*dr*∣

_{r=rcij}= 0. The cutoff,

*r*is set to

_{cij},*r*= 2.5σ

_{cij}*, and the total number of particles is*

_{ij}*N*= 4000. The thermodynamic glass transition temperature for this system, the Kauzmann temperature, is

*T*

_{K}= 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 *U _{ij}*(

*r*) =

*A*[1 + cos (π

*r*/

*r*)] is applied to all pairs of interacting particles within the cutoff,

_{c}*r*, 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

_{c}*T*= 0.466ϵ/

*k*

_{B}. 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.

### Microrheology

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 potential

Unless otherwise stated, Δ = 2.0σ and *r _{c}* = 2

^{1/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*= −

**F**⋅

**n**where

**n**=

**U**/∣

**U**∣. Here,

**F**= ∑

_{j}**f**

*, where*

_{j}**f**

*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.*

_{j}## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/29/eaba8766/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).