Heterogeneous-elasticity theory of instantaneous normal modes in liquids

Since decades, the concept of vibrational density of states in glasses has been mirrored in liquids by the instantaneous-normal-mode spectrum. In glasses instantaneous configurations are believed to be situated close to minima of the potential-energy hypersurface and all eigenvalues of the associated Hessian matrix are positive. In liquids this is no longer true, and modes corresponding to both positive and negative eigenvalues exist. The instantaneous-normal-mode spectrum has been numerically investigated in the past, and it has been demonstrated to bring important information on the liquid dynamics and transport properties. A systematic deeper theoretical understanding is now needed. Heterogeneous-elasticity theory has proven to be particularly successful in explaining many details of the low-frequency excitations in glasses, ranging from the thoroughly studied boson peak, to other anomalies related to the crossover between wave-like and random-matrix-like excitations. Here we present an extension of heterogeneous-elasticity theory to the liquid state, and show that the outcome of the theory agrees well to the results of extensive molecular-dynamics simulations of a model liquid at different temperatures. We find that the spectrum of eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\lambda )$$\end{document}ρ(λ) has a sharp maximum close to (but not at) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda =0$$\end{document}λ=0, and decreases monotonically with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\lambda |$$\end{document}|λ| on both its stable and unstable side. We show that the spectral shape strongly depends on temperature, being symmetric at high temperatures and becoming rather asymmetric at low temperatures, close to the dynamical critical temperature. Most importantly, we demonstrate that the theory naturally reproduces a surprising phenomenon, a zero-energy spectral singularity with a cusp-like character developing in the vibrational spectra upon cooling. This feature, known from a few previous numerical studies, has been generally overlooked in the past due to a misleading representation of the data. We provide a thorough analysis of this issue, based on both very accurate predictions of our theory, and computational studies of model liquid systems with extended size.


