Universal mechanism of shear thinning in supercooled liquids

Nonlinear rheology is commonly observed in complex fluids and soft materials¹. 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 \(\dot{\gamma }\) as \(\eta ,\,{\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\), 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⁹, the shear transformation zone theory¹⁰, and the elastoplastic model¹¹.

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¹². 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, \({\dot{\gamma }}^{-1}\), becomes comparable to τ_α0, that is, the onset shear rate is \({\dot{\gamma }}_{c} \sim {\tau }_{\alpha 0}^{-1}\). The theory then predicts that as the shear rate further increases, τ_α decreases as \({\dot{\gamma }}^{-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 \({\dot{\gamma }}_{c}\), which have been found to be orders of magnitude smaller than the theoretical prediction \({\tau }_{\alpha 0}^{-1}\), i.e., \({\dot{\gamma }}_{c}\ll {\tau }_{\alpha 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 semi-microscopic 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¹⁹, 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, \({\dot{\gamma }}_{c}\ll {\tau }_{\alpha 0}^{-1}\), for both fragile and strong glass formers.

Questions naturally arise within us. (i) Does the advection scenario of the sheared MCT entirely 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₀ is proportional to \({\tau }_{\alpha 0}^{-1}\), 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²⁷, and the growth of the dynamical heterogeneities is explained by the inhomogeneous MCT²⁸. This is in stark contrast with other glass formers, where the separation of fast- and slow-moving clusters of particles characterizes the dynamical heterogeneities²⁹. 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³⁰, the soft sphere (SS) model³¹, and the van Beest-Kramer-van Santen (BKS) model³². 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 \({\tau }_{\alpha },\,\eta \propto {\dot{\gamma }}^{-\nu }\) with ν ~ 0.7 ( < 1) and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\ll {\tau }_{\alpha 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³³, to account for the distortion effect. Our theoretical and numerical results resolve long-standing inconsistencies between the theory and the observations in experiments 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 spatial dimensions. The density is fixed at ρ = 2.0 where the dynamical transition point has been estimated as T_c ≃ 2.68 × 10^−624,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 \(\dot{\gamma }\) 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 \(\dot{\gamma }\). 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. ¹⁸, the two-dimensional SS (2DSS) model from ref. ¹⁷, and the BKS model from ref. ¹⁹. 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. ³⁴. 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.

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. a Plots of τ_α/τ_α0, η/η₀, or \({D}^{-1}/{D}_{0}^{-1}\) are shown as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). The temperature is T = 2.9 × 10⁻⁶ (GCM), 0.45 (KA), 0.275 (SS), 0.577 (2DSS), and 0.511 (BKS), all above the dynamical transition temperature T_c. The black line represents \({\tau }_{\alpha }/{\tau }_{\alpha 0},\,\eta /{\eta }_{0},\,{D}^{-1}/{D}_{0}^{-1}=1\). The blue line represents \(\propto {\dot{\gamma }}^{-0.7}\), while the red line refers to \(\propto {\dot{\gamma }}^{-1}\) (advection scenario). The vertical dotted line indicates the onset shear rate \({\dot{\gamma }}_{c}\). b \({\dot{\gamma }}_{c}\) is plotted against τ_α0 or η₀. The blue line presents \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-1.4}\) or \({\eta }_{0}^{-1.4}\), while the red line refers to \({\dot{\gamma }}_{c}={\tau }_{\alpha 0}^{-1}\) or \({\eta }_{0}^{-1}\) (advection scenario). The data for the SS model are extracted from ref. ¹⁸, the 2DSS model from ref. ¹⁷, and the BKS model from ref. ¹⁹.

Full size image

In panel (a) of Fig. 1, we plot τ_α/τ_α0 and η/η₀ against \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Here, τ_α, η, and \(\dot{\gamma }\) are normalized using the equilibrium values τ_α0 and η₀, and the onset shear rate \({\dot{\gamma }}_{c}\), respectively, in order to compare different systems. We observe that all the systems studied exhibit similar dependences on \(\dot{\gamma }\), which are not proportional to \({\dot{\gamma }}^{-1}\) (red line), but rather proportional to \({\dot{\gamma }}^{-\nu }\) with ν ~ 0.7 (blue line). We particularly emphasize that the GCM does not follow the \(\propto {\dot{\gamma }}^{-1}\) dependence of the advection scenario. For the GCM, we also plot data on \({D}^{-1}/{D}_{0}^{-1}\) which are indistinguishable from those on τ_α/τ_α0. This agreement demonstrates that the SE law in the form of τ_α ∝ D⁻¹ 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 \(\propto {\dot{\gamma }}^{-1}\) of the advection scenario.

In panel (b) of Fig. 1, we present the onset shear rate \({\dot{\gamma }}_{c}\) as a function of τ_α0 or η₀. For all the systems studied, we observe that \({\dot{\gamma }}_{c}\) is much smaller than \({\tau }_{\alpha 0}^{-1}\) or \({\eta }_{0}^{-1}\) (red line), and it follows \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) or \(\propto {\eta }_{0}^{-\delta }\) with δ ~ 1.4 (blue line). In particular, the GCM does not follow \({\dot{\gamma }}_{c} \sim {\tau }_{\alpha 0}^{-1}\) 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 \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) with ν ~ 0.7 and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\,(\ll {\tau }_{\alpha 0}^{-1})\) with δ ~ 1.4. This result suggests a universal mechanism of the shear thinning in supercooled liquids.

