Long-wavelength fluctuations and anomalous dynamics in two-dimensional liquids

TL;DR

This paper investigates how long-wavelength fluctuations, caused by Mermin-Wagner effects, influence the dynamics of two-dimensional liquids, revealing a decoupling between diffusion and relaxation at high temperatures through combined numerical and experimental methods.

Contribution

It demonstrates that long-wavelength fluctuations affect high-temperature liquid dynamics, causing a decoupling between diffusion and structural relaxation not previously recognized.

Findings

Long-wavelength fluctuations induce a D ∝ τ^(-κ) decoupling in 2D liquids.

This decoupling occurs even without transient solid-like responses.

Fluctuations have negligible effects only at extremely high temperatures or low densities.

Abstract

Long-wavelength Mermin-Wagner fluctuations prevent the existence of translational long-range order, in two-dimensional systems at finite temperature. Their dynamical signature, which is the divergence of the vibrational amplitude with the system size, also affects disordered solids and washes out the transient solid-like response generally exhibited by liquids cooled below their melting temperature. Through a combined numerical and experimental investigation, here we show that long-wavelength fluctuations are also relevant at high temperature, where the liquid dynamics does not reveal a transient solid-like response. In this regime, they induce an unusual but ubiquitous decoupling between long-time diffusion coefficient $D$ and structural relaxation time $\tau$, where $D\propto \tau^{-\kappa}$, with $\kappa > 1$. Long-wavelength fluctuations have a negligible influence on the relaxation dynamics only at extremely high temperatures, in molecular liquids, or extremely low densities, in colloidal systems.

Full text

I Introduction

The dimensionality of a system strongly influences the equilibrium properties of its solid phase Chaikin and Lubensky (1995). According to the Mermin and Wagner theorem Mermin and Wagner (1966), indeed, systems with continuous symmetry and short-range interactions lack true long-range translational order at finite temperature, in d $\leq$ 2 dimensions. This occurs as in small spatial dimensions the elastic response is dominated by the Goldstone modes, elastic excitations that in the limit of long-wavelength (LW) have a vanishing energy and a diverging amplitude. A signature of these LW fluctuations is a system size dependent dynamics, which arises as the system size provides a cutoff for the maximum wavelength. This dependence occurs in both ordered and disordered solids, as LW fluctuations are insensitive to the local order Illing et al. (2017). This dynamical signature of the LW fluctuations also appears in supercooled liquids, which are liquids cooled below their melting temperature without crystallization occurring. In particular, LW fluctuations affect the transient solid-like response observed in the supercooled regime, which we recap in Fig. 1 via the investigation of the mean square displacement (MSD), $\Delta r^{2}(t)$, and of the self-intermediate scattering function (ISF), $F_{s}(q,t)$, of the two-dimensional (2d) modified Kob-Andersen (mKA) system, a prototypical model glass former (see Materials and Methods). Notice that in the supercooled regime the ISF and the MSD develop a plateau revealing a solid-like response in which particles vibrate in cages formed by their neighbors Debenedetti and Stillinger (2001). However, this is a transient response as at longer times the ISF relaxes and the MSD enters a diffusive behavior. Flenner and Szamel Flenner and Szamel (2015) have demonstrated that this glass relaxation dynamics depends on the system size, the signatures of a transiet solid-like response disappearing in the thermodynamic limit. This trend can be appreciated in Fig. 1, where we compare the dynamics for two different system sizes, but see Ref. Flenner and Szamel (2015) for a full account. Following works have then demonstrated that this size dependence results from the LW fluctuations Illing et al. (2017); Vivek et al. (2017) by showing that the glassy features of the relaxation dynamics are recovered when the effect of LW fluctuations is filtered out Shiba et al. (2016); Hayato et al. (2018); Vivek et al. (2017); Illing et al. (2017).