Theory
We consider a system formed by N particles with a potential energy V ({r i (t)} i=1..N ) , which depends on the instantaneous positions r i (t) of the particles at time t.We assume that it can be expressed in terms of pairwise potentials φ(r ij (t)) as, with r ij (t) = |r i (t) − r j (t)| .In the INM procedure one considers virtual small displacements, u i (t) , around the instantaneous particle positions, r i (t) , obtained in a MD simulation.The second-order coefficients of a Taylor expansion of V form the Hessian matrix, (1) www.nature.com/scientificreports/ (α, β = 1, . . ., 3 ), which can be diagonalized to obtain the eigenvalues p , with p = 1, . . ., 3(N − 1) , together with the associated spectrum, ρ( p ) .(Note that 3 eigenvalues are trivially zero, due to the conditions i H αβ ij = 0.) If the system is close enough to a minimum, all the curvatures are positive and so are all the p .These are related to the square of the vibrational frequencies, p = ω 2 p , and one can consider the Hessian as the counterpart of the dynamical matrix of a solid.In this regime one can therefore identify the spectrum with the density of vibrational states (DOS) of the system, If, in contrast, the system is away from a minimum, one or more negative curvatures exist, and the non-trivial eigenvalues of the Hessian are both positive and negative.In this case the identification of p with the square of a "frequency" is obviously no longer possible.For the unstable part of the spectrum ( p ≤ 0 ) it has become common practice to define the positive quantity ω = −i √ , together with the corresponding "DOS" as, This function is usually plotted on the negative ω-axis, and the stable g(ω) , according to Eq. ( 3), on the positive ω-axis.This representation, although commonly used as an alternative to directly looking at the ρ( ) , is misleading as unfortunately hides many crucial details of the INM spectrum.Indeed, the former is obtained from the latter by multiplication by |ω| , a procedure, which strongly suppresses all low-frequency details of the spectrum ρ( ) 34 .Based on this observation, Sastry et al. 6 , and Taraskin and Elliott 7 showed that the INM spectrum of liquids, when plotted directly as ρ( ) exhibits a characteristic cusp-like maximum at = 0 , which cannot be seen in an analysis based on Eqs. ( 3) and (4).
In order to obtain a theoretical description of the spectrum ρ( ) , we start from the equations of motion of linear elasticity with a spatially fluctuating shear modulus, G(r) , similar to HET in glasses.HET has been rather successful in explaining several low-frequeny anomalies in the vibrational spectrum of glasses 35-37, 39, 42 , including the so-called Boson peak, an enhancement of the DOS at finite frequencies over the Debye prediction g(ω) ∝ ω 2 .We now aim to explore the application of this theory in the unstable regime, and compare the results with the INM spectrum of a simulated soft-sphere liquid (see Numerical Results).
We consider the equations of motion for the continuum version of the virtual displacements, u i (t) → u(r, t) .The time derivative of this quantity, v(r, t) = u(r, t) , is the Eulerian instantaneous local velocity field.The HET equation of motion is, where M(r) = K 0 + 4 3 G(r) is the longitudinal modulus.Here, K 0 is the bulk modulus, which is considered not to exhibit relevant spatial fluctuations.In contrast, the fluctuating shear modulus is written as G(r) = G 0 − �(r) , where G 0 is the spatial average of G(r) , and �(r) is a local deviation from the average, which we assume to follow a normal distribution with a standard deviation σ 2 given by In the HET theory for a glass 35-37, 39, 42 , the standard deviation σ is, of course, temperature independent.In the present treatment for a liquid, in contrast, one expects a temperature dependent σ , as the features of the visited regions of the PEL depend on T. By comparing with the simulation data we shall see that σ 2 is proportional to the temperature, a feature already observed earlier in a simulation of a glass-forming liquid at higher temperatures 41 .
To proceed further, we decompose the displacements into the longitudinal ( ∇ × u L = 0 ) and transverse ( ∇ • u T = 0 ) components, and transform to frequency space, leading to the stochastic Helmholtz equations, Here, v 2 L (r) = M(r)/ρ and v T (r) 2 = G(r)/ρ are the locally fluctuating squared sound velocities, longitudinal (L) and transverse (T), respectively, and ρ is the mass density.Considering Gaussian fluctuations for �(r) , we can use the replica trick to calculate the average spectrum 35,43,44 .
In three dimensions, the longitudinal and transverse Green's functions for Eqs.(7) are, where the v L,T (z) are the effective complex frequency-dependent sound velocities, We choose G 0 as the reference shear modulus, which is the average shear modulus in the stable glass.G 0 might be approximately correspond to G ∞ , the liquid high-frequency shear modulus, which appears in the Maxwell relation τ = η/G ∞ between the relaxation time τ and the shear viscosity η 45 .The self-consistent HET-SCBA equation for the self energy �(z) is, with Here V c is the coarse-graining volume, and the upper cutoff, k ξ , is inversely proportional to the correlation length, ξ , of the fluctuations of the local shear modulus G(r).
The real parts of v L,T represent the longitudinal and transverse sound velocities.At high frequencies these quantities relate to well-defined vibrational wave excitations.Their imaginary parts are proportional to the corresponding sound attenuation coefficients 36,39 .As pointed out in Ref. 39 the longitudinal Green's function G L (k, z) is related (within the one-phonon approximation) to the inelastic neutron and X-ray scattering laws.The imaginary part of the function �(z) is proportional to the depolarized Raman spectrum.At low and nega- tive these quantities lose their meaning.The real part of G 0 − �(z) becomes negative, and the imaginary part of �(z) dominates the spectrum, as will be demonstrated now.
We introduce the dimensionless variables, , and the (also dimensionless) disorder parameter, In terms of these quantities, Eq. ( 10) takes the form, where, By further defining the dimensionless squared frequency, z=zρ/G 0 k 2 ξ , and, correspondingly, ˜ = ρ/G 0 k 2 ξ , we can write, with The self-consistent Eqs.(15) and (16), must be solved numerically for the self energy, �(z) , and the Green's functions, GL,T (q, z) .The latter can be used in turn to calculate the dimensionless density of eigenvalues as, where ℑ{. . .} denotes the imaginary part, and we have defined the local Green's functions as (9) (13) �(z) = γ 1 0 dq q 4 2 GL (q, z) + 3 GT (q, z) , The numerical solution of the self-consistent equations provides the density of eigenvalues shown in Fig. 1, for the indicated values of the parameter γ (Eq. ( 12)).We found that to a good approximation the predicted spectra agree to the simulated ones (shown in Fig. 4), if we assume that the disorder parameter γ increases linearly with temperature.In Lennard-Jones units the proportionality constant turned out to be 3.5, i.e., γ = 3.5 T.