Distortion scenario

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 1, 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.

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 \(\dot{\gamma }{\tau }_{\alpha }\). Thus, assuming \(\dot{\gamma }{\tau }_{\alpha }\ll 1\) (we confirm that \(\dot{\gamma }{\tau }_{\alpha }\) is at most 10⁻¹ for the GCM and the KA model), ρ_eff can be described as

where ρ is the density in the unsheared (equilibrium) state, and \({b}_{\rho }={\left.\partial {\rho }_{{{{{{{{\rm{eff}}}}}}}}}/\partial (\dot{\gamma }{\tau }_{\alpha })\right\vert }_{\dot{\gamma }{\tau }_{\alpha } = 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 \({\tau }_{\alpha }(T,\dot{\gamma })\). Recall that τ_α0(T) of the unsheared (equilibrium) system follows the MCT power-law scaling, Eq. (1), close to T_c; \({\tau }_{\alpha 0}\propto {[T-{T}_{c}(\rho )]}^{-\gamma }\) 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 \({b}_{T}=-{\left.\partial {T}_{c}/\partial (\dot{\gamma }{\tau }_{\alpha })\right\vert }_{\dot{\gamma }{\tau }_{\alpha } = 0}=-{b}_{\rho }{\left.d{T}_{c}/d{\rho }_{{{{{{{{\rm{eff}}}}}}}}}\right\vert }_{{\rho }_{{{{{{{{\rm{eff}}}}}}}}} = \rho }\) 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 \({\tau }_{\alpha }(T,\dot{\gamma })\);

By using Eq. (5), we can make predictions for the onset shear rate and the thinning scaling as below. Since \({b}_{T}\dot{\gamma }{\tau }_{\alpha 0}\) becomes comparable to \(T-{T}_{c}(\rho )\propto {\tau }_{\alpha 0}^{-1/\gamma }\) at the onset \(\dot{\gamma }={\dot{\gamma }}_{c}\), we obtain

In addition, once the shear thinning builds up, \({b}_{T}\dot{\gamma }{\tau }_{\alpha } \, \gg T-{T}_{c}(\rho )\) holds. Thus, Eq. (5) leads to \({\tau }_{\alpha }\propto {({b}_{T}\dot{\gamma }{\tau }_{\alpha })}^{-\gamma }\), 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¹⁹. However, by replacing ρ_eff with E_eff, we arrive at the same self-consistent equation for \({\tau }_{\alpha }(T,\dot{\gamma })\), Eq (5). This results in obtaining the same formulations as in Eqs. (6) and (7) for the strong glass formers.

We thus derive expressions for the thinning exponents ν = γ/(γ + 1) in Eq. (7) and δ = (γ + 1)/γ in Eq. (6). 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 1. 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¹². The schematic MCT for unsheared (equilibrium) liquids, which is also known as the Leutheusser equation³⁵, 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^36,37. The schematic MCT has been extended to sheared liquids^33,38,39, which again retains the consequences of the k-dependent full-sheared MCT.

We start with the sheared \({F}_{2}^{(\dot{\gamma })}\) model proposed by Fuchs and Cates³³,

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 \(\dot{\gamma }=0\), this model predicts the MCT power-law scaling,

with γ ≃ 1.76, and the dynamical transition at λ_c = 4. The term \({(\dot{\gamma }t)}^{2}\) in the denominator in m(t) of Eq. (9) accounts for the advection effect by the shear flow. The fact that \(\dot{\gamma }\) 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 \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-1}\) and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-1}\), 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.

