Microscopic origin of liquid viscosity from unstable localized modes

blueIntroduction\colorblack – Viscosity is a fundamental transport coefficient in fluid dynamics [1] that reflects the ability of transferring momentum in the direction transverse to the motion of the liquid [2]. In the hydrodynamic regime, viscosity governs the flow of physical systems at very different energy scales, from honey [3] to electrons in ultra-clean metallic systems [4] and ultra-relativistic plasma of quarks and gluons [5].

From a macroscopic point of view, shear viscosity is defined by the linear response relation between shear stress and shear strain rate, that can be formalized using Green-Kubo formalism [6, 7]. This method allows for a direct estimate of the viscosity based on the time correlation function of shear stress, but it is limited to Newtonian fluids. Another established method to compute the viscosity is based on the Einstein-Stokes relation between the self-diffusion constant of a probe particle and the mobility of the fluid in which the latter is flowing [8, 9]. However, it is well-known that, as temperature approaches the glass transition, this relation breaks down hindering the applicability of this method.

Moreover, it is important to stress that both the Green-Kubo formula and the Einstein-Stokes relationship do not fully elucidate the microscopic physical mechanisms underlying momentum transport, and ultimately the viscosity itself. In particular, what is the direct link between microscopic dynamics and macroscopic viscosity in a fluid remains an open question. Using a hyperbole, the physical origin of liquid viscosity has not been fully disclosed yet.

In fact, this issue is part of a bigger problem that is the lack of an atomic-scale microscopic description of liquid dynamics and thermodynamics [10, 11], which is mostly due to the absence of a well-defined normal mode formalism for liquids, as achieved on the contrary for crystalline solids [12]. In this context of liquids, the existence of normal modes was already proposed by Zwanzig in 1967 [13], on the basis of Maxwell’s intuition [14] that liquids present solid characteristics at short time scales. These ideas were later formalized by Keyes and collaborators within the so-called instantaneous normal mode (INM) approach [15, 16], that is based on the diagonalization of the Hessian matrix in instantaneous liquid configurations.

Because of the absence of a well-defined equilibrium configuration in liquids, this procedure inevitably leads to the appearance of negative eigenvalues [17], or equivalently purely imaginary frequencies, that correspond to regions of the potential energy landscape with negative local curvature. Modes corresponding to positive eigenvalues are stable, while those with purely imaginary frequencies are called unstable INMs. The physical meaning of these unstable INMs remains matter of debate, contributing to the famous question of “what INMs are and are not” [18].

Partially answering this question, in a series of works [19, 20], Keyes proposed that unstable INMs might provide the fundamental blocks to achieve a microscopic derivation of the liquid self-diffusion constant, highlighting their physical nature as relaxational processes mediated by jumps over potential barriers. It was later clarified that only unstable delocalized INMs contribute to the self-diffusion constant [21, 22] (see also [23, 24, 25, 26, 27, 23, 28] for related works and different definitions of this “diffusive” subset of INMs), providing a first connection between microscopic excitations and macroscopic transport in liquids, but leaving the physical meaning of the localized unstable modes unclear.

In this Letter, we investigate whether a microscopic definition for another fundamental transport coefficient in liquids – shear viscosity– can be achieved using a normal mode approach. In doing so, we will reveal the so-far unknown physical significance of unstable localized modes and their essential role for the microscopic origin of liquid viscosity.

We anticipate that a connection between high-frequency macroscopic elastic moduli and INMs in liquids was already proposed by Stratt [29], alluding to a possible connection between “certain” INMs and liquid viscosity as well. In this Letter, we corroborate that such a relation indeed exists and that unstable localized INMs are the building blocks of liquid viscosity. Our results provide a microscopic definition of viscosity, together with a parameter-free prediction for its value as a function of temperature.

blueTheoretical model and methods\colorblack – According to the theoretical formalism presented in [30] based on nonaffine viscoelastic response [31], the shear viscosity derived from the imaginary part of the complex shear modulus is given by

where $\rho=N/V$ represents the number density, $N$ is the number of atoms and $V$ is the volume of the system. Additionally, $g(\omega)$ is the density of states (DOS); $m$ is the atomic mass; $\Gamma(\omega)$ is the affine force field correlator [31] (see SM, Appendix D, for more details); $\tilde{\upsilon}(0)$ is the zero-frequency limit of the spectral density, $\tilde{\upsilon}(\omega)=\int_{0}^{\infty}\upsilon(t)e^{-i\omega t}dt$, with $\upsilon(t)$ the memory kernel arising from the coupling between tagged particles and thermal environment [32, 33]. All our following results and conclusions will rely on the validity of Eq. (1).

