Universal mechanism of shear thinning in supercooled liquids

Soft glassy materials experience a significant reduction in viscosity $\eta$ when subjected to shear flow, known as shear thinning. This phenomenon is characterized by a power-law scaling of $\eta$ with the shear rate $\dot{\gamma}$, $\eta \propto \dot{\gamma}^{-\nu}$, where the exponent $\nu$ is typically around $0.7$ to $0.8$ across different materials. Two decades ago, the mode coupling theory (MCT) suggested that shear thinning occurs due to the advection. However, it predicts too large $\nu = 1$ (>$0.7$ to $0.8$) and overestimates the onset shear rate by orders of magnitude. Recently, it was claimed that a minute distortion of the particle configuration is responsible for shear thinning. Here we extend the MCT to include the distortion, and find that both advection and distortion contribute to shear thinning, but the latter is dominant. Our formulation works quantitatively for several different glass formers. We explain why shear thinning is universal for many glassy materials.


INTRODUCTION
Nonlinear rheology is commonly observed in complex fluids and soft materials [1].In particular, supercooled glassy liquids exhibit a significant reduction in viscosity η and relaxation time τ α when subjected to fast shear flow, a phenomenon known as shear thinning.Both experiments [2][3][4] and simulations [5][6][7][8] have shown that η and τ α follow a power-law scaling with the shear rate γ as η, τ α ∝ γ−ν , where the exponent ν remains consistent across different materials, typically around 0.7 to 0.8.Understanding shear thinning is crucial not only for the manufacturing and processing of materials but also for broader physical phenomena such as volcanoes and earthquakes.However, we have yet to understand the mechanism responsible for this nonlinear flow and its universal nature.Various theories have been proposed in the past, including the soft glassy rheology theory [9], the shear transformation zone theory [10], and the elastoplastic model [11].
Among these theories, the mode coupling theory (MCT) is a first-principles theory, which was originally developed to explain the equilibrium dynamics of supercooled liquids near the glass transition point [12].The theory describes the slow glassy dynamics in terms of the caging effect; particles are trapped in the cages formed by their neighbors until the structure undergoes reconfiguration at the equilibrium relaxation time τ α0 , which diverges at the dynamical transition point in the mean-field limit.Note that we denote equilibrium values by the subscript 0 throughout this article.In later years, the MCT has been generalized to sheared liquids [13,14].The sheared MCT explains that advection induced by the shear flow breaks the cages and accelerates the dynamics.Shear thinning begins when the timescale of shear, γ−1 , * hideyuki.mizuno@phys.c.u-tokyo.ac.jp becomes comparable to τ α0 , that is, the onset shear rate is γc ∼ τ −1 α0 .The theory then predicts that as the shear rate further increases, τ α decreases as γ−1 and thus the thinning exponent ν = 1.
Although the sheared MCT provides a qualitative explanation for the reduction in relaxation time and viscosity, its prediction of ν = 1 is larger than the observation of ν ∼ 0.7 to 0.8 in many previous works [2][3][4][5][6][7][8].Moreover, the theory overestimates the values of the onset shear rate γc , which have been found to be orders of magnitude smaller than the theoretical prediction τ −1 α0 , i.e., γc ≪ τ α0 , in experiments and simulations [2,7,15,16].These discrepancies between the theory and the observations have remained unaddressed for more than two decades.
Recently, Furukawa [17][18][19] has proposed a semimicroscopic theory to explain the shear thinning in supercooled liquids, which is distinct from the advection scenario of the sheared MCT.The theory claims that anisotropic distortion of the particles' configuration due to shear flow, rather than advection, is responsible for shear thinning.Although this distortion is tiny in dense glassy fluids [20][21][22][23], it reduces effective density for fragile glass formers [17,18] or effective activation energy for strong glass formers [19], which, in turn, induces a drastic acceleration of the dynamics and causes the shear thinning.Based on this distortion scenario, Furukawa succeeded in quantitatively explaining the observed small thinning exponent, ν < 1, and the small onset rate, γc ≪ τ −1 α0 , for both fragile and strong glass formers.Questions naturally arise within us.(i) Does the advection scenario of the sheared MCT fail to explain the shear thinning in supercooled liquids?(ii) Does the distortion scenario work universally in different types of supercooled liquids?If so, (iii) can we renovate the sheared MCT by integrating the distortion effect into the theory and reconciling the theory with the observations?
To answer these three questions, we must first assess the validity of the advection scenario of the sheared MCT.To do this, we need a model system that can serve as an ideal fluid for testing the mean-field theory of the glass transition.The Gaussian core model (GCM) is a promising candidate because its slow dynamics are better described by the equilibrium MCT than any other glass-forming liquid [24][25][26][27].Firstly, the GCM is a clean, glassy model that does not require size dispersity [24][25][26].The monatomic GCM exhibits slow glassy dynamics close to the dynamical transition point T c without being affected by unwanted crystallization.Secondly, the equilibrium relaxation time of the GCM follows the MCT power-law scaling, over a wider temperature window than other models of glass formers.The agreement of the exponent γ ≃ 2.7 with the MCT is quantitative.Even the transition temperature T c , routinely used as a fitting parameter, agrees quantitatively with the MCT prediction.Thirdly, the violation of the Stokes-Einstein (SE) law is very weak, and the diffusion constant D 0 is proportional to τ −1 α0 , which is again consistent with the MCT prediction.Lastly, although a dramatic increase in dynamic heterogeneities accompanies the slow dynamics, the statistics of particles' displacements remain nearly Gaussian [27], and the growth of the dynamical heterogeneities is explained by the inhomogeneous MCT [28].This is in stark contrast with other glass formers, where the separation of fastand slow-moving clusters of particles characterizes the dynamical heterogeneities [29].Therefore, if the sheared MCT has any prediction regarding the shear thinning, the GCM should be the first model to be compared with that.
In addition to the GCM, we investigate canonical glass formers such as the Kob-Andersen (KA) model [30], the soft sphere (SS) model [31], and the van Beest-Kramervan Santen (BKS) model [32].The KA and SS models are typical fragile glass formers, while the BKS model mimics the silica melt, a representative strong glass former.We find that the GCM and these different types of supercooled liquids share similar scaling laws in τ α , η ∝ γ−ν with ν ∼ 0.7 (< 1) and γc ∝ τ −δ α0 ≪ τ α0 with δ ∼ 1.4 (> 1).This result indicates that the mechanism of shear thinning is universal, and it can not be explained by the advection scenario of the sheared MCT.
In particular, the GCM does not adhere to the advection scenario, which dictates that the current sheared MCT fails to explain the shear thinning.We resolve this conundrum by incorporating the distortion effect into the diverging relaxation time and viscosity that the MCT prescribes.Our analysis of the resulting equation reveals that the thinning exponents of ν and δ can be formulated by simple equations, which produce values of ν ∼ 0.7 and δ ∼ 1.4.We also extend the schematic model of the sheared MCT, proposed by Fuchs and Cates [33], to account for the distortion effect.Our theoretical and numerical results resolve long-standing inconsistencies between the theory and the observations in experiments GCM, γ/ γc (advection scenario).The data for the SS model are extracted from Ref. [18], the 2DSS model from Ref. [17], and the BKS model from Ref. [19].and simulations, and establish a universal mechanism of the shear thinning in supercooled liquids.

