The Boson peak and the phonons in glasses

The Boson peak and the phonons in glasses arXiv:cond-mat/0312073v1 2 Dec 2003 S. Ciliberti∗ , T. S. Grigera† , V. Martín-Mayor∗∗ , G. Parisi∗ and P. Verrocchio∗∗ ∗ INFM UdR Roma1, Universitá di Roma “La Sapienza”, and Center for Statistical Mechanics and Complexity (SMC), P.le A. Moro 2, I-00185 Roma, Italy † Centro di Studi e Ricerche “Enrico Fermi”, via Panisperna 89/A, I-00184 Roma, Italy ∗∗ Departamento de Fisica Teorica I, Universidad Complutense de Madrid, Madrid 28040, Spain; Instituto de Biocomputación y Física de Sistemas Complejos (BIFI). Universidad de Zaragoza, 50009 Zaragoza, Spain. Abstract. Despite the presence of topological disorder, phonons seem to exist also in glasses at very high frequencies (THz) and they remarkably persist into the supercooled liquid. A universal feature of such a systems is the Boson peak, an excess of states over the standard Debye contribution at the vibrational density of states. Exploiting the euclidean random matrix theory of vibrations in amorphous systems. we show that this peak is the signature of a phase transition in the space of the stationary points of the energy, from a minima-dominated phase (with phonons) at low energy to a saddle-point dominated phase (without phonons). The theoretical predictions are checked by means of numeric simulations. INTRODUCTION X-ray and neutron scattering techniques allow to obtain very detailed physical insight into the high-frequency (0.1–10 THz) vibrational dynamics of supercooled liquids and glasses. Within this range of frequencies their spectra reveal several universal properties [1], related with the presence of sound-like excitations even for momenta p of the same order of magnitude of p0 , the first maximum of the static structure factor (typically corresponding to wave numbers of a few nm−1 ). This highfrequency sound is revealed as Brillouin-like peaks in the Thz region of the dynamic structure factor. An accessible quantity to experiments is the vibrational density of states (VDOS), g(ω ), whose most striking feature is the presence of an excess of states over the Debye ω 2 law in the “low” frequency region, (i.e.where the dispersion relation is linear, but still in the Thz region)[2]. This excess of states is seen as a peak when plotting g(ω )/ω 2 and has been named Boson peak (BP) 1 . The peak position ωBP usually shifts to lower frequency on heating [3], except for the case of silica [4]. In this material the shift is seen on lowering the density [5]. Due to its universality, the relevant physics underlying the Boson peak can be hopefully captured by some simple model. Furthermore, several recent numerical simu- 1 There exist alternative ways of defining the boson peak from experiments, for example as a peak in Raman scattering data or as a peak in the difference between the observed VDOS of the glass and that of the corresponding crystal lations have shown that a model of harmonic vibrations is wholly adequate to describe this frequency range [6] and that anharmonicity need not be invoked. Given the presence of well formed local structures (SiO2 tetrahedra, for instance) a natural approximation is to consider that the oscillation centers form a crystalline structure, the disorder in the atomic positions being mimicked by randomness in their interaction potential [7, 8] (disordered lattice models [9]). The main drawback of such a models is that they dramatically underestimate the scattering of sound waves [10]. A different approach studies vibrations around a topologically disordered [9] (liquid like) structure. It is followed by two different theories: modified mode-coupling theory [11] (which is not limited to harmonic excitations), and euclidean random matrix theory (ERMT) [12, 13, 14]. ERMT owes its name to the fact that it formulates the vibrational problem as random matrix problem [15]. The matrices involved are called Euclidean random matrices [16], and their study has required the development of new analytical tools. Both MCT and ERMT predict an enhanced scattering of sound waves as compared to disordered crystals. On the other hand, even within the harmonic framework the nature of the extra low frequency modes giving rise to the BP is still an open point. At a qualitative level, the frequency ωBP is close to the Ioffe-Regel [17] frequency ωIR , suggesting the possibility that the excess BP modes are localized [18]. However, numerical simulations have shown that the localization edge is at frequencies greater than ωBP and ωIR [19]. The Ioffe-Regel criterion signals rather a crossover to a region where the harmonic excitations are not longer propagating, due to the strong interaction with the disorder. We call these modes glassons (since they do not propagate but “diffuse”, they have also been called diffusons [19]). A large bump of glassons is generally found around the Ioffe-Regel frequency, due to the flattening of the dispersion relation. This can be considered as the glass counterpart of the van Hove singularity of crystals [8, 20]. All the recently proposed theoretical frameworks predict that this peak of glassons should move to lower frequencies when approaching an instability transition, where negative eigenvalues (imaginary frequencies) appear. The aim of the paper is showing that the ERM theory makes sharp predictions about the values of universal critical exponents describing the approach to this singularity and comparing them with numeric results. The emerging scenario describes the BP modes as given by the hybridization between the phonons and the low-energy tail of the glasson peak which softens when the system approaches the instability transition [21, 22]. THE EUCLIDEAN RANDOM MATRIX THEORY The starting approximation is that particles can only oscillate around fixed random positions, so that the position eq eq of particle i at time t is xi (t) = xi + ϕi (t); the xi are quenched equilibrium positions (whose distribution must be specified) and ϕi (t) are the displacements. Hence the Hamiltonian is 1,N H [x] = ∑ V (xi − x j ) ≃ i, j 1,N 1,3 1 µ ∑ Miµ , jν [xeq]ϕi ϕ νj (1) 2∑ i, j µ ,ν where the dynamical matrix M is an Euclidean Random Matrix: eq Miµ , jν [xeq ] ≡ − f µν (xeq i − x j ) + δi j N ∑ f µν (xeqi − xeqk), k=1 (2) with f µν (x) ≡ ∂µν V (x). In the one-excitation approximation the dynamic structure factor is S(1) (p, ω ) = kB T mω 2 ∑ n 2 ∑ p · en,ieip·xi eq i δ (ω − ωn ), (3) where en are the eigenvectors of the dynamical matrix and ωn its eigenfrequencies (= square root of eigenvalues). The overline means average over the disordered quenched positions, whose distribution P[xeq ] has to be specified. The density of states (VDOS) is obtained in the limit of large momenta: mω 2 (1) S (p, ω ). g(ω ) = lim p→∞ kB T p2 We can obtain S(1)(p, ω ) through the resolvent G(p, z):   eq eq 1 1 ip·(x j −xk ) Gµν (p, z) ≡ ∑ e N jk z − M j µ ,kν   p µ pν p µ pν ≡ GL (p, z) 2 + GT (p, z) δµν − 2 (5) p p separating the axial tensor in a longitudinal term and a transversal one. The dynamic structure factor is then: S(1) (p, ω ) = − (6) A transverse dynamic structure factor (not measurable in experiments) can be defined in an analogous way. However, a most important and general result is that for p → ∞ the resolvent becomes isotropic:   1 1 1 (z) = = δµν Tr [z − M]−1 . G∞ ∑ µν N j z − M j µ , jν N (7) So both longitudinal and transverse structure factors tend to a common limit (the VDOS, see eq. 4) at infinite momentum. 2 Leaving the potential V (r) unspecified and taking the simplified case P[xeq ] = 1/V N (V being the volume), one finds that:   1 Gµν (p, z) = , (8) z − ρ fˆ(0) + ρ fˆ(p) − Σ(p, z) µν The self-energy Σ(p, z) is a matrix with the standard form   p µ pν p µ pν . (9) Σµν (p, z) = ΣL (p) 2 + ΣT (p) δµν − 2 p p which vanishes at ρ = ∞ and that can be computed in a series expansion in 1/ρ . The main point is that the sum of all the infinite diagrams obtained recursively starting from this next-to-leading order result gives a self-consistent integral equation[22]: d3q V (q, p)Gλ σ (q, z)Vσ ν (q, p). (2π )3 µλ (10) where the vertices have the form Vµν (q, p) = ρ ( fˆµν (q)− fˆµν (p − q)). Let us remark that the self energy renormalizes the dispersion relations and gives a finite width to the Brillouin peaks:
2 2 0 ωL,T (p) = ωL,T (p) + Re ΣL,T (p, ωL,T (p)), ΓL,T (p) = Im ΣL,T (p, ωL,T (p))/ωL,T (p). (11) 1 Σµν (p, z) = ρ 2 (4) 2kB T p2 Im GL (p, ω 2 + i0+ ). ωπ Z Conseguently, both the dispersion relations saturate at the same value. However due to the broadening of the line, they are rather illdefined when ω ∼ ωIR The correlations between the equilibrium positions of the particles can be taken into account quite easily at the level of the superposition approximation in the above approach. The results derived above for the case without correlations are translated to the correlated case by replacing the functions f (x) by g(2) (x) f (x). In this way the usual power law divergence of the pair potential for |x| → 0 is balanced by the exponential behaviour of the pair distribution function, and this ensures the existence of the Fourier transform of the product f (x)g(2) (x). The phase transition 2ω Im G∞ (ω 2 + i0+), π d 3 q ˆ2 f (q)G(q, z) (2π )3 (13) R where G∞ (z), A ≡ (2π )−3 d 3 q fˆ2 (q) and the last term are matrices proportional to the identity. The solution of the above integral equation yields a VDOS which contains both the phonons, since g(ω ) ∝ ω 2 at ω → 0, and the extended but not propagating glasˆ sons, described √ by a semicircle with center at ω = ρ f (0) and radius 2 ρ A [20]. If we limit (for pedagogical purposes) to the case where the VDOS changes because of changes in the density, the key point is the existence of a phase transition in the space of the eigenvalues of the Hessian matrix. In fact G∞ (0) develops an imaginary part when ρ < ρc (ρc being a critical density), then the transition separates the stable phase (all positive eigenvalues) and the unstable phase (negative and positive eigenvalues). The glassons are the modes which move towards the negative zone of the spectrum (hybridizing the phonons) when approaching the transition. The order parameter is ϕ = −Im G∞ (i0+ ) which vanishes as ϕ ∼ |∆|β , with β = 1/2 and ∆ ≡ ρ − ρc The relation with the Boson Peak becomes quite clear when one writes down the VDOS in the stable phase arising from the theory without any reference to the control parameter, which then does not need to be ρ . In Z h(x) ∼ x2−γ const. (15) Let us note that the BP is indicated from a peak in the function g(ω )/ω 2 m not in g(ω ). Summarizing, according to ERMT the BP frequency moves linearly toward 0 when approaching the transition (from the stable side) and its height diverges as a power law whose exponent is η = 1/2. Eq. (4) shows that at the level of the onephonon approximation it can also be detected in the large p limit of the dynamic structure factor S(p, ω ). (12) where z = ω 2 + i0+ and G∞ (z) ≡ lim p→∞ G(p, z). Hence one have to solve the integral equation (10) in the p → ∞ limit: 1 = z− ρ fˆ(0)− ρ AG∞ (z)− ρ G∞ (z)  x≪1 , x≫1 (14) with ∆ defined in terms of an arbitrary control parameter.The ERMT (in the cactus approximation) predicts ρ = 1, γ = 3/2. Hence it exist a crossover frequency (in the region where the dispertion relation is still linear) between a ω 2 and a ω γ region. We shall identify that with the BP frequency ωBP . This implies that ωBP ∼ ∆ρ and 2 ∼ ∆−η , with g(ωBP , ∆)/ωBP g(ω , ∆) = ω γ h(ω ∆−ρ ), η = ρ (2 − γ ) . From equation (10) it is possible to derive a few analytic model-independent results about the arising of the Boson Peak. These results are expressed in form of scaling laws, whose exponents are predicted in this approximation. Simulations (see below) and experiments will allow to clarify the dependence of the exponents on the approximation. The VDOS can be obtained from g(ω ) = − fact, one has BOSON PEAK IN A GAUSSIAN MODEL In order to confirm that the saddle-phonon transition described by the Euclidean Random Matrix theory is not an artifact of the approximation involved (cactus resummation), we solved numerically the cactus equation for the case where f (p) has a Gaussian form and compare with direct numerical results for the same model[22]. The model is described by   p µ pν p µ pν , fˆµν (p) = fˆL (p) 2 + fˆT (p) δµν − 2 p p  3/2 2π 2 ˆfL,T (p) = ). (16) exp (−p2 /2σL,T σ02 This choice for fˆ(p) is mainly due to its simplicity. However the behaviour of the Boson peak close to the saddlephonon transition have to be independent of the details of the model. Moreover the superposition approximation takes fˆµν (p) = F [g(r)vµν (r)] yielding a fˆL,T (p) finite at p = 0, like in the Gaussian model. We shall consider various values of the density, which is here the control parameter, comparing the analytical (cactus) results with the numerical spectra and dynamic structure factor obtained from the method of moments [23]. In the high density regime the approximations used in deriving the integral equation (13) are quite good since the analytic solution reproduces the numerical spectrum (and in particular the Debye behaviour) rather accurately (fig.1). numerical analytical 0.008 0.5 0.007 0.4 g(ω)/ω2 g(ω) 0.006 0.005 0.3 0.2 0.1 0.004 0 0 2 4 6 8 10 ω/ω0 0.003 0.002 0.001 0 0 1 2 3 4 5 6 ω/ω0 7 8 9 10 11 FIGURE 1. The VDOS g(ω ) as a function of eigenfrequencies divided by the Debye behaviour ω 2 for ρ = 4 > ρc , both numerical (obtained via the method of moments) and analytical. In the inset, we show g(ω ) vs. ω . However, the crucial check regards the exponents of the transition. Figs. 2a and 2c show that the position of the BP is linear with respect to ∆ ≡ (ρ − ρc ) and that the height of BP diverges as ∆−1/2. This confirms the theoretical predictions ν = 1 and η = 1/2. In Fig. 2b we determine the value of γ by studying the fraction of unstable modes. In fact, in the region of parameters where ρ < ρc the fraction of unstable modes, defined as R0 fu = −∞ gλ (λ )d λ , is given by fu (∆) = Z ∞ 0 d ω ω γ g̃(ω /|∆|) ∼ |∆|1+γ . (17) We find numerically (Figs. 3b) that fu ∼ (ρc − ρ )5/2 , i.e. γ = 3/2. Finally, the order parameter ϕ vanishes as (ρc − ρ )β with β = 1/2 (Fig. 2d). Hence our analytic treatment based on euclidean random matrix theory describe quite well the vibrational features of simple topologically disordered systems[22]. The following step is understanding to what degree of accurateness ERMT could describe the high frequency properties of more realistic systems[21]. BOSON PEAK IN A FRAGILE GLASS Starting from the hypothesis that the Thz region of supercooled liquids and glasses can be described in terms of purely harmonic excitations, the origin of the Boson peak in glasses can be understood if we consider the ensemble of generalized inherent structures (GIS). For each equilibrium configuration the associated GIS is the nearest stationary point of the Hamiltonian. If we start from an equilibrium configuration at low temperature, the GIS is a local minimum, and it coincides with the more frequently used inherent structures (IS) [24](i.e. the nearest minimum of the Hamiltonian). On the contrary, if we start from high temperature, the GISs are saddle points. In the GIS ensemble there is a sharp phase transition separating these two regimes. It takes place in glass-forming liquids [25] at the Mode Coupling temperature [26] (TMC ), above which liquid diffusion is no longer ruled by rare “activated” jumps between ISs but by the motion along the unstable directions of saddles. Phonons are present in the spectrum of the VDOS in the low temperature phase (IS dominated) but are absent in the saddle phase. The key point is that the minima obtained starting from configurations below TMC and the saddles obtained starting above TMC join smoothly at TMC . Thus we can study GIS as a single ensemble parametrized by their energy [25]. Since this transition separates a phase where all the eigenvalues are positive from another one where even negative eigenvalues exist, we expect that ERMT is able to describe correctly this phenomenum. Hence we measured numerically the values of some exponents predicted by the theory in a simple model of a fragile glass [21]. We simulated a soft-spheres binary mixture [27] in the stable (phonon) phase with the Swap Monte Carlo algorithm [28], and computed the VDOS of the ISs obtained starting from equilibrium configurations at temperatures below TMC 3 . In Fig.3 we show that the theoretical predictions agree with the numerical data. Taking the IS’s energy as the 3 At very low T , where equilibrium is not achieved, runs were followed until eIS got very close to its asymptotic value 0.020 data A(ρ-ρc) 0.012 0.008 (a) 0.6 0.7 0.8 0.9 1 1.1 1.2 ρ (b) 0.004 0.000 0.3 0.30 0.35 0.4 ρ 0.45 0.5 0.55 0.35 data 1/2 C+D/(ρ-ρc) 0.25 data 1/2 E(ρc-ρ) 0.30 0.25 0.20 ϕ g(ωBP)/ωBP2 data 5/2 B/(ρc-ρ) 0.016 fu ωBP/ω0 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.15 0.20 0.15 0.10 0.10 (c) 0.05 0.00 0.6 0.7 0.8 0.9 1 1.1 1.2 ρ (d) 0.05 0.00 0.3 0.35 0.4 ρ 0.45 0.5 0.55 FIGURE 2. Numeric results. The critical density in the following fits has been fixed to ρc = 0.54 and capital letters are the fitting parameters. (a): The position of Boson peak as a function of the density near the critical point. ωBP vanishes linearly in ∆ = ρ − ρc . (b): The fraction of unstable modes vanishes as (ρc − ρ )2.5 , thus yielding γ = 3/2. (c): The height of the BP, defined 2 , diverges as ∆−η , with η ∼ 1/2. (d): The order parameter ϕ ≡ −Im G∞ (0) vanishes as (ρ − ρ )β , with β ∼ 1/2. by g(ωBP )/ωBP c CONCLUSIONS In summary, we have shown that the saddle (negative eigenvalues)-phonon (no negative eigenvalues) transition, a common feauture of vibrating topologically disordered systems, is well described by the euclidean random matrix theory. It provides a coherent scenario for the arising of a Boson Peak, which results from the hybridization ωBP/ω0 3.0000 Ec=1.74(1) 2.0000 1.0000 0.0000 0.0050 g(ωBP)/ωBP2 relevant parameter for describing the spectral properties, one has ∆ = ec − eIS , eIS being the energy of the ISs and ec the critical value. In fact, plotting g(ω )/ω 2 a peak is clearly identified, which is seen to grow in height and shift to lower frequency on rising the IS’s energy. Using all the spectra for which the peak position can be clearly identified, we find that the relationship between ωBP and the energy of the IS is linear (Fig. 3a). The energy at which ωBP becomes zero, ec , is found from a linear fit as ec = 1.74 ε (ε is the energy scale), quite close to the value where the GIS stop to be minima (IS) and become saddles [25]. As for the height of the peak (Fig. 3b), the results are compatible with a power-law divergence. Fixing ec at the value 1.74 ε arising from the linear fit of ωBP vs. eIS , a power-law fit yields an exponent β = 0.40(15), while if one fixes the exponent at β = 1/2, then the critical value is ec = 1.752(2) ε . Thus the numerical data are compatible with the theoretically predicted scaling, although we have not been able to work close to the critical point, and thus cannot get a great accuracy on the critical exponents or the critical point. β=0.4(1) 0.0045 0.0040 0.0035 0.0030 1.66 1.67 1.68 1.69 1.7 1.71 EIS FIGURE 3. Scaling of the position, ωBP , and height of the Boson peak near the saddle-phonon transition (energies and frequencies in units of ε and ω0 respectively). Top ωBP is linear in the control parameter eIS and vanishes at eIS = ec = 1.74(1) ε . Bottom The height of the Boson peak diverges as a power law with exponent β ∼ 0.4. Height and position of the BP were obtained by fitting a parabola to the peak of g(ω )/ω 2 . Reprinted with permission from Nature [21], Copyright (2003) Macmillan Magazines Limited of acoustic modes with high-energy modes that soften upon approaching the transition. Hence we applied the theory to describe the saddle-phonon transition and the BP in supecooled liquids, comparing the predicted scaling laws with the numeric results obtained for a simple fragile glass former. The agreement found is quite encouraging. The present discussion applies to experiments as long as one is in the regime where the inverse frequency is much larger than the structural relaxation time, when the harmonic approximation makes sense. We ex- pect that the saddle-phonon transition point of view will be able to bridge the realms of experiment and numerical studies of the energy landscape. As a matter of fact, a recent experiment on the poly(methyl methacrylate) (PMMA) glass gave the first experimental confirmation of the ERMT predictions[29]. 9. 10. 11. 12. ACKNOWLEDGMENTS V.M.-M. is a Ramón y Cajal research fellow (MCyT, Spain). P.V. was supported through the European Community’s Human Potential Programme under contract HPRN-CT-2002-00307, DYGLAGEMEM. 13. 14. 15. 16. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. Sette, F., Krisch, M., Masciovecchio, C., Ruocco, G. and Monaco, G. Science 280, 1550-1555 (1998); Ruocco, G. and Sette, F. J. Phys.: Cond. Matt. 13, 9141-9164 (2001); Courtens, E.; Foret, M., Hehlen, B., Vacher, R., Solid State Commun. 117, 187-200 (2001). Foret, M., Courtens, E., Vacher, R., Suck, J.-B., Phys, Rev. Lett. 77, 3831-3834 (1996); Benassi, P., Krisch, M., Masciovecchio, C., Mazzacurati, V., Monaco, G., Ruocco, G., Sette, F. and Verbeni, R., Phys, Rev. Lett. 77, 3835-3838 (1996); Pilla, O., Cunsolo, A., Fontana, A., Masciovecchio, C., Monaco, G., Montagna, M., Ruocco, G., Scopigno, T., and Sette, F. Phys, Rev. Lett. 85, 2136-2139 (2000); Rufflé, B., Foret, M., Courtens, E., Vacher, R. and Monaco, G., Phys. Rev. Lett. 90, 095502-1 - 095502-4 (2003). Sokolov, A. P., Buchenau, U., Steffen, W., Frick, B., Wischnewski, A., Phys. Rev. B 52, 9815-9818 (1995); Tao, N. J., Li, G., Chen, X., Du, W. M., Cummins, H. Z., Phys. Rev. A 44, 6665-6676 (1991); Engberg, D., Wischnewski, A., Buchenau, U., Börjesson, L., Dianoux, A. J., Sokolov, A. P., Torell, L. M., Phys. Rev. B 59, 4053-4057 (1999). Wischnewski, A., Buchenau, U., Dianoux, A. J., Kamitakahara, W. A., Zarestky, J. L., Phys. Rev. B 57, 2663-2666 (1998). Sugai, S. and Onodera, A., Phys. Rev. Lett. 77 4210-4213 (1996); Inamura, Y., Arai, M., Kitamura, N., Bennington, S. M., Hannon, A. C., Physica B 241-243, 903-905 (1998); Inamura, Y., Arai, M., Yamamuro, O., Inaba, A., Kitamura, N., Otomo, T., Matsuo, T., Bennington, S. M., Hannon, A. C., Physica B 263-264, 299-302 (1999); Jund, P. and Jullien, R., J. Chem. Phys. 113, 2768-2771 (2000). Horbach, J., Kob, W. and Binder, K., J. Phys. Chem. B 103, 4104-4108 (1999); Ruocco, G., Sette, F., Di Leonardo, R., Monaco, G., Sampoli, M., Scopigno, T. and Viliani, G. Phys. Rev. Lett. 84, 5788-5791 (2000). Schirmacher, W., Diezemann, G., Ganter, C., Phys. Rev. Lett. 81, 136-139 (1998). Taraskin, S. N., Loh, Y. L., Natarajan, G. and Elliott, S. R., Phys. Rev. Lett. 86, 1255-1258 (2001); Taraskin S. N. and Elliott, S. R., J. Phys.: Cond. Matt. 14, 3143-3166 (2002); Simdyankin, S. I., Taraskin, S. N., Elenius, M., 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. Elliott, S. R. and Dzugutov, M. Phys. Rev. B 65, 104302-1 - 104302-7 (2002). Elliott, S. R., Physics of Amorphous Materials, Longman, New York (1990). Martín-Mayor, V., Parisi, G. and P. Verrocchio, Phys. Rev. E, 62 2373-2379 (2000). Götze, W. and Mayr, M. R., Phys. Rev. E 61, 587-606 (2000). Mézard, M., Parisi G. and Zee, A., Nucl. Phys. B 559, 689-701 (1999). Martín-Mayor, V., Mézard, M., Parisi, G. and Verrocchio, P., J. Chem. Phys. 114, 8068-8081 (2001). Grigera, T. S., Martín-Mayor, V., Parisi, G. and Verrocchio, P., Phys. Rev. Lett. 87, 085502-1 - 085502-4 (2001). Mehta, M. L., Random Matrices, Academic Press, London (1991). Wu T. M., and Loring, R. F., J. Chem. Phys. 97, 85688575 (1992); Wan Y. and Stratt, R., J. Chem. Phys. 100, 5123-5138 (1994); Cavagna, A., Giardina I. and Parisi, G., Phys. Rev. Lett. 83, 108-111 (1999). Ioffe, A. F. and Regel, A. R., Prog. Semicond. 4, 237-291 (1960); Taraskin S. N. and Elliott, S. R., Phys. Rev. B 61, 12017-12030 (2000). Alexander, S., Phys. Rev. B 40, 7953-7965 (1989). Fabian, J. and Allen, P. B., Phys. Rev. Lett. 77, 3839-3842 (1996); Feldman, J. L., Allen, P. B. and Bickham, S. R., Phys. Rev. B 59, 3551-3559 (1999). Grigera, T. S., Martin-Mayor, V., Parisi, G. and Verrocchio, P., J. Phys.: Cond. Matt. 14, 2167-2179 (2002). Grigera, T. S., Martín-Mayor, V., Parisi, G. and Verrocchio, P, Nature 422, 289-292 (2003). Ciliberti, S., Grigera, T. S., Martín-Mayor, V., Parisi, G. and Verrocchio, P., J. Chem. Phys. 119, 8577-8591 (2003). Benoit, C., Royer E. and Poussigue, G., J. Phys.: Condens. Matter 4, 3125-3152 (1992). Rahman, A., Mandell, M. and McTague, J. P. J. Chem. Phys. 64, 1564-1568 (1976); Seeley, G., and Keyes, T., J. Chem. Phys. 91, 5581-5586 (1989); Madan, B. and Keyes, T., J. Chem. Phys. 98, 3342-3350 (1992); Keyes, T., J. Chem. Phys. 101, 5081-5092 (1994); Cho, M., Fleming, G. R., Saito, S., Ohmine, I. and Stratt, R. M., J. Chem. Phys. 100, 6672-6683 (1994); Bembenek, S. and Laird, B., Phys. Rev. Lett. 74, 936-939 (1995); Bembenek, S., and Laird, B., J. Chem. Phys. 104, 5199-5208 (1996). Angelani, L., Di Leonardo, R., Ruocco, G., Scala, A. and Sciortino, F., Phys. Rev. Lett. 85, 5356-5359 (2000); Broderix, K., Bhattacharya, K. K., Cavagna, A., Zippelius, A. and Giardina, I Phys. Rev. Lett. 85, 5360-5363 (2000); Grigera, T. S., Cavagna, A., Giardina, I. and Parisi, G., Phys. Rev. Lett. 88, 055502-1 - 055502-4 (2002); Götze, W. and Sjorgen, L., Rep. Prog. Phys. 55, 241-376 (1992); Kob, W. and Andersen, H. C., Phys. Rev. E 51, 4626-4641 (1995). Bernu, B., Hansen, J.-P., Hiwatari, Y., and Pastore, G., Phys. Rev. A 36, 4891-4903 (1987). Grigera, T. S. and Parisi, G., Phys. Rev. E 63, 045102-1 045102-4 (2001). Duval, E., Saviot, L., David, L., Etienne, S., Jal, J.F., Europhys. Lett. 63, 778-784 (2003).