To derive the memory kernel, we use the projection operator formalism [34, 35, 36],

where $\mathbf{F}$ is the true force acting on the tagged particle. The random force $\mathbf{F}_{r}=e^{-i(1-\mathcal{P})\mathcal{L}t}\mathbf{F}$, $-i\mathcal{L}$ is the Liouvillian operator corresponding to the unperturbed dynamics and $\mathcal{P}$ is the Mori projection operator along the direction of the velocity (see SM, Appendix F, for more details). Consequently, the zero-frequency limit of spectral density can be simply expressed as

Finally, we present two important remarks in applying Eq.(1) to liquid systems. First, the definition of the DOS in liquids is not straighforward. In our work, we will identify $g(\omega)$ with the INM density of states, allowing a direct connection to normal modes (that is not possible using the velocity auto-correlation function). Second, following this choice, it is necessary to define the integration range on $\omega$, that now takes values along both the real and imaginary axes. In the following, in order to find the correct range of integration, we will consider different choices and compare them directly with the simulation data. As already discussed, for the case of the self-diffusion constant [19], only unstable delocalized modes should be considered. It is therefore conceivable that a similar situation will hold for the shear viscosity.

In this study, two simulations systems are considered: (I) a $\mathrm{Cu_{50}Zr_{50}}$ metallic liquid with $T_{\mathrm{g}}=690$ K [37] and (II) a Kob-Andersen (KA) simulation model [38, 39]. Ten independent samples were obtained to reduce the fluctuation error of the simulations. More details about the models and simulations are provided in the Supplementary Material (SM). The results in the main text refer to the metallic liquid; analogous results for the KA model are presented in the SM and present the same features showing the universality of microscopic origin of viscosity.

We notice that, once the range of integration on $\omega$ in Eq.(1) is fixed, our procedure will not involve any free fitting parameter. Therefore, our results could be taken as a direct test of the validity of Eq.(1) that, to the best of our knowledge, has never been verified.

blueMicroscopic description of liquid dynamics\colorblack – We performed INM analysis of the $\mathrm{Cu_{50}Zr_{50}}$ metallic liquid for temperatures between $800$ K and $2200$ K. The corresponding INM DOS is shown in Fig. 1(a), following the common practice of representing imaginary frequencies along the negative axes [15]. The arrows indicate the directions along which $T$ increases. We notice that $g(\omega)$ is linear at low frequency (see SM, Fig. S3, for more details), in agreement with theoretical expectations (e.g., [19, 40, 41]) and experimental observations [42, 43].

Panels (b) and (c) of Fig. 1 display respectively the re-scaled affine force field versus frequency and the auto-correlation function of the random force versus time at different temperatures. From the memory kernel defined in Eq. (2), we introduce the friction coefficient

Using the Einstein’s relation $\gamma=k_{\mathrm{B}}T/{D}$, where $D$ is the long-range diffusion constant of the tagged particles, and comparing Eqs. (3) and (4), we get

We find that computing $\tilde{\upsilon}(0)$ in this way gives results compatible with those obtained using the projection operator technique, Eq.(3), as shown explicitly in Fig. 1(d). The data also follow the empirical Vogel–Fulcher–Tammann (VFT) equation [44, 45, 46], $\eta(T)=\eta_{0}\exp\left(\frac{B_{\mathrm{VFT}}}{T-T_{0}}\right)$, shown with a red line in the same panel.

As anticipated, INMs can be classified into stable and unstable modes depending on the sign of the corresponding eigenvalues. The latter can be further divided into localized and delocalized, depending whether the number of particles participating in the mode is intensive or extensive. Bembenek and Laird originally argued that only unstable delocalized modes are related to diffusion [21, 47], while unstable localized modes do not contribute to it. However, there are indications that unstable localized modes may play a fundamental role for other dynamical properties, such as structural relaxation [48, 49].

In order to quantify the degree of localization, we consider the participation ratio

where $\mathbf{e}_{i}^{\omega}$ is the contribution of particle $i$ to the normalized eigenvector corresponding to eigenfrequency $\omega$. $P_{\omega}$ takes values between $1/N$ and $1$, with the former corresponding to an extremely localized mode with only one particle participating, while the latter to a fully delocalized mode involving all particles (see SM, Fig. S11, for direct visualization of these two types of modes and Ref. [50] for details about $P_{\omega}$).