A transient solid-like response is only observed below the onset temperature, where the relaxation time exhibits a super-Arrhenius temperature dependence, in fragile systems (see SI Appendix Fig. S2), and dynamical heterogeneities (DHs) affect the relation between diffusion coefficient and relaxation time Debenedetti and Stillinger (2001). In the normal liquid regime that occurs at higher temperatures in molecular liquids, or at lower densities in colloidal systems, the MSD and the ISF do not exhibit plateaus possibly associated to a transient particle localization, and are system size independent. This is apparent in Fig 1. This suggests that the normal liquid regime is not affected by LW fluctuations. Is this true? And more generally, how far a system should be from the solid phase for the LW fluctuations to have a negligible influence on the relaxation dynamics? Here we show that, surprisingly, LW fluctuations affect the structural relaxation dynamics of 2d systems in their normal liquid regime. Specifically, they induce DHs, qualitatively distinct from that observed in the supercooled regime, and an unusual decoupling between the structural relaxation time $\tau$ and the long-time diffusion coefficient, $D$ (Materials and Methods). This decoupling has been previously observed, both in experiments Mishra and Ganapathy (2015) of colloidal systems as well as in numerical simulations Sengupta et al. (2013); Kim and Sung (2015), but its physical origin remained mysterious. Our results are based on the numerical investigation of the relaxation dynamics of two model glass-forming liquids and on the experimental study of a quasi-2d suspensions of ellipsoids Mishra and Ganapathy (2015). Numerically, we consider the three-dimensional (3d) Kob-Andersen binary mixture (KA) Kob and Andersen (1994) and its 2d variant (mKA) Brüning et al. (2009), as well as the harmonic model O’Hern et al. (2003) (Harmonic) in both 2d and 3d. Numerical details are in Materials and Methods. Details on the experimental systems are in Refs. Mishra et al. (2013); Mishra and Ganapathy (2015). Our results thus demonstrate that the Mermin-Wagner LW fluctuations are critical to rationalize the properties of 2d systems not only in the crystalline Mermin and Wagner (1966); Strandburg (1988), amorphous solid Illing et al. (2017) and supercooled phase Vivek et al. (2017), but also, surprisingly, in the high temperature normal liquid regime.

II Results

II.1 Heterogeneous dynamics in the 2d normal liquid regime

If the long-time relaxation dynamics of a liquid is not influenced by its energy landscape, as suggested by Fig. 1 and by analogous results for the 2d Harmonic model reported in SI Appendix Fig. S1, then at the relaxation time scale the displacements of the particles should be uncorrelated, and their van-Hove distribution function a Gaussian. This is observed in Fig. 2A, where we plot the dependence of the non-Gaussian parameter, a common measure of DHs Weeks et al. (2000) (see Materials and Methods), on the relaxation time, for 2d mKA and the 3d KA models. Indeed, above the onset temperature, where the relaxation time is that indicated by the dash vertical lines in Fig. 1 (and in the the following figures), $\alpha_{2}(\tau)\ll 1$, and assumes comparable values in 2d and in 3d.

However, the non-Gaussian parameter is mostly sensitive to deviations from the Gaussian behavior of the tails of the van-Hove distribution, which is populated by the particles with large displacements. A measure which conversely probes deviations arising from the particles with small displacements can be extracted from the long-time decay of the ISF. Indeed, if the van-Hove distribution is Gaussian at long time, then the ISF decays exponentially Hansen and McDonald (2005), while a stretched exponential relaxation $F_{s}(q,t)\propto\exp(-(t/\tau_{e})^{\beta})$ with $\beta<1$ reveals DHs, as commonly observed in the supercooled regime Sastry et al. (1998). Surprisingly, in 2d we observe $\beta<1$ also in the normal liquid regime, as in Fig. 2B. Consistently, the dynamical susceptibility at the relaxation time, $\chi_{4}(\tau)$, which is a proxy for the presence of correlated spatio-temporal motion (see Materials and Methods), grows in the 2d normal liquid regimes, as in Fig. 2C. The differences in the behavior of $\beta$ and of $\chi_{4}(\tau)$ in 2d and in 3d, which are apparent comparing the left and the right column of Fig. 2, signal the presence of DHs in the 2d normal liquid regime.