Numerical observations
We perform molecular dynamics (MD) simulations on the GCM which is subjected to shear flow in three spa-tial dimensions.The density is fixed at ρ = 2.0 where the dynamical transition point has been estimated as T c ≃ 2.68 × 10 −6 [24,26].We study a range of temperatures T near T c , so that our simulations explore supercooled states close to the dynamical transition.The shear rate γ is controlled over a wide range to cover the Newtonian to the strongly nonlinear regimes.From the MD trajectory data, we measure the relaxation time τ α , the diffusion constant D, and the viscosity η as a function of T and γ.For detailed information on MD simulations and calculations of τ α , D, and η, please see Methods.
In addition to the GCM, we conduct MD simulations on the KA model under shear flow and measure τ α and η.We also extract available data on the SS model from Ref. [18], the two-dimensional SS (2DSS) model from Ref. [17], and the BKS model from Ref. [19].Details of MD simulations on the KA model and system descriptions of the SS, 2DSS, and BKS models are provided in Methods.
Figure 1 presents the data on all the studied systems together.We provide τ α for the GCM and the KA model and η for the KA, SS, 2DSS, and BKS models.Note that the linear relation τ α ∝ η normally holds, as we confirm for the KA model in Methods and for the SS model in Ref. [34].Thus, τ α and η provide essentially the same information on dynamics, in general.However, we find that τ α ∝ η breaks at high shear rates in the GCM, as shown in Methods.This point requires further detailed investigation.In Fig. 1, we show τ α (not η) for the GCM.
In panel (a) of Fig. 1, we plot τ α /τ α0 and η/η 0 against γ/ γc .Here, τ α , η, and γ are normalized using the equilibrium values τ α0 and η 0 , and the onset shear rate γc , respectively, in order to compare different systems.We observe that all the systems studied exhibit similar dependences on γ, which are not proportional to γ−1 (red line), but rather proportional to γ−ν with ν ∼ 0.7 (blue line).We particularly emphasize that the GCM does not follow the ∝ γ−1 dependence of the advection scenario.For the GCM, we also plot data on D −1 /D −1 0 which are indistinguishable from those on τ α /τ α0 .This agreement demonstrates that the SE law in the form of τ α ∝ D −1 holds throughout the shear thinning regime, not just in equilibrium states [24,26].Thus, both the structural relaxation and the diffusion dynamics of the GCM do not follow the ∝ γ−1 of the advection scenario.
In panel (b) of Fig. 1, we present the onset shear rate γc as a function of τ α0 or η 0 .For all the systems studied, we observe that γc is much smaller than τ −1 α0 or η −1 0 (red line), and it follows γc ∝ τ −δ α0 or ∝ η −δ 0 with δ ∼ 1.4 (blue line).In particular, the GCM does not follow γc ∼ τ −1 α0 of the advection scenario.We thus conclude that the current sheared MCT fails to explain the shear thinning in supercooled liquids.In contrast, the GCM and the other systems share similar scaling behaviors of τ α ∝ γ−ν with ν ∼ 0.7 and γc ∝ τ −δ α0 (≪ τ −1 α0 ) with δ ∼ 1.4.This result suggests a universal mechanism of the shear thinning in supercooled liquids.TABLE I. Dynamical transition temperature Tc, mode coupling theory (MCT) scaling exponent γ, and predicted thinning exponents ν and δ.We present data for the Gaussian core model (GCM), the Kob-Andersen (KA) model, the soft sphere (SS) model, the two-dimensional SS (2DSS) model, and the van Beest-Kramer-van Santen (BKS) model.The values of Tc and γ are obtained from published papers cited in References.For 2DSS, we determine Tc and γ by fitting the MCT power-law scaling to the data on η0 versus T published in Ref. [17].After obtaining γ, we calculate ν and δ using the formulas ν = γ/(γ + 1) and δ = (γ + 1)/γ, respectively.As the advection scenario of the sheared MCT fails to explain the shear thinning, we turn our attention to the distortion scenario proposed by Furukawa [17][18][19].All the systems studied in Fig. 1, including not only the GCM but also the KA, SS, 2DSS, and BKS models, follow the MCT power-law scaling, Eq. ( 1), in the temperature regime above the critical temperature T c .We summarize the values of T c and γ in Table I, which are sourced from published papers.Note that γ in the current systems lies between 2 and 3.This fact motivates us to explain a universal mechanism behind shear thinning by incorporating the distortion effect into the MCT power-law scaling.

