Local Density Fluctuation Governs the Divergence of Viscosity underlying Elastic and Hydrodynamic Anomalies in a 2D Glass-Forming Liquid

TL;DR

This paper demonstrates that in a 2D glass-forming liquid, local density fluctuations govern the divergence of viscosity, linking microscopic relaxation to macroscopic flow properties, with long-wavelength fluctuations playing a crucial role.

Contribution

It reveals that relative particle displacement correlates with viscosity growth and highlights the impact of elastic and hydrodynamic fluctuations on long-wavelength dynamics in 2D glasses.

Findings

Relaxation time scales with viscosity.

Long-wavelength fluctuations cause logarithmic divergence in diffusivity.

Hydrodynamic fluctuations significantly influence 2D glass dynamics.

Abstract

If a liquid is cooled rapidly to form a glass, its structural relaxation becomes retarded, producing a drastic increase in viscosity. In two dimensions, strong long-wavelength fluctuations persist, even at low temperature, making it difficult to evaluate the microscopic structural relaxation time. This Letter shows that, in a 2D glass-forming liquid, relative displacement between neighbor particles yields a relaxation time that grows in proportion to the viscosity. In addition to thermal elastic vibrations, hydrodynamic fluctuations are found to affect the long-wavelength dynamics, yielding a logarithmically diverging diffusivity in the long-time limit.

Equations

$$\textstyle F_{s}^{\rm R}(k,t)=(1/N)\left\langle\sum_{j=1}^{N}\exp\{i\mbox{\boldmath$k$}\cdot\Delta\mbox{\boldmath$r$}_{j}^{\rm rel}(t)\}\right\rangle.$$

$$G(t)=\frac{V}{k_{\mathrm{B}}T}⟨{\sigma }_{xy}(t){\sigma }_{xy}(0)⟩,$$

$$D^{\prime }(t)={\int }_0^{t}Z(t)dt$$

$$D\sim \frac{k_{\mathrm{B}}T}{8\pi \eta }\ln \left(\frac{L}{{\sigma }_{11}}\right).$$

Full text

In many two-dimensional ordering phenomena, fluctuations at long wavelengths are so strong that perfect order is destroyed. For example, the transition between a liquid and a crystalline solid is continuous or nearly continuous Halperin and Nelson (1978); Zahn et al. (1999); Han et al. (2008); Shiba et al. (2009); Kapfer and Krauth (2015); Thorneywork et al. (2017). Recently, large-scale molecular dynamics (MD) simulations Shiba et al. (2016); Flenner and Szamel (2015) and colloidal experiments Illing et al. (2017); Vivek et al. (2017) have revealed that such long-wavelength fluctuations also exist in two-dimensional (2D) liquids that are rapidly cooled toward the glass transition. Although retaining a random amorphous structure, elastic vibrations appear as the rigidity emerges with the decrease in temperature. The excess of low-frequency phonons in two dimensions Shiba et al. (2016); Mizuno et al. (2017) leads to an elastic anomaly, where the amplitude of thermal vibrations diverge at long wavelengths. Even in the presence of these long-wavelength fluctuations, the microscopic structural relaxations in 2D and 3D supercooled liquids appear to be similar, once the effect of these fluctuations has been eliminated by introducing quantities that characterize the local switching between neighbor particles Shiba et al. (2016); Vivek et al. (2017); Illing et al. (2017).

Albeit a similarity of structural relaxation modality between 2D and 3D glass-forming liquids, it does not mean that the transport properties, a key to the nature of glass transitions, are similar in between. One problem lies in the relationship between the structural relaxation time and the viscosity. The glass transition is marked by a drastic increase in macroscopic viscosity with decreasing temperature, which is intimately related to the divergence of the microscopic structural relaxation time. As such, theoretical and computational studies have focused on the dynamical mechanism of growth in the microscopic structural relaxation time, most typically the $\alpha$-relaxation time $\tau_{\alpha}$ Stillinger (1995); Ediger et al. (1996); Hurley and Harrowell (1996); W. Götze (2009), defined as the decay time of the relaxation function for density fluctuations, i.e., the intermediate scattering function. However, in two dimensions, the strong density fluctuation diverges at long wavelengths and suppresses its intermediate plateau. Thus, $\tau_{\alpha}$ no longer represents the microscopic structural relaxation time. The first problem arises on how to define the relaxation time that represents viscous slowdown of the dynamics.