II.2 Decoupling between relaxation and diffusion

The heterogeneities observed in the normal liquid regime of 2d systems are qualitatively different from those observed in the supercooled regime, as they occur in the presence of small values of the non-Gaussian parameter. Hence, these hetereogeneities do not reflect the coexistence of particles with small and large displacements, a hallmark of the supercooled liquid dynamics Weeks et al. (2000). In supercooled liquids, this coexistence induces a breakdown of the inverse relationship between the diffusion coefficient and the structural relaxation time, leading to $D\propto\tau^{-\kappa}$ with $\kappa<1$, being $D$ and $\tau$ primarily affected by the particles with the large and small displacements, respectively. In 2d normal liquids, we therefore not expect to observe $\kappa<1$. Indeed, and for reasons that are currently mysterious, one finds $\kappa>1$ Sengupta et al. (2013); Kim and Sung (2015); Mishra and Ganapathy (2015). We illustrate the crossover from the normal liquid regime where $\kappa>1$, to the supercooled regime where $\kappa<1$, in Fig. 3 (black squares). The figure reports results from numerical simulations of the 2d mKA model and of 2d Harmonic discs, as well as results from experiments of quasi-2d colloidal ellipses. Reader is warned to take with care the precise location of the relaxation time at the onset temperature, dashed line, as well as the value of $\kappa$ in the supercooled regime. Indeed, as we have already seen in Fig. 2, the supercooled regime dynamics has finite-size effects. We will give more details on these size effects later on. Conversely, we stress that results in the normal liquid regime, which is our main focus, are robust.

We provide an insight on the physical origin of the observed heterogeneities, as well as of the $\kappa>1$ value, by investigating the dynamics using cage-relative (CR) measures (see Materials and Methods), where the displacement of a particle is evaluated with respect to average displacement of its neighbors. In Fig. 3 we compare the standard and the CR measures by plotting both $D$ vs. $\tau$ (black squares) and $D_{\rm CR}$ vs. $\tau_{\rm CR}$ (red circles). In the normal liquid regime we find $\kappa>1$ for the standard measures, and $\kappa\simeq 1$ for the CR measures. To better illustrate that there exists an extended range of $\tau$ where $\kappa>1$ for the standard measure, and $\kappa=1$ for the CR measures, we show in SI Appendix Fig. S3 the presence of a plateau region when the product $\tau_{\rm CR}D_{\rm CR}$ is plotted vs. $\tau_{\rm CR}$.

These results can be interpreted considering that CR measures have been suggested to filter out the effect of the LW fluctuations on the dynamics Vivek et al. (2017), in the supercooled regime, as they subtract correlated particle displacements. Our results suggest that the LW fluctuations are also responsible for the DHs in the normal liquid regime of 2d systems, and for the $\kappa>1$ value. In this respect, notice that this large $\kappa$ value indicates that the relaxation time is smaller than expected, given the value of diffusion coefficient, an effect possibly arising from the LW fluctuations. Indeed, the LW fluctuations have a large amplitude which, when larger than $\frac{2\pi}{q}$, promotes the relaxation of the system but not its diffusion. This interpretation is supported by the results of Fig. 3, from which we understand that the standard and the CR measures mainly differ in their estimation of the relaxation time, which is smaller for the standard measure (see SI Appendix, Fig. S4 for a direct comparison of the standard and CR MSD and ISF). This interpretation is also consistent with the investigation of the difference between the standard and the CR measures in 3d. Indeed Fig. 4 clearly shows that in 3d, where LW fluctuations play a minor role, there are no noticeable differences between the two measures, neither in the normal liquid nor in the supercooled regime. We finally notice that Fig. 3 reveals a small difference between the standard and the CR 2d diffusivities, that is not observed in 3d. This difference cannot originate from the LW fluctuations as it occurs at long-times, as in the diffusive regime there is no elastic response. Rather, they may originate from 2d hydrodynamic fluctuations Alder and Wainwright (1970); Hayato et al. (2019).