GCM
The shear flow distorts the particles' configuration in an anisotropic manner, which causes the effective density ρ eff to increase along the compression axis and decrease along the decompression axis.Since the distortion occurs on the timescale of τ α , the variation of ρ eff due to the shear is characterized by the strain γτ α .Thus, assuming γτ α ≪ 1 (we confirm that γτ α is at most 10 −1 for the GCM and the KA model), ρ eff can be described as where ρ is the density in the unsheared (equilibrium) state, and b ρ = ∂ρ eff /∂( γτ α )| γτα=0 can be positive or negative, depending on the direction of compression or decompression.
We next consider how the variation in ρ eff impacts the relaxation time τ α .In the case of the GCM, since the dynamics accelerate as the density increases [24,26], the direction along the compression axis, where b ρ is positive and ρ eff increases, contributes to the shear thinning.Although the dynamics become slow in the other direction (along the decompression axis), this does not prevent the shear thinning because an acceleration along the compression axis leads to a significant acceleration in overall dynamics.On the other hand, for the KA, SS, 2DSS, and BKS models, the dynamics speed up with decreasing the density, and the direction of the negative b ρ along the decompression axis causes the shear thinning.This behavior is contrary to that of the GCM.However, in both cases, a minute but finite variation in ρ eff commonly plays a crucial role in the shear thinning process.
Focusing on the direction of the positive b ρ for the GCM or that of the negative b ρ for the KA, SS, 2DSS, and BKS models, we can proceed with the following formulations for τ α (T, γ).Recall that τ α0 (T ) of the unsheared (equilibrium) system follows the MCT power-law scaling, Eq. ( 1), close to T c ; τ α0 ∝ [T − T c (ρ)] −γ where T c is a function of the density ρ.We assume that this power-law scaling remains valid under shear by replacing ρ in T c by ρ eff , i.e., In addition, applying Eq. ( 2) for ρ eff , we can approximate T c (ρ eff ) as where is a positive constant regardless of the system, since T c (ρ eff ) is a decreasing function of ρ eff for the GCM [24,26] whereas it is an increasing function of ρ eff for the KA, SS, 2DSS, and BKS models.Finally, using Eq. ( 4) for T c (ρ eff ) in Eq. ( 3), we arrive at a self-consistent equation for τ α (T, γ); By using Eq. ( 5), we can make predictions for the onset shear rate and the thinning scaling as below.Since b T γτ α0 becomes comparable to In addition, once the shear thinning builds up, b T γτ α ≫ T − T c (ρ) holds.Thus, Eq. ( 5) leads to τ α ∝ (b T γτ α ) −γ , giving a thinning scaling of Note that for strong glass formers like the BKS model, we need to consider the activation energy E eff instead of the density ρ eff in the above formulations [19].However, by replacing ρ eff with E eff , we arrive at the same selfconsistent equation for τ α (T, γ), Eq (5).This results in obtaining the same formulations as in Eqs. ( 6) and (7) for the strong glass formers.
These expressions are applicable to any system that remains above the dynamical transition temperature T c and follows the MCT power-law scaling given by Eq. (1).By substituting specific values of γ into these expressions, we can obtain values for ν and δ, which are summarized in Table I.For the present systems, we have values of γ ranging from 2 to 3, resulting in ν ∼ 0.7 and δ ∼ 1.4, which are quantitatively consistent with the observations in Fig. 1.
Therefore, we conclude that the distortion scenario universally works in different types of supercooled liquids, including the GCM and the fragile and strong glass formers.The exponents, ν = γ/(γ +1) and δ = (γ +1)/γ, are determined by the MCT exponent γ.This means that the power-law scaling in the shear thinning comes from the equilibrium MCT critical scaling near T c .The present systems show similar thinning exponents ν (∼ 0.7) and δ (∼ 1.4), which are generated by similar values of γ (∼ 2 to 3).