The transverse contribution
A simple but robust approximation, at least in the low-frequency region, consists in only considering the transverse contribution to the DOS.This choice is justified by the observation that the longitudinal dynamics contributes to the spectra less than 10 % in the small region, and that the easy-to-follow theoretical equations for the transverse dynamics become more complicated and nontransparent if the longitudinal contributions are included.
Neglecting the longitudinal contribution, the HET-SCBA equation, Eq. ( 13) becomes, which, defining the auxiliary function f (z) , can be also expressed as, It is now convenient to write Eq. ( 21) as, where we have defined, As the function GT (z) reaches a finite value for |z| → 0 , Eq. ( 22) makes clear that the low-frequency limit of the self-energy is also finite.We may now solve the quadratic Eq. ( 22), and obtain, with γ c = 1 4 .(We keep the solution with the plus sign which, if one goes into the glass with γ very small, recovers the Born approximation result 46 .)The unstable (i.e., liquid) situation is characterized by γ > γ c .In this case it is convenient re-write Eq. ( 24) as, The self consistent equations for the transverse contribution to the level density are, therefore, Eqs. ( 23) and ( 25).

The low-frequency limit
It is worth noting that the self energy �(z) satisfies the following properties for vanishing argument, where we have denoted, and by combining both equations, To finalize the list of properties of the self energy at small frequency, we observe that, besides the constant 0 , the leading term is linear, as one can conclude from the relation, The density of eigenvalues , considering the transverse contributions only, is therefore given by, We finally analyze the imaginary part of F ( ) in Eq. ( 23) at low frequency.For small , the argument of the logarithm becomes -1, therefore its real part vanishes while the imaginary part become constant and equal to π,

The zero-frequency value and the cusp
It is now straightforward to obtain the density of eigenvalues at = 0 as, Interestingly, this value is not approached linearly by our theory, but rather in a non-analytic way, producing a characteristic cusp already observed in 7 .Indeed, beside the constant value reached at ˜ =0, the leading term of the function F ( ˜ ) is proportional to ˜ 1/2 .In the low frequency expansion of Eq. ( 31) we can therefore safely substitute �( ˜ ) with 0 , obtaining, We see that, for ˜ small but > 0, while, for negative ˜ we find, ( 25) Summing up, for γ c =1/4 we can write, where the function h + (γ ) ( h − (γ ) ) must be used for the positive (negative) -regions, and are defined as, A non-analytical behavior of the density of level is therefore clear in Eq. (36), with an asymmetric square-root cusp centered at = 0 .The functions h ± (γ ) , determining both direction and sharpness of the cusp, are shown in Fig. 2. We see that for negative values of , we always obtain a upward 1/2 cusp.For > 0 , in contrast, we find a positive rise proportional to 1/2 in the region where γ is small (i.e., at low temperatures), whereas for large γ an upwards 1/2 cusp is obtained.The transition between these two regimes occurs at γ = 1.
The validity of the low-frequency approximation of Eqs. ( 36) and ( 37) is demonstrated in Fig. 3, where the "exact" solution of the self-consistent Eqs. ( 25) and (23)(full line) is compared with the approximated one (dashed line), for the two indicated values of γ .We conclude that the description of the low frequency level density with a cusp is valid and robust.
Prefactor of the 1/2 term in the density of levels ρ( ) (Eq. ( 36)) for positive (blue line) and negative (red line) values of .Only the contributions of the transverse modes to the level distribution are considered here.The prefactor is always negative (thus the cusp is upward) in the negative region.For > 0 the direction of the cusp depends on the value of γ .At γ close to γ c (=1/4), corresponding to low temperatures, the cusp is downward, on increasing γ (in our approximation, at γ=1) the cusp changes direction.The origin of the non-analytic singular behavior discussed above can be traced back to the sum rule j H αβ ij = 0 for the Hessian 7 , due to the global translation invariance.In our elastic model the sum rule is trans- formed to the double spatial derivative in the elastic wave Eq. ( 5), or in the Helmholtz Eqs.(7), again reflecting global translation invariance.Technically, the 1/2 singularity stems from vanishing values of the wave-number k in the integrand of the SCBA Eq. ( 10), corresponding to density and stress fluctuations of very large extent.
Concluding the presentation of the theoretical results, and specifically of the existence of a cusp in the vibrational levels spectrum, we emphasize that this result is an additional evidence of the intrinsic weakness of exceedingly crude approximations sometimes employed in the literature, including the prescription ρ( ) ∼ const (or, equivalently, g(ω) ∼ |ω|) 47,48 .Indeed, in considering this approximation many details of the liquid dynamics are plainly ignored.