II.3 Long-wavelength fluctuations in the 2d normal liquid regime

The above results suggest that in 2d the dynamics of normal liquids is strongly influenced by the LW modes. This is a counter-intuitive speculation as the relaxation dynamics is expected to be weakly dependent on the features of the underlying energy landscape Sastry et al. (1998); Debenedetti and Stillinger (2001) above the onset temperature. To prove the relevance of LW fluctuations in the relaxation dynamics of 2d normal liquid, we therefore perform two additional investigations.

First, we consider the effect of the microscopic dynamics on the value of $\kappa$ comparing the two limiting cases of underdamped and overdamped dynamics, corresponding to a Langevin dynamics (see Materials and Methods) with Brownian time $\tau_{B}$ respectively much larger and much smaller than the inverse Debye frequency $1/\omega_{D}$. If the LW modes are responsible for the observed behavior, then the value of $\kappa$ must depend on the microscopic dynamics, and $\kappa$ closer to one should be obtained for the overdamped dynamics. Indeed, recall that in the solid phase the contribution of modes with frequency $\omega$ to the MSD is $r_{u}^{2}(t,\omega)=\frac{2k_{b}T}{m\omega^{2}}[1-\cos(\omega t)]$, in the underdamped limit, and $r_{o}^{2}(t,\omega)=\frac{2k_{b}T}{m\omega^{2}}[1-e^{-\frac{t}{2\tau_{B}}}]$ in the overdamped one Wang and Uhlenbec (1945). Accordingly, while in the overdamped limit the contributions of the different modes have the same time dependence, this is not so in the underdamped limit. Specifically, in the underdamped limit the modes contribute ballistically to the mean square displacement up to a time $\sim 1/\omega$, which implies that the underdamped dynamics is more strongly affected by the LW modes. Figure 5 compares the relation between $D$ and $\tau$ obtained in the underdamped and in the overdamped limits, for the 2d mKA model, with $\tau_{B}\omega_{D}\simeq 100$ and $\simeq 0.01$ respectively, having assumed $\omega_{D}$ to be of the order of the time at which the ballistic regime ends in the Newtonian dynamics (no damping). The underdamped results are analogous to those of the previously considered Newtonian dynamics reported in Fig. 3$A$, as expected. Conversely, the overdamped results strongly differ in that $\kappa\simeq 1$ is essentially recovered in the normal liquid regime. This agrees with our theoretical argument, and clarifies that the decoupling behavior with $\kappa>1$ is the combined effect of the presence of LW modes, and of the specific microscopic dynamics through which the system explores its phase space.

As a second check we explicitly evaluate the relevance of the LW modes to the overall particle displacement, as a function of time. To this end, we project the normalized displacement $\widehat{\Delta{\mathbf{r}}}(t)=\Delta\mathbf{r}(t)/|\Delta\mathbf{r}(t)|$ of the particles at time $t$ on the eigenvectors $u_{i}(\omega_{i})$ of the Hessian matrix of the nearest inherent configuration (see Materials and Methods) of the $t=0$ configuration: $\widehat{\Delta\mathbf{r}}(t)=\sum_{i}\beta_{i}(t)\mathbf{u}_{i}(\omega_{i})$. $\beta^{2}_{i}(t)$ is therefore the relative contribution of mode $i$ to the overall displacement at time $t$. To asses the relevance of the LW modes we compute the contribution of the 0.75% modes with the longest wavelength, $W_{0.75\%}(t)=\sum_{i}\beta_{i}^{2}(t)$. Figure 6 compares $W_{0.75\%}(t)$ for the underdamped and the overdamped dynamics, both in the liquid and in the supercooled regime. In both regimes $0.75\%$ of the longest wavelength modes contribute more than $30\%$ of the overall displacement. This is a clear indication that LW fluctuations are extremely relevant, also in the normal liquid regime. Consistently with our previous finding, we also find the contribution of the LW modes to be more relevant for the underdamped rather than for the overdamped dynamics. In addition we stress that, in the underdamped regime, for low frequencies $\beta_{i}(t)$ has a transient oscillatory behavior in line with the presence of transient solid-like modes, also in the liquid regime. We finally mention that we have verified that these results do not depend on the system size, in the normal liquid regime. Conversely, in the supercooled regime $W_{0.75\%}(t)$ slightly increases with the system size. This system size dependence is consistent with that expected in the presence of LW fluctuations.