Schematic model of the sheared MCT
As we have seen so far, the distortion scenario accompanied by the MCT power-law scaling is successful in explaining the nontrivial values of exponents ν ∼ 0.7 and δ ∼ 1.4 in observations, through ν = γ/(γ + 1) and δ = (γ + 1)/γ.In the next step, we will integrate the distortion mechanism into the current sheared MCT to renovate the theory.
For this goal, we shall consider the schematic version of the MCT which drops wavenumber (k) dependences [12].The schematic MCT for unsheared (equilibrium) liquids, which is also known as the Leutheusser equation [41], has the mathematically same form as the k-dependent full MCT and preserves key characteristics of nontrivial slow dynamics and the dynamical transition, such as the power-law divergence of the relaxation time [42,43].The schematic MCT has been extended to sheared liquids [33,44,45], which again retains consequences of the k-dependent full sheared MCT.
We start with the sheared model proposed by Fuchs and Cates [33], with the memory kernel, where φ(t) represents a normalized intermediate scattering function, and the dot denotes the time derivative.λ is a parameter that contains information on the static structure factor and the temperature.In the equilibrium states with γ = 0, this model predicts the MCT powerlaw scaling, with γ ≃ 1.76, and the dynamical transition at λ c = 4.The term ( γt) 2 in the denominator in m(t) of Eq. ( 9) accounts for the advection effect by the shear flow.The fact that γ is scaled by t demonstrates that the advection and its resultant decoupling of the nonlinear coupling of density fields are responsible for shear thinning.As does the k-dependent full MCT, the model predicts τ α ∝ γ−1 and γc ∝ τ −1 α0 , i.e., values of exponents ν = 1 and δ = 1 (see solid lines of b λ = 0 in Fig. 2).These predictions contradict the numerical observations of ν < 1 and δ > 1 in Fig. 1, which are, however, correctly captured by the distortion scenario [17][18][19], as we have demonstrated in the previous section.
In the equilibrium MCT, the static structure factor S(|k|) (or λ in the schematic version) is an essential input parameter (k is wavevector, and |k| = k).In the sheared systems, S(|k|) (or λ) is distorted and replaced by a nonequilibrium function S NE (k) [20][21][22][23] (or λ NE ).So far, the sheared MCT has never taken S NE (k) into account based on the observation that the distortion of S(|k|) to S NE (k) is very small [46,47].However, we now understand from the distortion scenario that this tiny distortion is surely responsible for shear thinning and needs to be addressed in the theory.
Here, we propose to introduce the distortion effect into the schematic MCT model, Eq. ( 9), by modifying λ to λ NE as follows.The procedure is the same in formulating Eq. ( 2).The distortion occurs on the timescale of the structural relaxation time τ α , and the density field experiences the strain γτ α .Thus, assuming γτ α ≪ 1, the distorted parameter λ NE can be expressed as where b λ = −∂λ NE /∂( γτ α )| γτα=0 (> 0) is a model parameter which quantifies sensitivity to the shear flow.
In the present systems of the GCM and the KA, SS, 2DSS, BKS models, the situations all relate to large values of b λ .In such cases, the primary factor affecting shear thinning is distortion rather than advection.As a result, we update the sheared MCT to include the distortion effect, which resolves long-standing inconsistencies between the theory and the observations in experiments and simulations.