Coherent-potential approximation (CPA)
An effective-medium theory, which is superior, compared to the SCBA, is the coherent-potential approximation 49 .The CPA may be also obtained as a saddle-point approximation from a suitable field theory 38 .In the corresponding effective medium the total fluctuating shear modulus G(r) = G 0 − �(r) ≡ G r is replaced by a complex, frequency-dependent modulus, which obeys the self-consistent equation Here Q(z) = G 0 − �(z) , and �(z) is given by Eq. ( 11), and . . .P(G r ) denotes an average over the distribution density P(G r ).
Whereas the SCBA saddle-point relies on γ −1 being a large number, i.e. that the relative fluctuations of G(r) are small, this is not the case for the CPA.(A further advantage of the CPA is that one can use a distribution other than a Gaussian one.)In the small-disorder limit the CPA reduces to the SCBA 38 .As the function �(z) is the same for the SCBA and the CPA, the analytic properties of the spectrum near |z| = 0, which have been discussed in detail in the previous paragraph, are the same.In Fig. 1 we compare the result of the CPA with Gaussian distribution density P(G r ) ∼ exp{−(G r − G 0 ) 2 /2σ 2 } (dashed lines) with those of the SCBA.We observe that near = 0 and for positive the results are almost the same.For negative the CPA is obviously superior, as the spectrum has a tail (like the numerical data, see below), in contrast to the SCBA, which exhibits a sharp lower edge of the spectrum.