II.4 Stokes-Einstein relation

We finally investigate whereas the breakdown of the inverse proportionality between $D$ and $\tau$ implies that of the Stokes-Einstein (SE) relation, $D\propto(\eta/T)^{-1}$, where $\eta$ is the viscosity of the system. This is possible as $\eta$ and $\tau$ are generally related. To verify this possibility, we measure the viscosity as the Green-Kubo integral Hansen and McDonald (2005) of the shear stress autocorrelation function (see Materials and Methods), for both the 3d KA and the 2d mKA model, and study its dependence on the relaxation time. In 3d, we find $\eta/T\propto\tau\propto\tau_{CR}$, as illustrated in Fig. 7$A$ and previously reported in different systems Sengupta et al. (2013); Kawasaki and Kim (2017). In the 2d normal liquid regime $\eta/T\propto\tau^{\zeta}$ with $\zeta=1.68$ and no system size dependence. Since in this regime $D\propto\tau^{-\kappa}$ with $\kappa\simeq\zeta$, as in Fig. 3$A$, the SE relation holds. In the supercooled regime we observe a system size dependence, which we highlight also considering data from Ref. Flenner and Szamel (2019) for systems with up to 4 million particles: the value of the scaling index $\zeta$ increases with the system size, and saturates when the system size is large enough. The saturation occurs as the time needed for the modes with the longest wavelength to develop, $\propto 1/L$, becomes larger than the relaxation time. Since the asymptotic value is $\zeta>1$, while in the supercooled regime $\kappa<1$ (see Fig. 3), the SE relation breaks down. Hence, both in 2d and 3d the SE relation only breaks down in the supercooled regime. In the 2d normal liquid regime the SE does not break down as $\eta$ is determined from the shear stress which, being related to the interparticle distances, is not affected by correlated particle displacements, as the CR measures. Indeed, we show in Fig. 7$B$ (open symbols) that $\eta/T$ is proportional to the CR relaxation time.

III Discussion

In conclusions, our results indicate that LW fluctuations affect the structural relaxation of 2d liquids in the normal liquid regime, where the relaxation dynamics does not suggest a transient solid response, not even when evaluated using CR measures (see SI Appendix, Fig. S4). In this regime the LW fluctuations induce DHs qualitatively different from that observed in the supercooled regime, which is not associated to the coexistence of particles with markedly different displacements. This result allows to rationalize an open issue in the literature Sengupta et al. (2013); Mishra and Ganapathy (2015); Kim and Sung (2015), namely the physical origin of the decoupling between relaxation and diffusion $D\propto\tau^{-\kappa}$ with $\kappa>1$. In the main manuscript, we have presented numerical data for the 2d mKA model and the 2d Harmonic model, and their 3d counterpart for comparison, as well as an experimental data for a 2d suspension of hard ellipses. We note, however, that we have also observed consistent results in a binary system with inverse-power-law potential Perera and Harrowell (1998) and in monodisperse systems of Penrose kites Zhao and Mason (2015); Zong et al. (2018); Li et al. (2019) (see SI Appendix, Fig. S5). Thus, our findings appear extremely robust as they do not depend on whether the interaction potential is finite or diverging at the origin, attractive or purely repulsive, isotropic or anisotropic.