CONCLUSIONS
We have addressed questions (i) to (iii) raised in the Introduction.Firstly, we have observed that different types of systems, namely the GCM, and the KA, SS, and BKS models, exhibit similar scalings of τ α , η ∝ γ−ν with ν ∼ 0.7 (< 1) and γc ∝ τ −δ α0 with δ ∼ 1.4 (> 1), as in Fig. 1.The GCM does not follow τ α ∝ γ−1 and γc ∝ τ −1 α0 of the advection scenario, which dictates that (i) the current sheared MCT fails to explain the shear thinning.Next, we used the distortion scenario accompanied by the MCT power-law scaling and formulated the thinning exponents, ν = γ/(γ + 1) in Eq. ( 7) and δ = (γ + 1)/γ in Eq. ( 6), in terms of the MCT exponent γ.These formulations provide quantitatively correct values of ν ∼ 0.7 and δ ∼ 1.4 in the observations, thus concluding that (ii) the distortion scenario works universally in the GCM and the fragile and strong glass formers.Finally, we integrated the distortion effect into the schematic MCT model as in Eq. ( 12), which explains ν = γ/(γ + 1) and δ = (γ + 1)/γ.Consequently, (iii) we renovated the sheared MCT by accounting for the distortion effect.Our numerical and theoretical results (i) to (iii) have resolved the long-standing discrepancies between the theory and the observations in experiments and simulations, establishing a universal mechanism of shear thinning in supercooled liquids.
The thinning exponents ν and δ are determined by the MCT exponent γ.This indicates that the power-law scalings observed in shear thinning originate from the criticality of the equilibrium MCT near the dynamical transition point T c .All the systems studied in this work exhibit similar shear-rate dependences for τ α or η, which is due to their similar values of γ, ranging from 2 to 3. It would be interesting in future research to investigate systems with values of γ that differ significantly from this range.For instance, the harmonic spheres can display γ ≃ 5.3 at high packing fractions above ϕ = 0.8 [48], resulting in ν ≃ 0.84 and δ ≃ 1.19.
On the other hand, although macroscopic observables (τ α and η) follow similar shear-rate dependences across different systems, microscopic dynamics are ex-pected to be quite different.The equilibrium dynamics of the GCM differ significantly from those of typical liquids with short-ranged, harshly-repulsive potentials like the KA model [27].In the KA model, dynamics are described by the caging mechanism with hopping motions between local cages, whereas the GCM exhibits rather continuous motions that are not characterized by the standard caging mechanism.The most recent work [49] reported that the GCM also exhibits the caging dynamics at low densities, and upon increasing the density, a smooth variation occurs towards the non-caging dynamics.In addition, it was reported that dynamics are very different between fragile glass formers (SS model) and strong glass formers (BKS model) [50].Therefore, one would expect that microscopic dynamics under shear flow differ significantly between the GCM and the fragile and strong glass formers.
The present work focuses on the temperature regime above the dynamical transition temperature T c .In this regime, the shear thinning is closely related to the equilibrium MCT criticality.On the other hand, we expect a distinct behavior below T c .At the mean-field level, the equilibrium dynamics transition from non-activation to activation as the temperature decreases across T c .In finite dimensions, non-mean-field effects disrupt this transition, but we can still observe its remnants as a dynamical crossover in the KA model [51,52].In the future, it would be interesting to explore the nonlinear rheology below the dynamical transition.
Finally, it is commonly accepted that the viscosity is proportional to the relaxation time as η ∝ τ α (as shown in Methods for the KA model); the relaxation time is responsible for controlling the viscosity in glass-forming liquids.However, as shown in Methods, we have observed that this relationship does not apply in the GCM at high shear rates.This suggests that the shear modulus, measured as G = η/τ α , is dependent on the shear rate γ; in the GCM, as γ increases, so does G.Further analysis is required to investigate this matter in the future.