Numerical methods
In the following Sections we will compare the above theoretical predictions to the outcome of extensive computer simulations.Here, we give details about the atomic liquid model used, the molecular-dynamics (MD) simulations performed, and the numerical tools employed in the INM analysis.Below, all observables are expressed in Lennard Jones units by setting the fundamental quantities, m, σ , ǫ (mass, length, and energy, respectively), and the Boltzmann constant k B to 1, while σ √ m/ǫ is the unit of frequency.We have considered systems formed by equimolar binary mixtures of N soft-spheres of type A (small) and B (large), respectively, with N A = N B = N/2 .The spheres interact via the pair-wise potential, which is the product of the usual soft sphere potential ( ∝ r −12 ), cutoff and shifted at r c = 1.5 , and a tapering function that ensures continuity of both potential and forces at the cutoff.Here r = r ij is the distance between the particles i and j, ǫ AA = ǫ BB = ǫ AB = 1 , σAA = 0.90 , σBB = 1.09 ≃ 1.2 σAA , and σAB = ( σAA + σBB )/2 ≃ 1 .Masses are m A = m B = m = 1 .All simulations have been conducted by using the high-performance-computing simulation code LAMMPS 50 .
We have considered 7 system sizes, N = 10 3 × 2 n with n = 0, . . ., 7 ( n = 6 excluded), ranging from the situa- tion where the entire set of the Hessian eigenvalues can be calculated comfortably (small N), to that where we are only able to obtain the INM spectrum in a limited frequency range around ≃ 0 (large N), although with high statistical accuracy.We have fixed the simulation box sizes, L, by matching the number density ρ = N/L 3 = 1 , for all systems.We considered 24 values of temperature, T, in the range [0.08 : 1].Starting at T = 1 , the system was slowly annealed to T = 0.08 in steps of δT = 0.04 .At each temperature, a total of 12 × 10 6 time steps were performed, consisting of a 2 × 10 6 thermalization period, followed by the production run with the dumping of the system configurations used in the subsequent analysis.The runs were carried out in the (NVE) ensemble with a numerical integration time step δt = 0.002.
We have implemented the calculation of the Hessian matrix, Eq. ( 2), in JAX 51 , which can automatically differentiate native Python and NumPy functions.H αβ ij has been evaluated at each value of T and N on a number of independent istantaneous configurations (at randomly chosen times t), sufficient to reach the desired statistical accuracy, and subsequently diagonalized to extract the eigenvalues.We have employed the standard scipy.linalg 52inear algebra functions at small values of N, while for N ≥ 16000 we used RALEIGH 53 , a Python implementation of the block Jacobi-conjugated gradients algorithm for computing eigenvalues and eigenvectors of large-scale real symmetric and Hermitian problems in selected eigenvalue ranges.We have eventually estimated the INM spectra, ρ( ) , by constructing histograms of the obtained eigenvalues. , The work-flow we have employed is as follows.First, we integrate numerically the system equations of motion and collect at chosen time steps the coordinates of all particles (configurations).Note that the latter are instantaneous thermal configurations and not inherent configurations, i. e., in general they do not correspond to potential-energy-landscape local minima.The number of collected configurations varies at different system sizes N and temperatures T. We can now determine the Hessian matrix corresponding to each configuration, which is next (completely or partially) diagonalized.Finally, the obtained eigenvalues, , are collected and used to populate the (normalized) histograms ρ( ) .The latter should, therefore, be considered as averaged over the configurations statistical ensemble.

Numerical results
We now present a comparison of the simulated spectra of the soft-sphere liquid with the predictions obtained by the unstable version of HET-SCBA theory, together with an in-depth discussion of the data especially in light of the discussion of the zero-frequency value and the cusp above.In particular, we focus on the finite-size effects on the distribution of the eigenvalues of the Hessian matrix, and on the details of the shape of ρ( ) in different conditions, both in the limit → 0 and for intermediate values of .

Overall INM spectrum
We show in Fig. 4 the calculated INM spectra of the investigated soft-spheres liquids, for the indicated values of T, for N = 4000.The lowest considered temperatures are just above the arrest temperature (we estimate, T MCT ≃ 0.05), while the highest values of T encompass the high-fluidity regime, up to T ≃ 20 T MCT .We observe that the spectrum changes from an almost stable situation with a very small fraction of negative eigenvalues at low T (black symbols), to a distribution which is almost symmetric with respect to positive and negative 's.These T-dependent data are very similar to those predicted by the HET-SCBA shown in Fig. 1 for the indicated values of the disorder parameter, γ .As said above, we chose these values according to γ = 3.5T .While a more detailed comparison of the numerical and analytical ρ( ) curves would give a slight deviation from this linear law, the overall trend is given correctly.We already mentioned that in our earlier simulation 41 , in which the temperature dependence of the variance σ 3 ∝ γ was evaluated, an increase with temperature was observed at the highest considered temperatures.
We turn now to a comparison between the simulated ρ( ) curves (Fig. 4), and those predicted by HET-SCBA and HET-CPA (Fig. 1).The main difference of the simulated and predicted spectra concerns the unstable regime < 0 .In this regime the simulated data exhibit a tail, whereas the SCBA does not.Notably, the CPA spectrum does, indeed, also exhibit a tail, which shows the superiority of the CPA with respect to the SCBA.
Our mean-field HET-SCBA and HET-CPA model description therefore seems to satisfactorily reproduce the numerical INM spectra modifications with T, including the transition from a stable situation at T ≃ T MCT to a quite symmetric distribution at high T, while, as expected due to the normalization, the peak heights decrease with T. We also note that the maximum of the distribution is always found at M > 0 , in contrast to the intuitive expectation M = 0 .We will better discuss this point below.