It is natural to ask if, at high enough temperature or low enough density, the effect of LW fluctuations becomes negligible. The answer to this question is affirmative. Indeed, we do see in Fig. 3 that the difference between the standard and the CR measures, which is a proxy for the relevance of LW fluctuations, decreases as the relaxation time decreases. In the numerical models we have explicitly verified that the two measures coincide in this very high temperature limit. Interestingly, we have found that in this limit the system relaxes before the ballistic regime of the mean square displacement ends. This leads to $D\propto\tau^{-\kappa}$ and $\kappa=2$, in both 2d and 3d, as we discuss and verify in SI Appendix, Fig. S6. In colloidal systems particles perform independent Brownian motions in the low density limit, where LW fluctuations are therefore negligible. We expect the crossover density below which LW fluctuations are negligible to depend on the viscoelasticity of the solvent.

We conclude with two more remarks. First, it is established that CR measures remove the effect of LW fluctuations Vivek et al. (2017); Illing et al. (2017). Here we notice that CR measures filter out all correlated displacements between close particles, regardless of their physical origin. In particular, in the supercooled regime they suppress the effect of correlated particles displacements arising from DHs (see SI Appendix, Fig. S7 for the comparison of the four-point dynamical susceptibility between standard and CR measures). This has to be taken into account when using 2d systems to investigate the glass transition. We notice that it appears difficult to selectively suppress only the correlations arising from one of these two physical processes, as DHs are associated to the low frequency vibrational modes Widmer-Cooper et al. (2008). In this respect, perhaps one may consider that DHs in the supercooled regime are associated to localized modes, while LW fluctuations are signatures of extended modes.

Finally, we highlight that LW fluctuations are found in quasi-2d colloidal experiments of both spherical Illing et al. (2017); Vivek et al. (2017) and ellipsoidal Zhang and Cheng (2019) particles, as we have shown. However, we have found no clear evidence of LW fluctuations in our overdamped numerical simulations. Hence, the overdamped simulations do not fully describe the behavior of colloidal suspensions. This is not a surprise, as it is indeed well known Weitz et al. (1989); Zhu et al. (1992); Mason and Weitz (1995); Dhont (1996); Li and Raizen (2013) that, due to the presence of hydrodynamic interactions, the velocity autocorrelation function of colloidal systems does not decay exponentially as in the numerical simulations of the overdamped dynamics. The upshot of this consideration is that collective vibrations observed in colloidal systems, including the LW fluctuations, may stem from the hydrodynamic interparticle interaction. It would be of interest to better characterize these collective hydrodynamic induced modes.

IV *Materials and Methods *

IV.1 Model systems