MD simulations on GCM subjected to shear flow
We perform MD simulations on the mono-disperse GCM in three spatial dimensions [24][25][26][27].The particles interact via the potential, where ǫ and σ characterize energy and length scales, respectively.The interaction is truncated at r = 5σ.The mass of particles is m.We use σ, ǫ/k B (k B is Boltzmann constant), and τ = (mσ 2 /ǫ) 1/2 as units of length, temperature, and time, respectively.The number density is fixed at ρ = N/V = 2.0, where N = 4000 is the number of particles and V is the system volume.At ρ = 2.0, the dynamical transition temperature estimated by the standard power-law fitting for τ α0 is T c ≃ 2.68×10 −6 [24, 26].
To explore supercooled states, we study various temperatures ranging from T × 10 6 = 7.0 to 2.9 which is close enough to T c .
After the system was equilibrated at each temperature T , we applied a steady shear flow to drive the system into a nonequilibrium state [5][6][7][8].We integrated the SLLOD equations using the Lees-Edwards boundary condition, with the Nosé-Hoover thermostat to maintain the temperature [53].To cover the Newtonian to the strongly nonlinear regimes, we control the shear rate γ over a wide range from γ = 10 −8 to 10 −3 .Here we set the x axis along the flow direction and the y axis along the velocity gradient direction.The mean velocity profile v is thus given as where e µ (µ = x, y, z) is the unit vector along the µ axis.We note that γ ∼ 5 × 10 −4 is high enough that the relaxation time τ α reaches the timescale of vibrations, the so-called Einstein period [54].