Furthermore, a more intriguing problem lies in the potential role of hydrodynamic effects on the transport properties. For a liquid in two dimensions, a slow $t^{-1}$ decay (the so-called long-time tail) of the velocity and stress autocorrelation functions leads to a hydrodynamic anomaly that is marked by a logarithmic divergence of transport coefficients such as diffusivity and viscosity Ernst et al. (1970); Alder and Wainwright (1970); van der Hoef and Frenkel (1991); Hansen and McDonald (2006). How it alters the transport properties of deeply cooled liquids at the macroscopic level is a fully open issue, and the glass transition may be influenced by a mechanism that is different from the freezing of microscopic structural relaxation due to the cage effect. In fact, it is difficult to distinguish the characteristic time scales of the long-time tail and the microscopic structural relaxation, both of which become significantly large upon supercooling. Such long-wavelength fluctuations derived from hydrodynamics may even possibly prohibit the 2D glass transition. Therefore, in addition to the effect of elastic fluctuations that inhibit the existence of 2D crystals, it is important to reveal how macroscopic hydrodynamic fluctuations can influence the microscopic structural relaxation in 2D glass-forming liquids.

In this Letter, in order to address these issues, we perform large-scale MD simulations of a 2D glass-forming liquid to examine how the growth of various relaxation times is related to the divergence of macroscopic viscosity in the presence of long-wavelength fluctuations. The simulations are performed with a particular focus on how the dynamics depend on the system size, so that the anomalous enhancement of elastic and hydrodynamic responses can be characterized. Simulations are performed for a 2D variant of Kob–Andersen binary Lennard–Jones mixtures Kob and Andersen (1995), in which the composition is 65:35 Sengupta et al. (2013); Flenner and Szamel (2015, 2016, 2019). The mixture is annealed for a sufficiently long time (maximum of $4\times 10^{9}$ simulation steps) after rapid cooling to target temperatures in the range $0.4\leq T\leq 1.0$ in the presence of Langevin heat baths. The production runs are then performed as Newtonian ($NVE$) dynamics simulations to prevent the damping of long-wavelength fluctuations. The data presented in the remainder of this Letter are averaged over four or eight independent simulations sup .

First, we revisit three relaxation functions that have been considered in recent studies Flenner and Szamel (2015); Shiba et al. (2016); Vivek et al. (2017); Illing et al. (2017); Flenner and Szamel (2019). The first is the standard self-intermediate scattering function (SISF) F_{s}(k,t)=(1/N)\langle\sum_{j=1}^{N}\exp\{i\mbox{\boldmathk}\cdot[\mbox{\boldmathr}_{j}(t+t_{0})-\mbox{\boldmathr}_{j}(t_{0})]\}\rangle, with the wavevector set to |\mbox{\boldmathk}|=2\pi/\sigma_{11} so that its decay represents particle movement over a distance of the particle diameter $\sigma_{11}$. However, in two dimensions, long-wavelength elastic vibrations persist, and these enhance the mean-squared thermal displacement $\delta^{2}$ Shiba et al. (2016); Zhang and Cheng (2019). Figure  1 shows the SISF for a fixed temperature $T=0.4$, where the size-dependent behavior is in agreement with previous studies Flenner and Szamel (2015); Shiba et al. (2016). The plateau heights represent the Debye–Waller factor $f_{p}\sim\exp(-k^{2}\delta^{2}/2)$ Pastore et al. (1988); sup , and tend toward zero as the system size increases. Therefore, the $\alpha$-relaxation time (the decay time of the SISF) is strongly influenced by long-wavelength fluctuations and cannot represent the microscopic structural relaxation time.

The other two functions are the “neighbor-relative” SISF (also known as the “cage-relative” SISF Illing et al. (2017); Vivek et al. (2017); Tong and Tanaka (2018); Flenner and Szamel (2019); Li et al. (2019)) and bond relaxation function Yamamoto and Onuki (1998a); Shiba et al. (2012); Kawasaki and Onuki (2013); Flenner and Szamel (2016). These are also plotted in Fig. 1. The latter two functions represent relaxations in the sense that the long-wavelength thermal vibrations are removed. The former is defined by

[TABLE]