We present predictions for different b_λ values; b_λ = 0 (solid lines), 3 × 10⁻⁴ (dashed lines), and 10⁻² (dotted lines). Note that b_λ = 0 corresponds to the original MCT model with advection only. a τ_α is plotted as a function of \(\dot{\gamma }\). λ is fixed at λ_c − λ = 2 × 10⁻⁶ (λ_c = 4). The arrow indicates \(\dot{\gamma }={\tau }_{\alpha 0}^{-1}\). The red and blue lines respectively indicate \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-1}\) due to advection and \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) with ν = γ/(γ + 1) ≃ 0.64 ( < 1) due to distortion. b \({\dot{\gamma }}_{c}\) is plotted against τ_α0. The red and blue lines respectively indicate \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-1}\) due to advection and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ = (γ + 1)/γ ≃ 1.57 ( > 1) due to distortion.

Full size image

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^40,41. 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 \(\dot{\gamma }{\tau }_{\alpha }\). Thus, assuming \(\dot{\gamma }{\tau }_{\alpha }\ll 1\), the distorted parameter λ_NE can be expressed as

where \({b}_{\lambda }={\left.-\partial {\lambda }_{{{{{{{{\rm{NE}}}}}}}}}/\partial (\dot{\gamma }{\tau }_{\alpha })\right\vert }_{\dot{\gamma }{\tau }_{\alpha } = 0}\,( > 0)\) is a model parameter that quantifies sensitivity to the shear flow. Replacing λ by λ_NE in Eq. (9) while keeping the advection effect, we have

which accounts for the distortion effect in addition to the advection effect. In the case of b_λ = 0, the model reduces to the original model, Eq. (9), with the advection effect only. By setting b_λ > 0, the distortion effect is introduced, and it is increased by increasing b_λ.

Let us first discuss the distortion effect solely by considering the MCT model without the advection;

Since λ is replaced by \({\lambda }_{{{{{{{{\rm{NE}}}}}}}}}=\lambda -{b}_{\lambda }\dot{\gamma }{\tau }_{\alpha }\) [Eq. (11)], the MCT power-law scaling, Eq. (10), for τ_α0(λ) (in the absence of \(\dot{\gamma }\)) transforms to a self-consistent equation for \({\tau }_{\alpha }(\lambda ,\dot{\gamma })\);

This equation is essentially the same as Eq. (5). Solving Eq. (14), we obtain the same formulations of δ = (γ + 1)/γ in Eq. (6) and ν = γ/(γ + 1) in Eq. (7). We thus conclude that the distortion effect is correctly embedded in the MCT model, Eq. (12). The MCT exponent γ ≃ 1.76 provides specific values of δ = (γ + 1)/γ ≃ 1.57 ( > 1) and ν = γ/(γ + 1) ≃ 0.64 ( < 1) in this MCT model.

Figure 2 presents numerical solutions of the MCT model, Eq. (12), for three different values of b_λ: 0 (solid lines), 3 × 10⁻⁴ (dashed lines), and 10⁻² (dotted lines). Panel (a) shows τ_α as a function of \(\dot{\gamma }\) for λ_c − λ = 2 × 10⁻⁶, while panel (b) shows \({\dot{\gamma }}_{c}\) as a function of τ_α0. When b_λ = 0, the advection scenario is produced, with \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-1}\) and \({\dot{\gamma }}_{c} \sim {\tau }_{\alpha 0}^{-1}\). However, when b_λ = 10⁻², the distortion scenario is produced, with \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) and \({\dot{\gamma }}_{c} \sim {\tau }_{\alpha 0}^{-\delta }\) as described in Eqs. (7) and (6), respectively. At the intermediate value of b_λ = 3 × 10⁻⁴, the model predicts the crossover from the distortion-induced \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) to the advection-induced \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-1}\). The model also explains that a larger distortion effect with increasing b_λ results in a much slower \({\dot{\gamma }}_{c}\).

In the present systems of the GCM and the KA, SS, 2DSS, and 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.

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 \({\tau }_{\alpha },\,\eta \propto {\dot{\gamma }}^{-\nu }\) with ν ~ 0.7 ( < 1) and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ ~ 1.4 ( > 1), as in Fig. 1. The GCM does not follow \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-1}\) and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-1}\) of the advection scenario, which dictates that (i) the current sheared MCT entirely 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⁴², resulting in ν ≃ 0.84 and δ ≃ 1.19.