In order to identify the mobility edge separating localized and delocalized modes, we follow [21, 47] and perform a finite size scaling analysis and compute $P_{\omega}$ for different $N$. In Fig. 2(a), we show $P_{\omega}$ as a function of frequency for the unstable branch of INMs and different system sizes. From there, we can observe that $P_{\omega}$ becomes system dependent below a certain threshold frequency, that for the $\mathrm{CuZr}$ glass-forming system corresponds to $P_{\omega}\approx 0.37$. Intensive modes, with $P_{\omega}$ below the mobility edge, are localized, while extensive modes with $P_{\omega}\gtrapprox 0.37$ are delocalized.

The participation ratio as a function of frequency for different temperatures is shown in Fig. 2(b). By repeating the same finite size scaling analysis for different $T$, we found that the value $P_{\omega}\approx 0.37$, defining the mobility edge, is in first approximation independent of temperature (see SM, Fig. S10, for a direct proof of this statement). By then searching for the corresponding (imaginary) frequency (see horizontal dashed line in Fig. 2(b)), we have derived the mobility edge frequency that is plotted in Fig. 2(c) as a function of temperature. We find that the mobility edge frequency increases (in absolute value) by increasing temperature, consistent with the idea that unstable modes tend to be more localized towards the glass transition temperature [21].

Finally, in Fig. 2(d), we respectively show the fraction of stable $f_{\mathrm{s}}$, unstable $f_{\mathrm{u}}$, unstable localized $f_{\mathrm{u}}^{\mathrm{L}}$ and unstable delocalized $f_{\mathrm{u}}^{\mathrm{D}}$ INMs in the temperature range considered. We find that the fraction of stable modes slowly decreases by increasing temperature, while all the other fractions increases with temperature following a similar trend. Once the temperature reaches the glass transition, unstable delocalized modes disappear and all unstable modes become localized, consistent with previous results [22].

We have now independently derived all the physical quantities entering in Eq.(1) and analyzed in detail the different types of modes over which the integration therein could be in principle performed. Before moving to the computation of the liquid viscosity, we notice that, due to the limited size of the simulation box $L$, the integration in Eq.(1) will be always performed up to a lower cutoff frequency $|\omega_{\min}|=(2\pi v_{s})/L$, where the speed of sound $v_{s}$ is given by $\sqrt{B/\bar{\rho}}$, with $B$ the bulk modulus and $\bar{\rho}$ the mass density (see SM, Appendix H, for more details).

blueLiquid viscosity\colorblack – Traditionally, the viscosity can be computed using the Green-Kubo formula [6, 7],

where $\sigma_{\alpha\beta}$ is the $\alpha,\beta$ component of the stress tensor (with $\alpha,\beta=x,y,x$). In Eq.(7), $\langle\cdots\rangle$ indicates ensemble average. Alternatively, the viscosity can be obtained from non-equilibrium molecular dynamics using the linear response relation $\eta=\sigma_{\alpha\beta}/\dot{\varepsilon}$, where $\dot{\varepsilon}$ is the shear strain rate. We will indicate this second route as shear MD.

In Fig. 3(a), we compare the viscosity calculated by different methods: black symbols corresponds to the Green-Kubo formula Eq.(7), red crosses to shear MD and blue symbols are the results obtained from Eq.(1), without free parameters, where the integration is performed over only the unstable localized INMs (ULINMs). For comparison, we also show the experimental data from Refs. [51]) (purple diamonds) and [52] (green crosses). At high temperature (above $T\approx 1000$ K, $1.5$ times the glass transition temperature), Green-Kubo formula and shear MD are in good agreement. However, when temperature drops, the shear viscosity becomes strain rate dependent and the linear response formula (hence, the shear MD method) is not anymore applicable, giving unreliable predictions for $\eta$.

Most importantly, we find that the viscosity computed via the Green-Kubo formula is in perfect agreement with the results of Eq.(1) using only ULINMs, and both data are well described by the VFT equation. This is the main result of this work. Same results are obtained for the Kob-Anderson liquid model and are presented in the SM, Fig. S1, confirming the excellent agreement found for the CuZr metallic system and the universality of our conclusions.

For completeness, in Fig. 3(b) we show a comparison between the viscosity obtained from Green-Kubo and the viscosity obtained from the normal mode formula Eq.(1) using different integration regimes or, in other words, different types of INMs as the microscopic building blocks for viscosity. We find that only the estimate using ULINMs is in agreement with the Green-Kubo data. All other choices give results that are almost an order of magnitude away from the correct value of $\eta$, especially at high temperature. This suggests that unstable localized modes are the microscopic origin of liquid viscosity, and that other excitations are unrelated to this macroscopic transport property.