We here introduce the neighbor-relative displacement \Delta\mbox{\boldmathr}_{i}^{\rm rel}(t)=(1/N_{\rm n.n.})\sum_{j\in n.n.}[\Delta\mbox{\boldmathr}_{i}(t)-\Delta\mbox{\boldmathr}_{j}(t)] (also known as the “cage-relative” displacement Mazoyer et al. (2009); Illing et al. (2017); Vivek et al. (2017)), where the summation is over $N_{\rm n.n.}$ initially neighboring pairs of particles, indicating the changes in the relative positions. A similar displacement was considered in previous studies on 2D melting for the same purpose of eliminating long-wavelength fluctuations Bedanov et al. (1985); Zahn and Maret (2000); Han et al. (2008). The bond relaxation function $F_{\rm B}(t)$, in contrast, does not involve displacements of the particles, but simply characterizes the proportion of initially neighboring pairs that have survived after a certain time sup . From the observation that neither function is strongly dependent on the system size, contrary to the standard SISF, the effect of long-wavelength fluctuations is marginal for these relaxation functions, as expected from their definitions.

Along with our aim to relate these microscopic relaxation functions with the macroscopic viscosity, Fig. 1 also shows the stress relaxation function (or “dynamic modulus” Flenner and Szamel (2019))

[TABLE]

where $\sigma_{xy}(t)$ is the off-diagonal stress tensor (the data are normalized with respect to the instantaneous shear modulus $G_{\infty}=G(0)$). $G(t)$ exhibits a stretched plateau modulus and no system size dependence at a low temperature of $T=0.4$, as shown in Fig. 1. In a recent paper, the plateau was found to become unclear, rendering it difficult to evaluate the plateau modulus, for higher temperatures at the onset of slow dynamics, $T\geq 0.7$ Flenner and Szamel (2019).

Next, we define relaxation times and compare them with the transport coefficients —the $\alpha$-relaxation time $\tau_{\alpha}$, neighbor-relative relaxation time $\tau_{\rm R}$, and bond relaxation time $\tau_{\rm B}$ can be defined as the decay times of the standard SISF, neighbor-relative SISF, and bond relaxation functions sup . For this purpose, we refer to the Stokes–Einstein (SE) relation $D\eta/T={\rm const.}$ between the diffusivity $D$ and the viscosity $\eta$. This relation holds in normal liquids at high temperatures, but is violated in the deeply supercooled regime Tarjus and Kivelson (1995); Liu et al. (2006); Shi et al. (2013); Sengupta et al. (2013); Kawasaki and Onuki (2013). In simulation studies of glass-forming liquids, because the microscopic structural relaxation time is expected to grow proportionally with the viscosity, the left-hand quantity $D\eta/T$ (the so-called SE ratio) is often replaced by the product of the diffusivity and the relaxation time $D\tau_{\alpha}$, and its $T$-dependence is usually examined. This assumption may break down in two dimensions because $\tau_{\alpha}$ is robustly suppressed by the long-wavelength fluctuations. Thus, we calculate the temperature-dependence of the generalized SE ratio $D\tau$ to examine how the three relaxation times change with respect to the diffusivity $D$ as the temperature decreases. We further compare $D\tau$ with the original SE ratio $D\eta/T$ This is done by explicitly calculating the shear viscosity $\eta$ via the Green–Kubo formula $\eta=\int G(t)\ dt$ Sengupta et al. (2013); Flenner and Szamel (2019), for which we require an error estimate because of the slow convergence of this integral sup .

In Fig.2 (a)–(c), the original SE ratio $D\eta/T$ is shown for system sizes of $N=1000$, 16000, and 256000 as functions of the inverse temperature; the generalized SE ratios $D\tau_{\alpha}$, $D\tau_{R}$, and $D\tau_{B}$ are also shown. Although the diffusivity $D$ has a logarithmic dependence on the system size $N$, the generalized SE ratios can be meaningfully compared across different values of $N$ because of their similar temperature dependence sup . Firstly, the standard SE ratio increases as the temperature decreases, as in 3D systems. However, $D\tau_{\alpha}$ exhibits system-size dependence and nonmonotonic temperature dependence and is clearly decorrelated from $D\eta/T$, which is consistent with the results of a previous study Sengupta et al. (2013). In contrast, $D\tau_{R}$ collapses to the standard SE ratio for all system sizes, showing that the neighbor-relative relaxation time $\tau_{R}$ grows in proportion to the viscosity, satisfying $\tau_{R}\sim\eta/T$ (see also Fig. 2(d)). This relation provides an alternative to $\tau_{\alpha}\sim\eta/T$ Tarjus and Kivelson (1995); Yamamoto and Onuki (1998b); Shi et al. (2013); Kawasaki and Kim (2017), and thus $\tau_{\rm R}$ clearly takes the role of the microscopic structural relaxation time. We also find that the generalized SE ratio for bond relaxation is preserved, i.e., $D\propto\tau_{B}^{-1}$, indicating that the bond relaxation function is a descriptor for the 2D diffusive motion in a similar manner to other 3D supercooled liquids Kawasaki and Onuki (2013); Kawasaki and Kim (2017).