On the other hand, although macroscopic observables (τ_α and η) follow similar shear-rate dependencies across different systems, microscopic dynamics are expected 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²⁷. 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⁴³ 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)⁴⁴. Therefore, one would expect that microscopic dynamics under shear flow differ significantly among 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 dynamic crossover in the KA model^45,46. 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 \(\dot{\gamma }\); in the GCM, as \(\dot{\gamma }\) 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 \(\tau ={(m{\sigma }^{2}/\epsilon )}^{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⁶ = 10.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\(\acute{{{{{{{{\rm{e}}}}}}}}}\)-Hoover thermostat to maintain the temperature⁴⁷. To cover the Newtonian to the strongly nonlinear regimes, we control the shear rate \(\dot{\gamma }\) over a wide range from \(\dot{\gamma }=1{0}^{-8}\) to 10⁻³. 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 \(\dot{\gamma } \sim 5\times 1{0}^{-4}\) is high enough that the relaxation time τ_α reaches the timescale of vibrations, the so-called Einstein period⁴⁸.

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 F_s(k, t),

and the mean-squared displacements 〈Δr²(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, \(\dot{\gamma }\int\nolimits_{0}^{t}{y}_{i}(s)ds{{{{{{{{\bf{e}}}}}}}}}_{x}\)⁵.

Figure 3a displays F_s(k, t) for \({{{{{{{\bf{k}}}}}}}}={k}_{\max }{{{{{{{{\bf{e}}}}}}}}}_{x}\), \({k}_{\max }{{{{{{{{\bf{e}}}}}}}}}_{y}\), and \({k}_{\max }{{{{{{{{\bf{e}}}}}}}}}_{z}\), where \({k}_{\max }\simeq 8.4\) is the wavenumber at which the static structure factor takes a maximum. In Fig. 3b, we report 〈Δr²(t)〉 by separating x, y, and z components;

We observe that both F_s(k, t) and 〈Δr²(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 \(\dot{\gamma }=1{0}^{-3}\).

a \({F}_{s}({{{{{{{\bf{k}}}}}}}}={k}_{\max }{{{{{{{{\bf{e}}}}}}}}}_{\mu },t)\) and b 〈Δμ²(t)〉 are plotted as a function of t. The temperature is T × 10⁶ = 3.0. The red lines present values in the equilibrium state with \(\dot{\gamma }=0\), while black lines give values in the sheared states with different \(\dot{\gamma }\); from right to left, \(\dot{\gamma }=1{0}^{-7}\), 10⁻⁶, 10⁻⁵, 10⁻⁴, and 10⁻³. We plot data for μ = x, y, and z directions or components, represented by solid, dashed, and dotted lines, respectively.

Full size image

Relaxation time of GCM

From the relaxation behavior of F_s(k, t), we calculate the relaxation time τ_α as

Figure 4a shows the shear-rate \(\dot{\gamma }\) dependence of \({\tau }_{\alpha }(T,\dot{\gamma })\) for various temperatures T. The figure demonstrates the shear thinning behavior in the GCM, which is characterized by a power-law scaling \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) with an exponent ν ≃ 0.73 (blue line). The value of ν ≃ 0.73 is obtained through ν = γ/(γ + 1) with γ = 2.7 as shown in Table 1. We note that at the high shear rate \(\dot{\gamma }\gtrsim 5\times 1{0}^{-4}\), τ_α reaches the timescale of vibrations⁴⁸, and consequently it deviates from the scaling behavior of \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\).

a \({\tau }_{\alpha }(T,\dot{\gamma })\) is plotted as a function of \(\dot{\gamma }\). Symbols of different colors represent values at different T; from bottom to top, T × 10⁶ = 10.0 (yellow), 7.0 (purple), 5.0 (green), 4.0 (cyan), 3.4 (orange), 3.2 (blue), 3.0 (red), and 2.9 (black), all of which are above T_c × 10⁶ ≃ 2.68. Closed symbols indicate the equilibrium values τ_α0(T). For T × 10⁶ = 2.9, black line indicates \({\tau }_{\alpha }={\tau }_{\alpha 0}\propto {\dot{\gamma }}^{0}\), whereas the blue line indicates \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.73 [Eq. (7)]. b \({\dot{\gamma }}_{c}\) and \({\tau }_{\alpha 0}^{-1}\) are plotted as a function of (T − T_c)/T_c. The lines present the power-law scalings, \({\tau }_{\alpha 0}^{-1}\propto {(T-{T}_{c})}^{\gamma }\) with γ ≃ 2.7 [Eq. (1)] and \({\dot{\gamma }}_{c}\propto {(T-{T}_{c})}^{\gamma \delta }\) with γδ ≃ 3.7 [Eq. (21)]. c \({\dot{\gamma }}_{c}\) is plotted against τ_α0. The blue line presents the scaling relation of \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ ≃ 1.37 [Eq. (6)].

Full size image

Next, for each temperature T, we measure the onset shear rate \({\dot{\gamma }}_{c}\) at which the shear thinning starts. In Fig. 4b, we plot \({\dot{\gamma }}_{c}\) as a function of (T − T_c)/T_c, and compare it to \({\tau }_{\alpha 0}^{-1}\) as the sheared MCT predicts \({\dot{\gamma }}_{c1} \sim {\tau }_{\alpha 0}^{-1}\). It is observed that \({\dot{\gamma }}_{c}\) is considerably (orders of magnitude) smaller than \({\tau }_{\alpha 0}^{-1}\). 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 \({\dot{\gamma }}_{c}\) against τ_α0 in Fig. 4c, showing that the data are well fitted by \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ ≃ 1.37 (blue line). The value of δ ≃ 1.37 is obtained through δ = (γ + 1)/γ with γ = 2.7 as in Table 1. This result suggests that \({\dot{\gamma }}_{c}\) follows a power-law scaling,

with γδ ≃ 3.7. This scaling indeed works close to T_c, as confirmed in Fig. 4b (see the line for circles).

We then present a scaled plot of \({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\) versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) in Fig. 5. Note that in Fig. 5, we exclude data at the high shear rates \(\dot{\gamma }\ge 5\times 1{0}^{-4}\) at which τ_α reaches the timescale of vibrations. The data collapse onto a single curve regardless of temperature, which establishes

where ν ≃ 0.73, and \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ ≃ 1.37. Figure 1 in the main text presents data on τ_α/τ_α0 versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at the lowest T = 2.9 × 10⁻⁶ in (a), and \({\dot{\gamma }}_{c}\) versus τ_α0 in (b).

\({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\) is plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Symbols of different colors represent values at different T; T × 10⁶ = 10.0 (yellow), 7.0 (purple), 5.0 (green), 4.0 (cyan), 3.4 (orange), 3.2 (blue), 3.0 (red), and 2.9 (black). Black and blue lines indicate τ_α/τ_α0 = 1 and \({\tau }_{\alpha }/{\tau }_{\alpha 0}\propto {(\dot{\gamma }/{\dot{\gamma }}_{c})}^{-\nu }\) with ν ≃ 0.73, respectively. The vertical dotted line indicates the location of the onset shear rate \({\dot{\gamma }}_{c}\).

Full size image

Diffusion constant of GCM

The diffusion constant D is determined by observing the diffusive behavior of 〈Δr²(t)〉 in the long-time limit, which can be quantified as

We present data on the inverse diffusion constant \({D}^{-1}(T,\dot{\gamma })\) in Fig. 6a and \({D}^{-1}(T,\dot{\gamma })/{D}_{0}^{-1}(T)\) in Fig. 6b, which are counterparts of Figs. 4a and 5 for \({\tau }_{\alpha }(T,\dot{\gamma })\), respectively. It is clear that \({D}^{-1}(T,\dot{\gamma })\) follows the same power-law behavior as that of \({\tau }_{\alpha }(T,\dot{\gamma })\) in Eq. (22). In Fig. 1 in the main text, we plot data on \({D}^{-1}/{D}_{0}^{-1}\) versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at T = 2.9 × 10⁻⁶, which are indistinguishable to those on τ_α/τ_α0. These observations demonstrate that the SE law in the form of τ_α ∝ D⁻¹ holds throughout the shear thinning regime, not just in equilibrium states^24,26.

a \({D}^{-1}(T,\dot{\gamma })\) is plotted as a function of \(\dot{\gamma }\), and b \({D}^{-1}(T,\dot{\gamma })/{D}_{0}^{-1}(T)\) is plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Symbols of different colors represent values at different T; T × 10⁶ = 10.0 (yellow), 7.0 (purple), 5.0 (green), 4.0 (cyan), 3.4 (orange), 3.2 (blue), 3.0 (red), and 2.9 (black). In a, closed symbols indicate the equilibrium values \({D}_{0}^{-1}(T)\). Black line indicates \({D}^{-1}={D}_{0}^{-1}\propto {\dot{\gamma }}^{0}\), whereas blue line indicates \({D}^{-1}\propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.73. In b, the vertical dotted line indicates the location of \({\dot{\gamma }}_{c}\).

Full size image

Viscosity of GCM

We measure the viscosity η as a function of T and \(\dot{\gamma }\). We calculate the shear stress σ_xy as⁴⁹

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 \(\eta ={\sigma }_{xy}/\dot{\gamma }\).

We present data on the viscosity \(\eta (T,\dot{\gamma })\) in Fig. 7a and \(\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) in Fig. 7b, which are counterparts of Figs. 4a and 5 for \({\tau }_{\alpha }(T,\dot{\gamma })\), respectively. Note that the equilibrium values η₀(T) are obtained by averaging values of \(\eta (T,\dot{\gamma })\) in the Newtonian regime. Although data on τ_α/τ_α0 in Fig. 5 (and \({D}^{-1}/{D}_{0}^{-1}\) in Fig. 6b) are well fitted by \(\propto {(\dot{\gamma }/{\dot{\gamma }}_{c})}^{-\nu }\) with ν ≃ 0.73 as mentioned in Eq. (22), it does not work for η/η₀ as we can observe in Fig. 7b.

a \(\eta (T,\dot{\gamma })\) is plotted as a function of \(\dot{\gamma }\), and b \(\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) is plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Symbols of different colors represent values at different T; T × 10⁶ = 10.0 (yellow), 7.0 (purple), 5.0 (green), 4.0 (cyan), 3.4 (orange), 3.2 (blue), 3.0 (red), and 2.9 (black). In a, closed symbols indicate the equilibrium values η₀(T). Black line indicates \(\eta ={\eta }_{0}\propto {\dot{\gamma }}^{0}\), whereas blue line indicates \(\eta \propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.73. Note that the blue line, which is well fitted to data on τ_α and D⁻¹, does not work for η. In b, the vertical dotted line indicates the location of \({\dot{\gamma }}_{c}\).

Full size image

In addition, Fig. 8 shows a comparison between η/η₀ and τ_α/τ_α0 as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). At low shear rates of \(\dot{\gamma }/{\dot{\gamma }}_{c}\lesssim 1{0}^{2}\), we observe that η/η₀ coincides with τ_α/τ_α0, which confirms that η is proportional to τ_α. However, this linear relation η ∝ τ_α is systematically violated at high shear rates of \(\dot{\gamma }/{\dot{\gamma }}_{c}\gtrsim 1{0}^{2}\). These high shear rates correspond to the power-law scaling regime of \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\).

a \({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T),\,\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) and b \(\left[{\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\right]/{\dot{\gamma }}^{-\nu },\,\left[\eta (T,\dot{\gamma })/{\eta }_{0}(T)\right]/{\dot{\gamma }}^{-\nu }\) are plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Circles and squares represent values of τ_α and η, respectively. The temperature is T × 10⁶ = 2.9. The black line indicates τ_α/τ_α0, η/η₀ = 1, whereas the blue line indicates \({\tau }_{\alpha }/{\tau }_{\alpha 0},\,\eta /{\eta }_{0}\propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.73. The vertical dotted line indicates the location of \({\dot{\gamma }}_{c}\).

Full size image

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. ³⁴. The sheared MCT also formulates η ∝ τ_α^13,14,33,40. Thus, it is considered that τ_α and η provide essentially the same information on dynamics in supercooled liquids. However, this assumption does not hold true for the GCM at high shear rates. This result suggests that the shear modulus measured as G = η/τ_α is dependent on \(\dot{\gamma }\); G increases as \(\dot{\gamma }\) 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³⁰. The KA model is composed of large (A) and small (B) particles of equal masses, m_A = m_B = m. The particles interact via the LJ potential,

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 \(\tau ={(m{\sigma }^{2}/\epsilon )}^{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,50,51. We study at various temperatures ranging from T = 0.8 to 0.45 close to T_c. The shear rate \(\dot{\gamma }\) is controlled over a wide range of \(\dot{\gamma }=1{0}^{-6}\) to 10⁻¹. Note that \(\dot{\gamma } \sim 1{0}^{-1}\) is high enough that τ_α reaches the timescale of vibrations (Einstein period)⁴⁸.

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 \(\dot{\gamma }\) and identify the onset shear rate \({\dot{\gamma }}_{c}\). We present results for the larger particles (A) below, but similar results were obtained for the smaller particles (B). Figure 9 presents τ_α versus \(\dot{\gamma }\) in (a), and the data on \({\dot{\gamma }}_{c}\) and \({\tau }_{\alpha 0}^{-1}\) in (b) and (c). Figure 10 then plots \({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\) as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\), where we exclude data at \(\dot{\gamma }=1{0}^{-1}\) at which τ_α reaches the timescale of vibrations. We also measure the viscosity η and present data on \(\eta (T,\dot{\gamma })\) and \(\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) in Fig. 11. Furthermore, we compare the relaxation time and the viscosity in Fig. 12. Figures 9, 10, 11, and 12 for the KA model are counterparts of Figs. 4, 5, 7, and 8 for the GCM, respectively.

a \({\tau }_{\alpha }(T,\dot{\gamma })\) is plotted as a function of \(\dot{\gamma }\). Symbols of different colors represent values at different T; from bottom to top, T = 0.8 (yellow), 0.7 (purple), 0.6 (green), 0.55 (cyan), 0.5 (orange), 0.48 (blue), 0.46 (red), and 0.45 (black), all of which are above T_c ≃ 0.435. Closed symbols indicate the equilibrium values τ_α0(T). For T = 0.45, black line indicates \({\tau }_{\alpha }={\tau }_{\alpha 0}\propto {\dot{\gamma }}^{0}\), whereas blue line indicates \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.71 [Eq. (7)]. b \({\dot{\gamma }}_{c}\) and \({\tau }_{\alpha 0}^{-1}\) are plotted as a function of (T − T_c)/T_c. The lines present the power-law scalings, \({\tau }_{\alpha 0}^{-1}\propto {(T-{T}_{c})}^{\gamma }\) with γ ≃ 2.4 [Eq. (1)] and \({\dot{\gamma }}_{c}\propto {(T-{T}_{c})}^{\gamma \delta }\) with γδ ≃ 3.4 [Eq. (21)]. c \({\dot{\gamma }}_{c}\) is plotted against τ_α0. The blue line presents the scaling relation of \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\) with δ ≃ 1.42 [Eq. (6)].

Full size image

\({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\) is plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Symbols of different colors represent values at different T; T = 0.8 (yellow), 0.7 (purple), 0.6 (green), 0.55 (cyan), 0.5 (orange), 0.48 (blue), 0.46 (red), and 0.45 (black). Black and blue lines indicate τ_α/τ_α0 = 1 and \({\tau }_{\alpha }/{\tau }_{\alpha 0}\propto {(\dot{\gamma }/{\dot{\gamma }}_{c})}^{-\nu }\) with ν ≃ 0.71, respectively. The vertical dotted line indicates the location of the onset shear rate \({\dot{\gamma }}_{c}\).

Full size image

a \(\eta (T,\dot{\gamma })\) is plotted as a function of \(\dot{\gamma }\), and b \(\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) is plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Symbols of different colors represent values at different T; T = 0.8 (yellow), 0.7 (purple), 0.6 (green), 0.55 (cyan), 0.5 (orange), 0.48 (blue), 0.46 (red), and 0.45 (black). In a, closed symbols indicate the equilibrium values η₀(T). Black line indicates \(\eta ={\eta }_{0}\propto {\dot{\gamma }}^{0}\), whereas blue line indicates \(\eta \propto {\dot{\gamma }}^{-\nu }\) with ν ≃ 0.71. In b, the vertical dotted line indicates the location of \({\dot{\gamma }}_{c}\).

Full size image

a \({\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T),\,\eta (T,\dot{\gamma })/{\eta }_{0}(T)\) and b \(\left[{\tau }_{\alpha }(T,\dot{\gamma })/{\tau }_{\alpha 0}(T)\right]/{\dot{\gamma }}^{-\nu },\,\left[\eta (T,\dot{\gamma })/{\eta }_{0}(T)\right]/{\dot{\gamma }}^{-\nu }\) are plotted as a function of \(\dot{\gamma }/{\dot{\gamma }}_{c}\). Circles and squares represent values of τ_α and η, respectively. The temperature is T = 0.45. Black line indicates τ_α/τ_α0, η/η₀ = 1, wheares blue line indicates \({\tau }_{\alpha }/{\tau }_{\alpha 0},\,\eta /{\eta }_{0}\propto {(\dot{\gamma }/{\dot{\gamma }}_{c})}^{-\nu }\) with ν ≃ 0.71. The vertical dotted line indicates the location of \({\dot{\gamma }}_{c}\).

Full size image

In Figs. 9 and 10, we can see that the relaxation time follows a power-law scaling of the form \({\tau }_{\alpha }\propto {\dot{\gamma }}^{-\nu }\) where ν ≃ 0.71, and the onset shear rate \({\dot{\gamma }}_{c}\) exhibits the scaling behavior \({\dot{\gamma }}_{c}\propto {\tau }_{\alpha 0}^{-\delta }\ll {\tau }_{\alpha 0}^{-1}\) with δ ≃ 1.42. The values of ν and δ are obtained using ν = γ/(γ + 1) and δ = (γ + 1)/γ with γ = 2.4 as shown in Table 1.

In addition, Figs. 11 and 12 show that η is proportional to τ_α as η ∝ τ_α, following the same scaling law as that of τ_α;

Differently from the case of the GCM, the linear relation η ∝ τ_α is kept even at the high shear rates of \(\dot{\gamma }/{\dot{\gamma }}_{c}\gtrsim 1{0}^{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 η/η₀ versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at T = 0.45 in (a), and \({\dot{\gamma }}_{c}\) versus τ_α0 in (b).

System description of SS model

The SS model is a binary mixture composed of large (L) and small (S) particles in three spatial dimensions³¹. The particles interact via the inverse power-law potential,

where α and β denote L or S, and σ_αβ = (σ_α + σ_β)/2 with σ_α being the diameter of particle α. The mass and size ratios are m_L/m_S = 2 and σ_L/σ_S = 1.2, respectively. σ_S, ϵ/k_B, and \(\tau ={({m}_{S}{\sigma }_{S}^{2}/\epsilon )}^{1/2}\) are employed as units of length, temperature, and time, respectively. The number density is set to be ρ = N/V = (N_L + N_S)/V = 0.8, and the compositions of the two species are the same as N_L/N = N_S/N = 0.5. At ρ = 0.8, the MCT power-law fitting for τ_α0 estimates T_c ≃ 0.267^52,53. Note that refs. ^52,53 have studied the case at ρ = 0.742 and estimated T_c = 0.198. This value is converted to T_c ≃ 0.267 at ρ = 0.8 since one dimensionless coupling constant Γ = ρT^−1/4 determines the states of the SS model^5,31.

Reference ¹⁸ has studied the SS model under shear flow at T = 0.306, 0.285, 0.275, and 0.267, and reported the viscosity \(\eta (T,\dot{\gamma })\) as a function of \(\dot{\gamma }\). Figure 1 in the main text presents data on η/η₀ versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at T = 0.275 in (a), and \({\dot{\gamma }}_{c}\) versus η₀ at T≥0.275 > T_c in (b).

System description of 2DSS model

The 2DSS model is a binary mixture composed of large (L) and small (S) particles in two spatial dimensions^5,7. The particles interact via the inverse power-law potential as described in Eq. (27). The mass and size ratios are m_L/m_S = 2 and σ_L/σ_S = 1.4, respectively. σ_S, ϵ/k_B, and \(\tau ={({m}_{S}{\sigma }_{S}^{2}/\epsilon )}^{1/2}\) are employed as units of length, temperature, and time, respectively. The number density is set to be ρ = N/V = (N_L + N_S)/V = 0.8, and the compositions of the two species are the same as N_L/N = N_S/N = 0.5. We estimate T_c ≃ 0.534 by the MCT power-law fitting on data η₀ versus T reported in ref. ¹⁹.

Reference ¹⁹ has studied the 2DSS model under shear flow at T = 1.43, 0.85, 0.665, 0.577, and 0.526, and reported the viscosity \(\eta (T,\dot{\gamma })\) as a function of \(\dot{\gamma }\). Figure 1 in the main text presents data on η/η₀ versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at T = 0.577 in (a), and \({\dot{\gamma }}_{c}\) versus η₀ at T≥0.577 > T_c in (b).

System description of the BKS model

The BKS model is often used for amorphous and supercooled silica (SiO₂)³². Si and O ions interact via the potential,

where α and β denote Si or O. The values of the partial charges q_α and the constants A_αβ, B_αβ, and C_αβ are found in refs. ^32,54. The units of length and time are set to be 2.84Å and 1.98 × 10⁻¹³s, respectively. The temperature is measured in units of 6973.9K. The mass density is fixed at 2.37 g/cm³, which corresponds to the number density ρ = N/V = (N_Si + N_O)/V = 1.632. The dynamical transition temperature was estimated as T_c ≃ 3330K = 0.4775^55,56.

Reference ¹⁹ has studied the BKS model under shear flow at T = 0.614, 0.511, 0.47, 0.429, and 0.39, and reported the viscosity \(\eta (T,\dot{\gamma })\) as a function of \(\dot{\gamma }\). Figure 1 in the main text presents data on η/η₀ versus \(\dot{\gamma }/{\dot{\gamma }}_{c}\) at T = 0.511 in (a), and \({\dot{\gamma }}_{c}\) versus η₀ at T ≥ 0.511 > T_c in (b).