Self-intermediate scattering function and mean squared displacements of GCM
We employ two measurements to study the dynamics of particles; the self-intermediate scattering function and the mean squared displacements ∆r 2 (t) , where r i = (x i , y i , z i ) is the position of particle i, denotes the ensemble average, and we subtract from the total displacement of each particle, the contribution resulting from the advective transport by the mean shear flow, γ t 0 y i (s)dse x [5].components; (19) We observe that both F s (k, t) and ∆r 2 (t) show a drastic acceleration of the dynamics due to the shear flow.In addition, both data are isotropic, showing little dependence on x, y, and z directions or components, even at the highest shear rate γ = 10 −3 .are plotted as a function of (T − Tc)/Tc.The lines present the power-law scalings, τ −1 α0 ∝ (T − Tc) γ with γ ≃ 2.7 [Eq. ( 1)] and γc ∝ (T − Tc) γδ with γδ ≃ 3.7 [Eq.(21)].(c) γc is plotted against τα0.The blue line presents the scaling relation of γc ∝ τ −δ α0 with δ ≃ 1.37 [Eq.( 6)].

Relaxation time of GCM
From the relaxation behavior of F s (k, t), we calculate the relaxation time τ α as Figure 4(a) shows the shear-rate γ dependence of τ α (T, γ) for various temperatures T .The figure demonstrates the shear thinning behavior in the GCM, which is characterized by a power-law scaling τ α ∝ γ−ν with an exponent ν ≃ 0.73 (blue line).The value of ν ≃ 0.73 is obtained through ν = γ/(γ + 1) with γ = 2.7 as shown in Table I.
In Fig. 4(b), we plot γc as a function of (T − T c )/T c , and compare it to τ −1 α0 as the sheared MCT predicts γc1 ∼ τ −1 α0 .It is observed that γc is considerably (orders of magnitude) smaller than τ −1 α0 .Note that this figure also confirms the MCT power-law scaling, Eq. ( 1), close to the critical temperature T c (see the line for squares), in keeping with previous works [24,26].In addition, we display γc against τ α0 in Fig. 4(c), showing that the data are well fitted by γc ∝ τ −δ α0 with δ ≃ 1.37 (blue line).The value of δ ≃ 1.37 is obtained through δ = (γ + 1)/γ with γ = 2.7 as in Table I.This result suggests that γc follows a power-law scaling, with γδ ≃ 3.7.This scaling indeed works close to T c , as confirmed in Fig. 4(b) (see the line for circles).

Diffusion constant of GCM
The diffusion constant D is determined by observing the diffusive behavior of ∆r 2 (t) in the long-time limit, which can be quantified as We present data on the inverse diffusion constant D −1 (T, γ) in Fig. 6(a) and D −1 (T, γ)/D −1 0 (T ) in Fig. 6(b), which are counterparts of Figs.4(a) and 5 for τ α (T, γ), respectively.It is clear that D −1 (T, γ) follows the same power-law behavior as that of τ α (T, γ) in Eq. (22).In Fig. 1 in the main text, we plot data on D −1 /D −1 0 versus γ/ γc at T = 2.9 × 10 −6 , which are indistinguishable to those on τ α /τ α0 .These observations demonstrate that the SE law in the form of τ α ∝ D −1 holds throughout the shear thinning regime, not just in equilibrium states [24,26].