Thus far, we have seen that the diffusivity and viscosity are linked to time scales associated with local particle motion that is irrelevant to the long-wavelength fluctuations. However, this is not the end of our discussion, and we further investigate the hydrodynamic effects on the diffusivity by thoroughly examining the dependence on the system size. As shown in Fig. 3(a) for $T=0.4$, the mean-squared displacements (MSDs) \langle|\Delta\mbox{\boldmathr}(t)|^{2}\rangle=(1/N)\langle\sum_{i=1}^{N}|\mbox{\boldmathr}_{i}(t+t_{0})-\mbox{\boldmathr}_{i}(t_{0})|^{2}\rangle exhibit linear growth and remain dependent on the system size in the long-time limit. For temperatures $0.4\leq T\leq 1.0$, we further estimate the diffusivity $D$ by fitting D=\langle|\Delta\mbox{\boldmathr}(t)|^{2}\rangle/4t in the long-time region 10\leq\langle|\Delta\mbox{\boldmathr}(t)|^{2}\rangle\leq 20 for different system sizes. The result is shown as a function of the box length $L$ in Fig. 3(b). The diffusivity $D$ grows logarithmically with system size, even at low temperatures. This size dependence cannot arise from a simple superposition of elastic vibrational fluctuations (the so-called “Mermin–Wagner fluctuations” Illing et al. (2017)), but should be attributed to a different origin.

As a statistical measure for the motion of individual particles, the velocity autocorrelation function (VACF) Z(t)=(1/d)\langle\mbox{\boldmathv}(t)\cdot\mbox{\boldmathv}(0)\rangle Hansen and McDonald (2006) contains information regarding the delay in the viscoelastic responses of liquids. Importantly, VACF is related to the diffusivity via the Green–Kubo formula $D=\int_{0}^{\infty}Z(t)\ dt$, and should provide clues as to the system-size dependence. However, calculating the full resolution of VACF for a glass-forming liquid is a difficult task, because random motions of the caged particles blur the slow process of diffusion. In a high-density 3D liquid at a moderately high temperature, the VACF exhibits a negative correlation in the intermediate time regions Alder et al. (1970); Williams et al. (2006); Peng et al. (2016), which can be attributed to the velocity reversal caused by elastic vibrations of tagged particles Williams et al. (2006). Hence, for the present dense 2D liquid, we focus on high temperatures to demonstrate the crossover from transient elastic response to long-time hydrodynamic decay over a full time range. Figure 4(a) shows the VACF for different system sizes $N$ at $T=1.0$. Negative correlations exist for all $N$, indicating backward motion originating from elastic vibrations (see the plot for $N=256000$ in the inset). However, VACF does not simply decay from negative values to zero, but becomes positive over a longer time range for sufficiently large system sizes, such as $N\geq 256000$. By examining a much larger system size ($N=4096000$), the long-time limiting behavior is found to be consistent with the hydrodynamic $t^{-1}$ tail that appears in normal 2D liquids Ernst et al. (1970); Alder and Wainwright (1970); van der Hoef and Frenkel (1991); Hansen and McDonald (2006); Isobe (2008); Choi et al. (2017). Because the kinematic viscosity is large, i.e., $\nu=\eta/(nm)\gg D$, the analytical expression can be simplified to $Z(t)=(k_{\rm B}T/8\pi\eta)t^{-1}$. Both the magnitude and power-law exponent of VACF coincide with this expression. Therefore, the system size dependence of diffusivity $D$ is attributed to the purely hydrodynamic origin in the long-time limit where the transient elastic response vanishes.

Notably, the VACF itself exhibits system-size dependence at times before the $t^{-1}$ power-law tail. For each $N$, the VACF exhibits a systematic decrease before becoming uncorrelated. In Fig. 4(b), we show the size dependence of the finite-time diffusivity

[TABLE]

in which the long-time limit yields the long-time diffusivity $D=\lim_{t\rightarrow\infty}D^{\prime}(t)$. This finite-time diffusivity $D^{\prime}(t)$ is in good agreement with the diffusivity $D$ evaluated from MSD in the long-time limit. At the same time, it exhibits size dependence in an earlier time range, before converging to the long-time diffusivity $D$. Therefore, quite reasonably, the VACF $Z(t)$ itself is affected by the hidden hydrodynamic long-time tail, which has been difficult to find in simulations of 2D glass-forming liquids Perera and Harrowell (1998); Liu et al. (2006); Sengupta et al. (2013); Flenner and Szamel (2019).