Finite-size effects
We now quantify the effect due to the finite size of the simulated systems on the ρ( ) .In Fig. 5 we show the details of ρ( ) in the very low-| | region, for the indicated system sizes N, at the lowest investigated temperature, T = 0.08 (corresponding to the black symbols in Fig. 4).We observe a strong N-dependence of the spectral shape, with a slope for < 0 which consistently increases by increasing N. Importantly, the theory predicts a derivative For the investigated soft sphere system, T MCT ≃ 0.05 .All quantities are expressed in LJ units.These data are the numerical counterpart of those shown in Fig. 1.The number of unstable modes decreases significantly upon cooling and, as a consequence, distributions which are quite asymmetric around = 0 at high T, becomes consistently increasingly asymmetric by decreasing T, as expected.These data are discussed at length in the main text.which diverges at =0 (see, for instance, Eq. ( 36)), which cannot be directly observed in our simulation.Our data, however, point in the correct direction, allowing us to qualitatively conclude that our simulation is consistent with the theory in the thermodynamic limit, N → ∞.
To better emphasize this important point, we have extracted the main features of ρ( ) close to = 0 by a linear fit of the data of Fig. 5 restricted to the negative (unstable) -branch, and up to the position of the distribution maximum, always localized at positive values of , as already noted above.We show in Fig. 6 (top) the obtained values for ρ( = 0) , represented as a function of N −1/2 , and for the inverse derivative 1/ρ ′ (0) shown as a function of N −1/3 (bottom).
In the N → ∞ limit, we find that ρ(0) extrapolates to a finite value ≃ 1.686 × 10 −2 (see Eq. ( 32)), while ρ ′ (0) diverges, as expected from the presence of the cusp, Eq. ( 36), for γ close to criticality.HET-SCBA therefore truly succeeds in identifying quantitatively the zero-energy singularity in the simulated energy spectrum.We now demonstrate that it also provides quantitative information about the cusp-like character of the singularity.Here, the possibility of analyzing quite large system sizes is crucial.

The cusp at =0
We are now in the position to analyze quantitatively the shape of the simulated ρ( ) curves at small values of | | on both the stable and unstable regions, and its dependence on temperature or, equivalently, on the disorder parameter, γ .We show in Fig. 7 the details of the ρ( ) at small | | values for the largest investigated system size N =128,000, at the indicated low temperatures.
We conclude that the simulations performed on a soft sphere model confirm both qualitatively and quantitatively the prediction of the INM made on the basis of the HET extended to the liquid case.In particular, we find that the critical behavior of the distributions ρ( ) in the vicinity of the zero-energy spectral singularity is very satisfactorily described by the theory.HET-SCBA seems indeed to adequately grasp the behavior of a realistic model of an atomic fluid.

The slope of ρ( ) at intermediate frequency
An interesting detail of the spectra obtained in numerical simulations, which cannot be easily extracted from analytical calculations, are the strength and sign of the slope of the ρ( ) curves at small but finite | | (thus away from the cusp region), in both the stable and unstable regions.We have estimated both slopes by a single nonlinear fit around = 0 to an empirical function of , ρ( ) = [e −α (c 0 + c 1 ) + e α (c 2 + c 3 )]/(e −α + e α ) .The values of the slopes are directly provided by c 1 and c 3 .We have checked that fitting procedures based on different model functions, including linear models, provide very similar results.
We plot in the left panel of Fig. 9 the estimated slopes as a function of T for N = 4000 (main panel) and 128000 (inset).In the unstable region (circles) the slope is always positive, and consistently decreases upon heating.In the stable region, in contrast, the slope is negative.Interestingly, it increases (in absolute value) on cooling, goes through a maximum at T ≃ 0.2 and eventually decreases approaching T MCT .(Very similar results have been reported in 34 ).The data for large N in the inset provide substantially the same picture, with the additional information that the slope of the stable region becomes positive at T ≃ 0.1, as already noted above.
Unfortunately we do not have an intuition of the mechanism behind the observed behavior of the stable region slope, which we postpone to future work.As more work is needed to extract from the theory a clear explanation  of the data shown in the right panel of Fig. 9. Here, we plot the difference S = S s − S u of the stable and unsta- ble slopes, on a double-logarithmic scale as a function of T for both sets of data.As already observed by some of us 34 , the data follow a T −2/3 law, in the entire investigated temperature range.A satisfactorily explanation of this observation is lacking.

