Modeling the instantaneous normal mode spectra of liquids as that of unstable elastic media.

We study the instantaneous normal mode (INM) spectrum of a simulated soft-sphere liquid at different equilibrium temperatures T We find that the spectrum of eigenvalues [Formula: see text] has a sharp maximum near (but not at) [Formula: see text] and decreases monotonically with [Formula: see text] on both the stable and unstable sides of the spectrum. The spectral shape strongly depends on temperature. It is rather asymmetric at low temperatures (close to the dynamical critical temperature) and becomes symmetric at high temperatures. To explain these findings we present a mean-field theory for [Formula: see text], which is based on a heterogeneous elasticity model, in which the local shear moduli exhibit spatial fluctuations, including negative values. We find good agreement between the simulation data and the model calculations, done with the help of the self-consistent Born approximation (SCBA), when we take the variance of the fluctuations to be proportional to the temperature T More importantly, we find an empirical correlation of the positions of the maxima of [Formula: see text] with the low-frequency exponent of the density of the vibrational modes of the glasses obtained by quenching to [Formula: see text] from the temperature T We discuss the present findings in connection to the liquid to glass transformation and its precursor phenomena.

The investigation of the potential energy surface (PES) of a liquid (made up of N particles with positions at a time instant t) and the corresponding instantaneous normal modes (INMs) of the (Hessian) matrix of curvatures has been a focus of liquid and glass science since the appearance of Goldstein’s seminal article (1) on the relation between the PES and the liquid dynamics in the viscous regime above the glass transition (2–27).

The PES has been shown to form a rather ragged landscape in configuration space (8, 28, 29) characterized by its stationary points. In a glass these points are minima and are called “inherent structures.” The PES is believed to contain important information on the liquid–glass transformation mechanism. For the latter a complete understanding is still missing (28, 30, 31). The existing molecular theory of the liquid–glass transformation is mode-coupling theory (MCT) (32, 33) and its mean-field Potts spin version (28, 34). MCT predicts a sharp transition at a temperature , where T is the temperature of structural arrest (glass transition temperature). MCT completely misses the heterogeneous activated relaxation processes (dynamical heterogeneities), which are evidently present around and below T and which are related to the unstable (negative-λ) part of the INM spectrum (28, 30).

Near and above T, apparently, there occurs a fundamental change in the PES. Numerical studies of model liquids have shown that minima present below T change into saddles, which then explains the absence of activated processes above T (2–24). Very recently, it was shown that T is related to a localization–delocalization transition of the unstable INM modes (25, 26).

The INM spectrum is obtained in molecular dynamic simulations by diagonalizing the Hessian matrix of the interaction potential, taken at a certain time instant t:with . For large positive values of the eigenvalues λ (, N being the number of particles in the system) they are related to the square of vibrational frequencies , and one can consider the Hessian as the counterpart of the dynamical matrix of a solid. In this high-frequency regime one can identify the spectrum with the density of vibrational states (DOS) of the liquid via

For small and negative values of λ this identification is not possible. For the unstable part of the spectrum () it has become common practice to call the imaginary number and define the corresponding DOS as

This function is plotted on the negative ω axis and the stable , according to [2], on the positive axis. However, the (as we shall see, very interesting) details of the spectrum near λ = 0 become almost completely hidden by multiplying the spectrum with . In fact, it has been demonstrated by Sastry et al. (6) and Taraskin and Elliott (7) already 2 decades ago that the INM spectrum of liquids, if plotted as and not as according to [2] and [3], exhibits a characteristic cusp-like maximum at λ = 0. The shape of the spectrum changes strongly with temperature. This is what we find as well in our simulation and what we want to explore further in our present contribution.

In the present contribution we demonstrate that the strong change of the spectrum with temperature can be rather well explained in terms of a model, in which the instantaneous harmonic spectrum of the liquid is interpreted to be that of an elastic medium, in which the local shear moduli exhibit strong spatial fluctuations, which includes a large number of negative values. Because these fluctuations are just a snapshot of thermal fluctuations, we assume that they are obeying Gaussian statistics, the variance of which is proportional to the temperature.