The above results suggest that the origin of the logarithmic divergence of diffusivity is the $t^{-1}$ long-time tail and that it should exist even at low temperature. The long-time tail $Z(t)=(k_{\rm B}T/8\pi\eta)t^{-1}$ no longer involves the diffusivity $D$ in its expression and is free from the self-consistency problem leading to a faster than $t^{-1}$ decay van der Hoef and Frenkel (1991). Therefore, the diffusivity simply diverges as

[TABLE]

Hydrodynamic theories also predicts the $t^{-1}$ tail in the shear stress autocorrelation function, which may cause logarithmic divergence of shear viscosity Ernst et al. (1970). The shear stress autocorrelation functions is calculated for the same temperature but still such a power-law tail is not clear with the system size $N=4096000$ sup .

It is worth noting that the analytic expression Eq. (4) brings about increase of SE ratio $D\eta/T\sim\ln L$ over all the temperature ranges. This ratio increases by an order of magnitude if the system length $L$ increases by five orders of magnitude. While a common expectation seems to be that the effect of long-time tail is marginal in 2D glass transitions Perera and Harrowell (1998), it is clear that the long-time tail also causes violation of SE relation, indicating that it may influence the 2D glass transition.

Finally, we address the relevance of our results to recent studies. A recent simulation study on the same 2D Kob–Andersen liquid reported that the mean-square of the neighbor-relative displacement (i.e., cage-relative MSD) asymptotically approaches the normal MSD in the long-time limit Flenner and Szamel (2019). Together with our observation of the finite size effects in the MSD, their result implies that the neighbor-relative displacement is also system-size dependent, finally approaching the linear behavior of the usual MSD and giving rise to the same diffusivity. However, the neighbor-relative relaxation time $\tau_{\rm R}$ is short enough such that it remains virtually unaffected by the long-wavelength fluctuations. Therefore, our present results fully justify the usage of the equivalent cage-relative SISF in recent colloidal experiments Illing et al. (2017); Vivek et al. (2017).

We note that our results rely on the use of the $NVE$ ensemble to conserve the total momentum. Both the thermal vibrations and the hydrodynamic fluctuations at long wavelengths are suppressed by using the Brownian dynamics Flenner and Szamel (2015) or specialized Monte Carlo algorithms Berthier et al. (2019). It has long been assumed that glassy dynamics are unaffected by the choice of ensembles, as in 3D liquids Gleim et al. (1998), but this is not the case in two dimensions for quantities that are affected by long-wavelength elastic and hydrodynamic fluctuations, including the standard SISF and the MSD.

In conclusion, the dynamics of a 2D glass-forming liquid are covered by hydrodynamic power-law correlations that lead to the logarithmic divergence of diffusivity, in addition to the recently revealed long-wavelength elastic fluctuations arising from the emerging rigidity of the liquid. Both the elastic and hydrodynamic fluctuations persists at long-wavelengths to produce a concerted effects on the transport properties. Moreover, it is found that the relaxation time defined from the relative displacement between neighbors ($\tau_{\rm R}$) grows in proportion to $\eta/T$ as the temperature decreases, implying that local density fluctuations govern the drastic increase in viscosity. The combined elastic and hydrodynamic anomalies are expected to be relevant to both the existence of the 2D glass transition Berthier et al. (2019) and the dynamical features of 2D crystal melting, although further clarification is required in forthcoming studies.

Acknowledgements

We thank Hajime Tanaka, Kunimasa Miyazaki, and Patrick A. Bonnaud for fruitful discussions. This work was supported by JSPS KAKENHI Grant Numbers JP18K13513 (H. S.), JP19K03767 (T. K.), and JP18H01188 (K. K.). H. S. was also financially supported by Building of Consortia for the Development of Human Resources in Science and Technology, Ministry of Education, Culture, Sports, and Technology (MEXT), Japan. The numerical calculations were performed on the Cray XC50-LC at the Institute for Materials Research, Tohoku University, Japan, the NEC LX406Rh-2 at the Research Center of Computational Science, Okazaki Research Facilities, National Institutes of Natural Sciences, Japan, and the SGI ICE XA and HPE SGI 8600 at the Institute for Solid State Physics, University of Tokyo, Japan.