In 2d, we investigated the mKA model Brüning et al. (2009) and the harmonic model O’Hern et al. (2003). The mKA model is a mixture of $65\%$ $A$ and $35\%$ $B$ particles. The potential between particles $\alpha$ and $\beta$ is $U_{\alpha\beta}(r)=4\epsilon_{\alpha\beta}[(\sigma_{\alpha\beta}/r_{\alpha\beta})^{12}-(\sigma_{\alpha\beta}/r_{\alpha\beta})^{6}+C_{\alpha\beta}]$, when $r_{\alpha\beta}$ is smaller than the potential cutoff $r^{c}_{\alpha\beta}$, and zero otherwise. Here, $\alpha,\beta\in$ $\{A,B\}$. The interaction parameters are given by $\sigma_{AB}/\sigma_{AA}=0.8$, $\sigma_{BB}/\sigma_{AA}=0.88$, $\epsilon_{AB}/\epsilon_{AA}=1.5$, and $\epsilon_{BB}/\epsilon_{AA}=0.5$. The potential is truncated and shifted at $2.5\sigma_{\alpha\beta}$. $C_{\alpha\beta}$ guarantees $U_{\alpha\beta}(r^{c}_{\alpha\beta})=0$. The number density is $\rho=1.2$. Length, energy and time are recorded in units of $\sigma_{AA}$, $\epsilon_{AA}$ and $\sqrt{m\sigma_{AA}^{2}/\epsilon_{AA}}$, respectively. For this model, we consider $N=2000$ (if not otherwise stated) and $N=10000$. The Harmonic model O’Hern et al. (2003) is a $50:50$ mixture of $N=3000$ particles with interaction potential $U_{\alpha\beta}(r)=0.5\epsilon(1-r/\sigma_{\alpha\beta})^{2}$, for $r<\sigma_{\alpha\beta}$ and $U_{\alpha\beta}(r)=0$ otherwise. The size ratios are $\sigma_{AB}/\sigma_{AA}=1.2$ and $\sigma_{BB}/\sigma_{AA}=1.4$, and the number density is $\rho=0.699$. The units for energy, length and time are $\epsilon$, $\sigma_{AA}$, and $\sqrt{m\sigma_{AA}^{2}/\epsilon}$, respectively. In 3d, we simulated the KA model Kob and Andersen (1994), which consists of $N=3074$ with $80\%$ $A$ and $20\%$ $B$ particles, as well as the Harmonic model, with $N=3000$. All the results are averaged over at least $4$ independent runs. We show the data for $A$ particles if the system is a binary mixture.

IV.2 Numerical details

We have performed simulations in different thermodynamic ensembles, integrating the equation of motion via a Verlet algorithm. For the mKA and the KA models, we have performed first simulations in the NVT ensemble to thermalize the system, using a Nosé-Hoover thermostat, and then in the NVE ensemble for the production runs. For the mKA model, we have also considered a Langevin dynamics. In this case the equation of motion are $m\mathbf{\ddot{r}}_{i}=-\mathbf{\nabla}_{i}\sum_{j(\neq i)}U(r_{ij})-\gamma\mathbf{\dot{r}}_{i}+\mathbf{\eta}_{i}(t)$, where $\mathbf{r}_{i}$ is the position of the $i$th particle, $r_{ij}$ is the inter-particle distance, $U(r)$ is the potential, $\gamma$ is the friction coefficient and $\mathbf{\eta}_{i}$ is a random noise. The random noise satisfies $\langle\mathbf{\eta}_{i}(t)\mathbf{\eta}_{j}(t^{\prime})\rangle=2k_{B}T\gamma^{-1}\delta_{ij}\delta_{t-t^{\prime}}\mathbf{1}$, with $\mathbf{1}$ the unit tensor. In the Langevin model, the velocity autocorrelation function decays exponentially with a Brownian timescale $\tau_{B}=m/\gamma$, in the dilute limit. If $\tau_{B}$ is the smallest timescale of the problem, this also occurs at finite densities, and the system is overdamped. If $\tau_{B}$ is very large, conversely, the system is underdamped and behaves as in the NVE ensemble, where the decay of velocity autocorrelation is induced by the interparticle collisions.

For the Harmonic models, we have performed simulations in the NVT+NVE ensemble. At very low temperature numerical noise causes a small temperature drift in long simulations, in the supercooled regime. In this case the production runs are also carried in the NVT ensemble.

We performed both Newtonian and Langevin dynamics with the GPU-accelerated GALAMOST package Zhu et al. (2013). All reported data are taken after equilibrating the system for at least $20\tau$.

Results of Fig. 6 are obtained by projecting the normalized particle displacement at time $t$ on the modes of the inherent structures of the $t=0$ configurations. We have obtained these modes minimizing the energy of the $t=0$ configurations using the conjugate gradient method, and then diagonalizing their Hessian matrix.

IV.3 Calculation details