Relations to experimental data
The INM procedure encompasses the snapshot of a liquid, which in principle is not accessible to measurement.However, at high frequencies the dynamics of a liquid is known to be governed by vibrational modes, and the corresponding vibrational DOS, g(ω) , is given by Eq. ( 3).Within the one-phonon approximation g(ω) can be shown 39 to be equal to the Fourier transform of the velocity autocorrelation function.However, in a liquid the ω = 0 value of the latter is known to be equal to the diffusion coeffient 54 .As discussed above, because g(ω) is  The slopes on the stable (squares) and unstable (circles) regions around = 0 have been determined by a simultaneous non-linear fit of the two regions, as discussed in the text.In the unstable region the slope is always positive, and consistently increases upon cooling.In the stable region, in contrast, the slope is always negative, increasing (in absolute value) on lowering T, going through a maximum at T ≃ 0.2 and eventually decreasing.)Left, Inset: Same that in the main panel, for N = 128000 .The behavior is similar but the decreased impact of finite-size effects indicates that the slope of the stable region becomes positive at T ≃ 0.1.Right: Slopes differences S = S s − S u for both values of N, together with the ∝ T −2/3 dependence (dashed lines).All data are discussed at length in the main text.3), which is finite at = ω 2 = 0 , the "DOS" g(ω = 0) is always zero.This can be traced to the breakdown of the one-phonon approximation for liquids, and, equivalently, the breakdown of the concept of a "liquid DOS" 48,55 .At high frequencies the liquid dynamics is similar to that of an amorphous material, which is well described by heterogeneous elasticity theory, and can be measured by inelastic Raman, neutron and X-Ray techniques 39 .The instable part of the spectrum ( < 0 ) has been related in the literature to the mobility, viz. the diffusion constant of the liquid 1,4 .Indeed, a one-to-one correspondence of the fraction of unstable modes to the liquid diffusivity has been established 25,27 .

Conclusions
Melting of a crystal is accompanied by the suppression of sharp features of the vibrational spectrum.For instance, the Van Hove singularities in the density of states, g(ω) , (associated to the critical points of the dispersion, ω(q) ) disappear, and are substituted by the appearance of a broad distribution of vibrational modes, including both stable and unstable excitations.Interestingly, the reverse process can be very different.This is the case for systems that can be supercooled, i.e., where it is possible to avoid crystallization by a succession of equilibrium liquid states extending below the melting temperature.In this process, the high fraction of unstable modes contained in the spectrum of the Hessian matrix of the liquid vanishes at T MCT , which is associated to a dynamical phase transition, and marks the boundary where a saddle-dominated PEL liquid-like kinetics gives way to a minima-dominated flow, eventually leading to glass formation.The limiting shape of the vibrational spectrum approaching T MCT turns out to be surprising, with the appearance of a zero-energy spectral singularity, now unrelated to the ω(q) , with a cusp-like character, and an origin which can be traced back to the most obvious of the symmetries, translational invariance.Despite the generality of this phenomenon also discovered, for instance, in Anderson hamiltonians with offdiagonal disorder 56 , where a similar cusp in the electronic DOS at E = 0 is observed, an analytic description of this phenomenon until now has not been provided.Here, we have addressed in depth the characterization of the zero-energy singularity, by putting together insight coming from theoretical advances and extensive computer simulations.We have presented a generalization to the liquid state of the heterogeneous elastic theory, by integrating in the elastic medium description of the instantaneous normal modes spectrum an explicit temperature dependence of the spatial shear modulus fluctuations.Specifically we have been able not only to demonstrate the existence of a disorder-induced zero-energy spectral singularity, but also fully characterize the disorder dependent spectral shape, in both stable and unstable energy regions.We have next tested the theoretical predictions against extensive Molecular Dynamics simulations of an atomic glass-forming liquid, confirming the relevance of the developed mean-field theory for realistic systems.
In a nutshell, we have demonstrated that focusing to the original INM spectrum ρ( ) reveals much more interesting information than by transforming it to the ω = √ spectrum, which completely suppress most relevant zero-energy features of the spectra.Our findings also pose serious constraints on the allowed shape of ρ( ) employed in theoretical developments, sometimes at variance with superficial recent proposals.