Viscosity of GCM
We measure the viscosity η as a function of T and γ.We calculate the shear stress σ xy as [55] where v i = (v ix , v iy , v iz ) is the velocity of particle i, r ij = (x ij , y ij , z ij ) denotes the vector r i − r j = (x i − x j , y i − y j , z i −z j ), and r ij = |r ij |.The viscosity is then obtained through η = σ xy / γ.
In many previous works [5][6][7][8][17][18][19], it has been assumed that the relaxation time controls the viscosity in glass-forming liquids, and that η ∝ τ α .This assumption has been confirmed for the KA model below and for the SS model in Ref. [34].The sheared MCT also formulates η ∝ τ α [13,14,33,46].Thus, it is considered that τ α and η provide essentially the same information on dynamics in supercooled liquids.However, this assump- tion does not hold true for the GCM at high shear rates.This result suggests that the shear modulus measured as G = η/τ α is dependent on γ; G increases as γ gets larger.Further detailed investigation is required in the future to understand this point better.

MD simulations on KA model subjected to shear flow
We perform MD simulations on a binary Lennard-Jones (LJ) mixture, the KA model, in three spatial dimensions [30] where α and β denote A or B, and the parameters are set to be ǫ AA = ǫ, ǫ AB = 1.5ǫAA , ǫ BB = 0.5ǫ AA , σ AA = σ, σ AB = 0.8σ AA , σ BB = 0.88σ AA .The interaction is truncated at r = 2.5σ αβ .We employ σ, ǫ/k B , and τ = (mσ 2 /ǫ) 1/2 as units of length, temperature, and time, respectively.The number density is fixed at ρ = N/V = 1.2, and the number of particles is N = N A + N B = 4000 with N A = 3200 and N B = 800.At ρ = 1.2, the standard power-law fitting for τ α0 estimates the dynamical transition temperature to be T c ≃ 0.435 [30,35,36].We study at various temperatures ranging from T = 0.8 to 0.45 close to T c .The shear rate γ is controlled over a wide range of γ = 10 −6 to 10 −1 .Note that γ ∼ 10 −1 is high enough that τ α reaches the timescale of vibrations (Einstein period) [54].
We analyze the KA model in the same way as we do for the GCM.At each temperature T , we measure the relaxation time τ α as a function of γ and identify the onset shear rate γc .We present results for the larger particles (A) below, but similar results were obtained for the smaller particles (B). Figure 9 presents τ α versus γ in (a), and the data on γc and τ −1 α0 in (b) and (c).Figure 10 then plots τ α (T, γ)/τ α0 (T ) as a function of γ/ γc , where we exclude data at γ = 10 −1 at which τ α reaches the timescale of vibrations.We also measure the viscosity η and present data on η(T, γ) and η(T, γ)/η 0 (T ) in Fig. 11.Furthermore, we compare the relaxation time and the viscosity in Fig. 12. Figures 9, 10 In Figs. 9 and 10, we can see that the relaxation time follows a power-law scaling of the form τ α ∝ γ−ν , where ν ≃ 0.71, and the onset shear rate γc exhibits the scaling behavior γc ∝ τ −δ α0 ≪ τ −1 α0 , with δ ≃ 1.42.The values of ν and δ are obtained using ν = γ/(γ + 1) and δ = (γ + 1)/γ with γ = 2.4 as shown in Table I.
In addition, Figures 11 and 12 show that η is proportional to τ α as η ∝ τ α , following the same scaling law as that of τ α ; η(T, γ) η 0 (T ) Differently from the case of the GCM, the linear relation η ∝ τ α is kept even at the high shear rates of γ/ γc 10 2 .This result indicates that the relaxation time controls the viscosity in the KA model, as we normally expect and suppose.Figure 1 in the main text presents data on τ α /τ α0 and η/η 0 versus γ/ γc at T = 0.45 in (a), and γc versus τ α0 in (b).

Figure 3 (FIG. 3 .
Figure 3(a) displays F s (k, t) for k = k max e x , k max e y , and k max e z , where k max ≃ 8.4 is the wavenumber at which the static structure factor takes a maximum.In Fig. 3(b), we report ∆r 2 (t) by separating x, y, and z