Inspired by Adam-Gibbs theory [53], and motivated by recent results on the relation between liquid diffusion and excess entropy [54], we construct a theoretical model to quantitatively correlate ULINMs with macroscopic viscosity. Following the ideas in [22], we make a similar attempt to link ULINMs to configurational entropy, and then viscosity. First, in Fig. 4(a) we find that the configurational entropy $S_{\mathrm{conf}}$ displays a linear relationship with $\ln(f_{\mathrm{u}}^{\mathrm{L}})$. Here, the configurational entropy is estimated by subtracting the vibrational entropy from the total entropy obtained by thermodynamic integration, as detailed in [55]. Second, we observe that the viscosity follows a direct non-Arrhenius relation with the configurational entropy, $\eta\propto A\exp\left(B/\left(TS_{\text{conf }}\right)^{\lambda}\right)$, with $\lambda=2.10$. We notice the analogy with the connection between the $\alpha$ relaxation time and $S_{\text{conf }}$, usually rationalized using random first-order transition theory [53, 56, 57]. By combining these two observations, we conclude that $\eta\propto A\exp\left(B/\left(T\left[a\ln\left(f_{\mathrm{u}}^{\mathrm{L}}\right)+b\right]\right)^{\lambda}\right)$, as demonstrated using the simulation data in Fig. 4(c). This result provides a further link between viscosity and ULINMs and validate their intimate relation, providing a bridge between microscopic and macroscopic physics in liquids.

blueConclusions\colorblack – In this Letter, we have performed MD simulations combined with instantaneous normal mode analysis to disclose the microscopic origin of viscosity in liquids. Our results indicate that unstable localized INMs are the microscopic excitations responsible for this process. Moreover, our work provides a parameter-free formula to predict from normal modes the viscosity of liquids and its temperature dependence, defying the famous Landau argument that is “impossible to derive any general formulae giving a quantitative description of the properties of a liquid” [58].

More in general, our results contribute to answering the question of what unstable INMs are and are not. Unstable delocalized modes are the building blocks for liquid self-diffusion, while the complementary ULINMs are the microscopic facilitators underlying liquid viscosity.

Interestingly, this answer complements several previous analyses. In particular, in Ref.[49] it was shown that localized soft modes are the origin of irreversible structural reorganization. Here, we give a microscopic definition of these localized soft modes and we propose that they should be identified with ULINMs. This is also consistent with the results of [48] that showed that ULINMs act as facilitators of the liquid dynamics, located at the boundary regions between mobile and immobile clusters generated by dynamical heterogeneity. Moreover, our findings are aligned with the idea that, in glass-forming liquid, excitations are localized and relaxation is hierarchical [59]. Finally, viscosity being governed by localized modes is also reflected in the recent results of [60], where it was found that the change in the local barrier energy caused by the fundamental rearrangement or excitation of a few particles (hence, localized modes) determines the change in the activation energy and the dynamics of the supercooled liquids.

In summary, it was repeatedly recognized that “the viscosity of liquids is a subject which has hitherto been without any general theoretical basis” [61]. By identifying the microscopic carriers of momentum transport, our work suggests a possible microscopic and quantitative theory of liquid viscosity, complementing the previous approaches [1] and providing a new path towards a complete understanding of liquid dynamics.

blueAcknowledgments \colorblack – We thank A. Zaccone for collaboration at an early stage of this work and for valuable comments and suggestions. We would like to thank J. Douglas, M. Pica Ciamarra, J. Dyre, T. Keyes, J. Moon and A. Lemaitre for useful comments and discussion. This work was financially supported by the Strategic Priority Research Program (grants nos. XDB0620103 and XDB0510301) and the Youth Innovation Promotion Association of Chinese Academy of Sciences, and the National Natural Science Foundation of China (grant no. 12072344). The numerical calculations in this study were carried out on the ORISE Super-computer. M.B. acknowledges the support of the Shanghai Municipal Science and Technology Major Project (Grant No.2019SHZDZX01) and the sponsorship from the Yangyang Development Fund.

Supplementary Material

In this Supplementary Material (SM), we provide further details about the numerical and analytical computations presented in the main text. Moreover, we show further analysis on a Kob-Andersen model to corroborate our findings.