Figure 1 .
Figure1.INM spectra, ρ( ) , of our unstable-elasticity model, calculated in the self-consistent Born approximation (SCBA) (full lines) and in the coherent-potential approximation (CPA) for a Gaussian distribution (dashed lines), as described in the main text.We have considered different values of the disorder parameter γ , all of them on the liquid side ( γ > γ c =1/4), as indicated in the legend.For better comparison with the numerical data of Fig.4we converted the eigenvalues to LJ units according to = 40 ˜ (SCBA) and = 33 ˜ (CPA).

Figure 3 .
Figure 3. Test of the low frequency approximation for the function F ′′ ( ) .The dashed lines are approximate solution from Eqs. (36) and(37), while the full lines are obtained by a numerical solution of the self-consistent Eqs.(25) and(23).The red lines refer to γ =0.5, the green to γ =2.

Figure 4 .
Figure 4. Simulated INM spectra, ρ( ) , for N = 4 × 10 3 at the indicated values of T.For the investigated soft sphere system, T MCT ≃ 0.05 .All quantities are expressed in LJ units.These data are the numerical counterpart of those shown in Fig.1.The number of unstable modes decreases significantly upon cooling and, as a consequence, distributions which are quite asymmetric around = 0 at high T, becomes consistently increasingly asymmetric by decreasing T, as expected.These data are discussed at length in the main text.

Figure 5 .
Figure 5.Effect of the finite system size on the low-INM spectral shape.Simulated INM spectra, ρ( ) , at the lowest investigated temperature, T = 0.08 , for the indicated values of N = 10 3 × 2 n , with n = 0, . . ., 7 ( n = 6 excluded).At the smaller values of N the ρ( ) intensities increase continuously on the unstable side, crossing = 0 and eventually reaching maxima at positive values of .The positions of the maxima, however, continuously shift toward = 0 on increasing N.At the highest investigated values of N the cusp-like non- analytical behavior starts to develop, as discussed into details in the main text.

Figure 6 .
Figure 6.Finite-size scaling of values and derivatives of the INM spectra at = 0 .The data have been determined by fitting linearly the data of Fig. 5 around = 0 .Top: Values of the INM spectra ρ( = 0) at T = 0.08 , as a function of N −1/2 .Data are indicated by the symbols, the solid line is a linear guide for the eyes.Bottom: Inverse of the derivative of the INM spectra at = 0 , represented as a function of N −1/3 .The details of the observed scaling with the system size are discussed in the main text.

Figure 7 .
Figure 7. Details of the low frequency region of the INM spectra.We plot ρ( ) for N=128000, at the four lowest values of T. All details of the data, both in the unstable and stable -regions, are discussed in the main text.

Figure 8 .
Figure 8. Same data of Fig. 7, now plotted as a function of sgn( )| | 1/2 .The consequences of this representation are discussed in details in the main text.

Figure 9 .
Figure 9. Local slopes of ρ( ) close to = 0 , both in the stable and unstable regions.Left, Main panel: Local slope for N = 4000 as a function of temperature, T.The slopes on the stable (squares) and unstable (circles) regions around = 0 have been determined by a simultaneous non-linear fit of the two regions, as discussed in the text.In the unstable region the slope is always positive, and consistently increases upon cooling.In the stable region, in contrast, the slope is always negative, increasing (in absolute value) on lowering T, going through a maximum at T ≃ 0.2 and eventually decreasing.)Left, Inset: Same that in the main panel, for N = 128000 .The behavior is similar but the decreased impact of finite-size effects indicates that the slope of the stable region becomes positive at T ≃ 0.1.Right: Slopes differences S = S s − S u for both values of N, together with the ∝ T −2/3 dependence (dashed lines).All data are discussed at length in the main text.