The MSD is $\langle\Delta r^{2}(t)\rangle=\frac{1}{N}\langle\sum_{i=1}^{N}\Delta\mathbf{r}_{i}(t)^{2}\rangle$, where $\Delta\mathbf{r}_{i}(t)=\mathbf{r}_{i}(t)-{\mathbf{r}}_{i}(0)$ is the displacement of particle $i$ at time $t$. Its long time behavior defines the diffusion coefficient $D=\lim_{t\to\infty}\frac{\langle r^{2}(t)\rangle}{2dt}$, with $d$ the dimensionality. The ISF is $F_{s}(q,t)=1/N\langle\sum_{j=1}^{N}e^{i\mathbf{q}\cdot\Delta\mathbf{r}_{j}(t)}\rangle$ with $q=|\mathbf{q}|$ the wavenumber of the first peak of the static structure factor. The relaxation time $\tau$ is such that $F_{s}(q,\tau)=c$, where $c$ is a small constant. We fix $c=e^{-1}$ as customary, but have checked that this choice does not affect our results. The non-Gaussian parameter is defined as $\alpha_{2}(t)=\frac{\langle\Delta x(t)^{4}\rangle}{3\langle\Delta x(t)^{2}\rangle^{2}}-1$ with $\Delta x(t)$ the displacement in the $x$ coordinate Weeks et al. (2000). Finally, the four-point dynamical susceptibility $\chi_{4}(t)$ is defined as the variance of the self-intermediate scattering function: $\chi_{4}(t)=N[\langle\hat{F}_{s}(q,t)^{2}\rangle]-\langle\hat{F}_{s}(q,t)\rangle^{2}]$, with $\hat{F}_{s}(q,t)=1/N\sum_{j=1}^{N}e^{i\mathbf{q}\cdot\Delta\mathbf{r}_{j}(t)}$.

The CR quantities, including CR-MSD, CR-ISF and CR four-point susceptibility are defined by replacing the standard displacement with the CR one, which is defined as $\Delta\mathbf{r}_{i}^{\rm CR}(t)=\mathbf{r}_{i}(t)-\mathbf{r}_{i}(0)-1/N_{j}\sum_{j=1}^{N_{j}}[\mathbf{r}_{j}(t)-\mathbf{r}_{j}(0)]$ with $N_{j}$ the number of neighbors of particle $i$ evaluated at time [math]. We use the Voronoi construction to identify the neighbors of a particle.

The viscosity is computed via the Green-Kubo Hansen and McDonald (2005) integral of the shear stress autocorrelation function, $\eta=\frac{V}{k_{B}T}\int_{0}^{\infty}dt\langle S_{ab}(t)S_{ab}(0)\rangle$, where the stress tensor is $S_{ab}(t)=\sum_{i=1}^{N}mv_{ia}v_{ib}-\sum_{i=1}^{N}\sum_{j>1}^{N}r_{ija}dU(r_{ij})/dr_{ijb}$, with $a,b=x,y,z$ indicating the Cartesian components. To reduce numerical noise, we evaluate the viscosity integrating a numerical fit of the stress autocorrelation function. In the normal liquid regime, $G_{p}/(1+(t/\tau_{p})^{\beta_{p}})$ well describes the stress autocorrelation function, at all times. In the supercooled regime we use a piece-wise approximation, using the above functional form to describe the short time behavior, and a stretched exponential function $G_{e}\exp(-(t/\tau_{e})^{\beta_{e}})$ to describe the long time tails. In any case, the integral is dominated by the long-time decay.

Acknowledgements.

MPC and YWL acknowledge support from the Singapore Ministry of Education through the Academic Research Fund (Tier 2) MOE2017-T2-1-066 (S) and from the National Research Foundation Singapore, and are grateful to the National Supercomputing Centre (NSCC) of Singapore for providing computational resources. KZ acknowledges the support from the National Natural Science Foundation of China (21573159 and 21621004). ZYS acknowledges the support from the National Natural Science Foundation of China (21833008, 21790344). TGM acknowledges financial support from UCLA.