Evidence for a characteristic change in the liquid configurations in the temperature range above T has been obtained in recent simulation studies of the low-frequency vibrational spectrum of glasses, which have been rapidly quenched from a certain parental temperature . If is decreased from high temperatures toward T, the low-frequency exponent of the vibrational DOS of the daughter glass (quenched from to T = 0) changed from Debye-like to with s > 2. In our numerical investigation of the INM spectra we show a correlation of some details of the low-eigenvalue features of these spectra with the low-frequency properties of the daughter glasses obtained by quenching from the parental temperatures.

The stochastic Helmholtz equations (Eq. of an elastic model with spatially fluctuating shear moduli can be readily solved for the averaged Green’s functions by field theoretical techniques (35–37). Via a saddle point approximation with respect to the resulting effective field theory one arrives at a mean-field theory (self-consistent Born approximation [SCBA]) for the self-energy of the averaged Green’s functions. The SCBA predicts a stable spectrum below a threshold value of the variance. Restricted to this stable regime, this theory, called heterogeneous elasticity theory (HET), was rather successful in explaining several low-frequency anomalies in the vibrational spectrum of glasses, including the so-called boson peak, which is an enhancement at finite frequencies over the Debye behavior of the DOS (37–41). We now explore the unstable regime of this theory and compare it to the INM spectrum of our simulated soft-sphere liquid.*

We start Results by presenting a comparison of the simulated spectra of the soft-sphere liquid with those obtained by the unstable version of HET-SCBA theory. We then concentrate on some specific features of the INM spectra, namely, the low-eigenvalue slopes and the shift of the spectral maximum from λ = 0. Both features are accounted for by HET-SCBA. In particular, we find an interesting law for the difference between the slopes of the unstable and the stable parts of the spectrum, which behaves as , which, again, is accounted for by HET-SCBA.

In the end we compare the shift of the spectral maximum with the low-frequency exponent of the DOS of the corresponding daughter glasses and find an empirical correlation. We discuss these results in connection with the saddle to minimum transformation near T.

Results

In Fig. 1 we have plotted the INM spectrum of our soft-sphere liquid (Materials and Methods) for nine temperatures. The lowest ones are just above the glass transition temperature (in Lennard–Jones units [LJ]; Materials and Methods), and the highest one is in the high-fluidity regime . The mode-coupling temperature of the soft-sphere liquid has been reported (26, 42) to be , corresponding to the light violet curves in Fig. 1. We observe that the spectrum changes from an almost stable one at low T to a spectrum, which is almost symmetric with respect to positive and negative λ. In Fig. 1 we show the spectrum predicted by the HET-SCBA (Materials and Methods). We see that HET-SCBA satisfactorily explains the development of the INM spectrum: it changes from a stable one at T to a completely symmetric one at high T, while due to the normalization, the peak height decreases with temperature. Obviously, T qualitatively marks the cross-over from a strongly asymmetric to rather symmetric INM spectrum.

Fig. 1.

(A) INM spectra for different temperatures as indicated in the legend. The eigenvalues and temperatures are in LJ units (Materials and Methods). (B) INM spectra of our unstable elasticity model, calculated in SCBA disorder parameters, assuming for the same values of T as for A. We converted the SCBA units (Materials and Methods) to LJ units.

An interesting detail of the spectra concerns the slope of the curves at small on the stable and unstable sides. In Fig. 2 we plot the positions of these slopes against temperature for both the numerical and the SCBA calculations. We see that they follow the same trends, as can already be seen in Fig. 1. In Fig. 2 we plot the difference of these slopes double-logarithmically against T. Both for the numerical and for the SCBA data we observe a law.

Fig. 2.

(A and C) Small- slopes of the unstable (black squares) and stable (red circles) INM spectra of Fig. 1 as a function of temperature. The lines are guides for the eye. (B and D) Difference of the slopes in double-logarithmic representation. The dashed lines indicate a behavior. (A and B) Slopes of the simulated data and (C and D) slopes of the SCBA spectra.

There is a further interesting detail in which the numerical and HET-SCBA curves differ: the numerical data meet the maximum always with a finite slope, while the SCBA data meet the λ = 0 axis always with a rectangular slope according to (parabolic shape), which implies an inverse square root singularity of the slope at This nonanalytic and singular behavior can be traced to the sum rule of the Hessian (Eq. (7), due to global translation invariance. At a cross-over temperature the positive parabola turns into a negative parabola (7). In our elastic model the sum rule is transformed to the double spatial derivative in the elastic wave Eq. or in the Helmholtz Eq. , reflecting global translation invariance.

Technically, the singularity stems from vanishing values of the wave number k in the integrand of the SCBA Eq. , corresponding to density and stress fluctuations of very large extent. We convinced ourselves that the singularity at λ = 0 disappears if an infrared cutoff k0 is introduced at the lower end of the integral in [8], which mimics the effect of a finite sample size. (The simulated sample is rather small, 1,024 particles.) So, the finite slope at λ = 0 (as compared to the λ = 0 singularity of the SCBA) might be a finite-size effect.

Both in the simulated and the SCBA spectra the maximum is not exactly at λ = 0. In the simulations the maximum position continuously decreases (Fig. 3) and becomes negative at a cross-over temperature . In the HET-SCBA theory, at low enough temperatures, there is also a maximum at due to the upward parabola at λ = 0. For , there is on both sides of the spectrum a downward parabola, and the maximum with the cusp is always at λ = 0.

Fig. 3.

(A) Maxima of , as a function of temperature. (B) Exponent of the daughter glasses, quenched from the parental temperatures against the position of the maxima of the INM spectra at the temperatures, from which the glasses are quenched to . The line indicates a behavior. (C) The exponent vs. . The dotted line is a guide for the eye.

We now turn to the discussion of stable (daughter) glasses, quenched toward 0 from a given parental temperature . We mentioned in the Introduction that it has been shown in ref. 43 that above a certain parental temperature the low-frequency DOS of the glass exhibits Debye-like behavior with , and at lower parental temperatures a continuous increase of s toward the value of is observed.

It should be emphasized that these low-frequency modes are not wave-like modes but diffusive quasi-localized ones (44), which are only visible in very small model samples that do not support elastic standing waves (45).

In Fig. 3 we show the temperature dependence of . In Fig. 3 we have plotted (taken from ref. 43) against the position of the INM-spectral maximum of the corresponding parental liquid. We observe a positive correlation, which looks like ( –2) . At present we have an explanation neither for this empirical law nor for the behavior of the slopes of the INM spectra, recorded in Fig. 2 .

As the non-Debye exponents are associated with nonphononic quasi-localized excitations, we conclude, in agreement with the reasoning in ref. 43, that these modes arise from incomplete relaxation due to the presence of immobile clusters in the parental liquid. The saddle to minimum cross-over in the PES near T with decreasing temperature is here reflected by the increase of the exponent s from 2 to higher values. This increase might be due to a delocalization–localization transition of the modes near λ = 0, inherited from the transition, which occurs in the parent liquid (26).

Discussion

We have demonstrated that looking at the original INM spectrum reveals much more interesting information than transforming it to the spectrum. In studying the INM spectrum of a soft-sphere liquid, we found that it can be well modeled by a spectrum of an unstable elastic medium, in which the shear-elastic modulus has strong spatial thermal fluctuations with variance proportional to the temperature. We found that the maxima of the INM spectra occur at positive values of λ below a cross-over temperature and at negative values above T. In the HET-SCBA model also a finite low-temperature maximum occurs, but not a maximum position in the unstable regime.

In the light of our findings we can sketch the following scenario for the INM spectrum of a liquid with decreasing temperature. At large temperatures the INM spectrum is symmetric: there are as many positive as negative curvatures in the PES. In our model this corresponds to a width of the distribution of shear moduli, which is much larger than the average value μ0 and then leads to the symmetric spectrum. Lowering the temperature toward T, the spectrum becomes increasingly asymmetric. The unstable modes support less and less delocalized states (26). Below T the INM spectrum has a maximum at finite λ, the position of which increases monotonically.

We are convinced that our model description of the INM spectrum will be helpful in the future in order to obtain a better theoretical understanding of the temperature regime of dynamical heterogeneities between T and T.

At the end we would like to comment on a recent paper, which has appeared in PNAS (46) and which deals with a very similar subject. Building on the assumption that the dynamic of the liquid is described by a single-particle local Langevin-type dynamics, the authors of ref. 46 claim that the INM DOS would have a universal shape and, specifically, that at small frequency, , meaning that would be constant in this regime. It is clear from our results that is not constant at small λ and not universal, as it strongly changes with temperature. The authors of ref. 46 claim to have good agreement with recently evaluated INM spectra of Zhang et al. (27). We divided the data of ref. 27 by for the stable (unstable) regime and found not a constant but nonuniversal curves, peaked at λ = 0, similar to ours in Fig. 1. Further, we reanalyzed the Langevin model of these authors and found that a factor is missing in their result for the spectrum. This error occurred probably because they mistook the Green’s function (they do not distinguish between longitudinal and transverse degrees of freedom) for the velocity for that for the displacements , which explains the missing factor . The mathematically correct result gives a density of states (for positive ), which would lead to a singularity of , which is neither observed in the simulations nor predicted by HET-SCBA.

Materials and Methods

Simulation.

For simulating a soft-sphere liquid (which has also the ability of forming a glass) (43, 47), we consider a binary mixture of large and small soft spheres in three dimensions (48). Indicating with r the position of the particle i, with , two particles i, j interact via a pure repulsive potential , where . The potential readswhere we have set the prefactor (LJ). Moreover, we impose a cutoff to the potential at in such a way that for . The coefficient α and β are chosen in a way such that has continuous first and second derivatives at r = r.

σ takes the value σ for the large particles and σ for the small ones, with and . We consider particles (N = N) that are enclosed in a three-dimensional box of side where periodic boundary conditions are employed. The expression for L guarantees , with ρ the particle density.

The normal modes are obtained by means of the diagonalization of the Hessian matrix (Eq. , which can be related to the potential as

For the diagonalization we used the gsl-GNU libraries.

SCBA for the Unstable Elasticity Model.

The equations of motion of linear elasticity with a spatially fluctuating shear modulus can be written aswith the longitudinal modulus . K is the bulk modulus, which is not considered to exhibit spatial fluctuations. We now divide the displacements into longitudinal ones with and transversal ones with and go into frequency space:

with and . ρ is the mass density. We now write . μ0 is the average of . Assuming Gaussian fluctuations of , one can use replica theory to calculate the average spectrum (35–37). From a saddle point approximation of the resulting effective action one obtains an effective medium theory, in which the fluctuating quantity is replaced by a complex, frequency-dependent self-energy with , which obeys the following self-consistent equation [SCBA (37–40)]:

with . γ is proportional to the variance of the fluctuations, which are taken to be proportional to the temperature according to , where γ is in SCBA units (see below) and T is in LJ units.

The spectrum is then given by

The ultraviolet cutoff is given by , where is the correlation length of the fluctuations and C is of order unity. can be estimated to be of the order of one or two interparticle distances. In all calculations we used the value 3.166 for the ratio , which is a realistic value for a soft-sphere glass (39). Further, we used as the unit for squared frequencies. For strong fluctuations (high temperatures) the value of μ0 becomes irrelevant, so equivalently, we measure squared frequencies in , where K is the bulk modulus of the liquid. In Figs. 1 and 2 we converted these units to LJ units by multiplying λ with 800 and, correspondingly, by dividing the spectra by this factor.