VOLUME 80, NUMBER 6 P H Y S I C A L R E V I E W L E T T E R S 9 FEBRUARY 1998 High Frequency Sound Waves in Vitreous Silica R. Dell'Anna,1 G. Ruocco,1 M. Sampoli,2 and G. Viliani3 1Universitą di L'Aquila and Istituto Nazionale di Fisica della Materia I-67100, L'Aquila, Italy 2Universitį di Firenze and Istituto Nazionale di Fisica della Materia I-50139, Firenze, Italy 3Universitą di Trento and Istituto Nazionale di Fisica della Materia I-38050, Povo, Trento, Italy (Received 22 July 1997) We report a molecular dynamics simulation study of the sound waves in vitreous silica in the mesoscopic exchanged momentum range. The calculated dynamic structure factors are in quantitative agreement with recent experimental inelastic neutron and x-ray scattering data. Analysis of the longitudinal and transverse current spectra allows us to discriminate between opposite interpretations of the existing experimental data in favor of the propagating nature of the high frequency sound waves. [S0031-9007(98)05310-1] PACS numbers: 61.43.Bn, 61.43.Fs, 63.50.+x In the last 20 years great effort has been devoted to The energy of these excitations extends both below and understanding the microscopic dynamics in topologically above the boson peak region, located at E 4 meV disordered systems. The need originates from the experi- [7]. Recently this result has been confirmed in the whole mental observation that a large variety of glasses 1800­300 K temperature range [12]. At variance with this exhibits common dynamical and thermodynamical prop- interpretation, analogous room temperature IXS data com- erties, which are anomalous with respect to those of the bined with nonkinematic inelastic neutron (INS) spectra corresponding crystals [1]. In particular, around 5 K the [13] were interpreted as an indication of the existence of a thermal conductivity has a plateau and in the 10­30 K crossover from sound waves to localized acoustic modes, temperature range the reduced specific heat CP T T3 which should be responsible for the boson peak. goes through a maximum, indicating the existence of In order to clarify the propagating nature of low fre- excess vibrations in the low frequency range with respect quency modes around the boson peak region, this Letter to the Debye model. In a more extended range, below presents a molecular dynamics (MD) and normal mode 10 meV, inelastic incoherent neutron [2] and Raman analysis (NMA) study of the high frequency dynamics in [3] scattering reveal a strong broadband, usually called vitreous SiO2 in the Q 1.4 5.1 nm21 range. The reli- boson peak, which can again be related to an excess of ability of our investigation is validated by the good agree- the vibrational density of states. At present, a unique ment of the calculated S Q, E with the available INS [8] microscopic description of these excess modes does not and IXS [7] data. The detailed knowledge of the spectral exist, and their theoretical interpretations include local- shape allows us to discriminate between the two proposed ized vibrational states [4], soft anharmonic vibrations [5], interpretations [7,13] and to establish the propagating na- relaxational motions [2], and intrinsically diffusive [6] ture of both longitudinal and transverse modes giving rise and propagating [7] modes. to the boson peak. Moreover, the widths of the excitations Among disordered systems, vitreous silica is a typical calculated in the harmonic approximation (i.e., via NMA) strong glass former and has widely been studied, both ex- are in quantitative agreement with both MD and experi- perimentally [2,3,7,8] and numerically [9]. The tempera- mental data, indicating that in the examined Q range their ture analysis of the inelastic neutron-scattering intensity origin is structural rather than dynamical. [2] suggested the coexistence of sound waves, below Standard microcanonical MD simulations with periodic 4 meV, with a second class of local harmonic excita- boundary conditions were carried out for systems consist- tions, which become increasingly anharmonic towards the ing of 648 and 5184 particles enclosed in cubic boxes lowest frequencies. However, kinematic reasons did not of length L 2.139 nm and L 4.278 nm, respectively (and do not yet) allow the use of the neutron-scattering (mass density r 2.20 g cm3). We used the two- and technique to directly study the acoustic excitations in three-body interaction potential proposed by Vashista the boson peak energy region. Only recently it became et al. [14] which, compared with the simpler- to-deal two- possible to utilize the inelastic x-ray scattering (IXS) with body potentials, gives a better description of the struc- meV energy resolution [11] to measure the dynamic struc- tural properties of SiO2. Details of this study will be ture factor, S Q, E , in the first pseudo-Brillouin zone. discussed in a subsequent paper. The electrostatic long Using this technique, the high frequency dynamics of range interactions were taken into account by the tapered y-SiO2 at T 1050 K has been studied in the 1 6 nm21 reaction field method, and the equations of motion were momentum transfer range. The measured S Q, E clearly integrated with "leap-frog" algorithm using a time step indicates the existence of collective propagating excita- of Dt 0.5 fs. A b-cristobalite crystal was melted at tions up to Q 4 nm21, i.e., up to an energy of 13 meV. 15 000 K and equilibrated for a long time so that the 1236 0031-9007 98 80(6) 1236(4)$15.00 © 1998 The American Physical Society VOLUME 80, NUMBER 6 P H Y S I C A L R E V I E W L E T T E R S 9 FEBRUARY 1998 initial state had no effect on the final configuration. The give identical S Q, E in their common Q range, indicat- system was then cooled and thermalized through differ- ing negligible size effects. ent temperatures ranging from 5000 K down to 300 K. Figure 1 compares in a very satisfactory way INS data The temperature dependence of the total internal energy (squares) of Ref. [8] measured at 5.4 meV (corresponding and of the self-diffusion constants DSi, DO are studied to the position of the boson peak at T 1100 K) with to recognize the structural arrest, which is found at about the neutron-weighted dynamic structure factor SN Q, E 2700 K, a value definitely higher than the experimental calculated by our MD approach at the same temperature. glass-transition temperature [15] Tg 1450 K . This We have also calculated the longitudinal current spec- can be an effect, at least partially, of the quench rate [16], tra CL Q, E , defined as CL Q, E E2S Q, E Q2, at which in our MD simulation is much higher 1013 K s selected small Q values (QMIN 1.47 nm21 for the sys- than the experimental value 102 K s . Two different tem with L 42.78 A). These spectra, reported in Fig. 2 techniques were adopted to study the dynamic of the sys- for some low-Q values, have been fitted by the same tem in the glassy state: (i) (MD) the finite temperature model function [damped harmonic oscillator (DHO)] al- (T 1500, 1000, 600, and 300 K) atomic dynamics has ready used [7] to analyze the experimental spectral shape. been followed for 60 ps by recording the instantaneous In Fig. 3 we report as full circles the best fit values of configurations every 0.01 ps; (ii) (NMA) the dynamical the excitation energies V Q , corresponding to the po- matrix has been calculated in a "glassy" (metastable) en- sitions of the current spectra maxima, and of the ex- ergy minimum, reached by the steepest-descent method citation widths G Q . The form of the current spectra [17], and then either diagonalized to obtain a complete set and the positions of the current spectra maxima are in of eigenvalues and eigenvectors or used in the application qualitative agreement with those of [20]. The linear dis- of the spectral moments method [18,19]. persion observed for V Q gives a longitudinal sound We investigated the dynamical correlation in y-SiO2 velocity which is slightly higher than that compatible with in terms of neutron-weighted SN Q, E and x-ray­ IXS data [7]; this could be due to the thermal history and weighted SX Q, E dynamic structure factors. The quench rate of the numerical model, which can produce presence of two distinct atomic species requires the a tempering of the vitreous structure. The evident linear calculation of the partial dynamic structure factors Q dependence of the current spectra maxima positions, as Z dt already observable from the raw data of Fig. 2, indicates Sab Q, E e2ivt the propagating nature of the excitations up to 20 meV, 2 2p * + a range which includes the boson peak. Na X Nb X eiQ?rj t e2iQ?ri 0 3 p p , This finding is not in contrast with earlier INS re- sults by Buchenau et al. [21], which demonstrate that i[a j[b Na Nb pure plane-wave eigenvectors cannot explain the Q de- where summations over i and j run over all atoms of type pendence of the S Q, E at E 4 meV. Indeed, as it was a and b, respectively, and rj t is the position vector shown also in other simulated glasses [22], the eigenvec- of atom j at time t. Sab Q, E is easily evaluated from tors of propagating modes present a random uncorrelated the MD trajectories. In the NMA one introduces the component superimposed to the plane wave. displacements ui from the equilibrium positions Xi and expands them in terms of the normal modes; in the Q ? ui ! 0 limit the one excitation contribution can be calculated as X d v 2 v Na X Nb X e2Wij S l ab Q, E kBT p l v2 i[a j[b NaNbmamb 3 Q ? ei l Q ? ej l eiQ? Xi2Xj , (1) where the Bose-Einstein factor has been approximated by its classical limit,Wij is the Debye-Waller factor, and vl and ei l are the eigenvalues and eigenvectors of the dy- namical matrix of the system. SN Q, E and SX Q, E are then obtained by adding the partial dynamic structure fac- tors weighted by the atomic species concentrations and by the corresponding neutron or x-ray scattering lengths. An average over the independent directions of Q, i.e., Q 6l, 6m, 6n 2p L, with l, m, n integer numbers, is also performed. The reliability of the harmonic approxi- mation is probed by the good agreement existing between FIG. 1. Q dependence of the inelastic neutron scattering from vitreous silica at 1.3 THz. The circles represent the data of the S Q, E calculated via NMA and via MD at finite tem- Ref. [8] at T 1104 K, the full line the results of the present perature. Moreover, the systems with 648 and 5184 atoms NMA analysis multiplied by an arbitrary factor. 1237 VOLUME 80, NUMBER 6 P H Y S I C A L R E V I E W L E T T E R S 9 FEBRUARY 1998 FIG. 3. (a) Excitation energy V Q from the DHO model for the IXS data at T 1050 K (open circles) [7] and for the present NMA on the 0 K configuration. (b) As in (a) but for the full width at half maximum of the excitations G Q . quantitatively extrapolates over two decades in exchanged momentum down to the Brillouin light scattering (BLS) range Q 0.02 nm21 , suggesting a common structural origin of the S Q, v linewidth, at T , 100 K the G Q measured by BLS shows a strong temperature dependence [23]. Consequently, a structural origin of G Q at high Q values is consistent with available data, while at small Q the dynamical effects could be predominant. At present FIG. 2. X-ray longitudinal (line plus open circles) and trans- we have no explanation either for such different behavior verse (line plus solid circles) current spectra obtained by the in different Q ranges or for the G Q Q2 law observed NMA at different Q values. Open and solid triangles refer to at T . 300 K in the Q 0.02 5 nm21 range. longitudinal and transverse spectra, respectively, obtained by MD (216 SiO In Ref. [13], the localized nature of the high frequency 2 units at 300 K). Full lines are the DHO fits of NMA longitudinal current spectra. excitations was derived from the line-shape analysis of the IXS data. However, the presence of a strong elastic contribution and the poor statistics did not allow one There is quantitative agreement between the G Q to discriminate between different line-shape models. A values obtained from our simulation and that measured detailed line-shape analysis can be performed on the by IXS for temperature ranging from 300 to 1800 K [12], simulated CL Q, E spectra. In the inset of Fig. 4 we both having a G Q ~ Q2 behavior. The agreement of the report an example of fits of our calculated longitudinal G Q values found in our harmonic approximation with current spectra with the expression proposed in Ref. [13]. those measured in the real and simulated system strongly Reliable fits can be obtained only leaving the important indicates that the excitations linewidth, in the investigated vIR parameter free. As reported in Fig. 4, the crossover Q range, is likely to be related to the structural disorder frequency vIR that according to the localization model rather than to anharmonicity or other dynamical effects. should indicate the localization edge and therefore the As suggested in Ref. [7], the structural origin of the boson peak energy is strongly Q dependent and assumes high frequency excitations widths is also supported by values definitely higher than that of the boson peak. If the temperature independence of G Q derived from IXS we impose the same value of the vIR parameter for the data. The picture becomes more complicated at small different Q spectra, the fits result unacceptably bad. Q. Although the G Q Q2 law found at temperature To investigate the nature of the propagating excitations above 300 K by IXS and MD/NMA in the nm21 range giving rise to the boson peak, we have also calculated the 1238 VOLUME 80, NUMBER 6 P H Y S I C A L R E V I E W L E T T E R S 9 FEBRUARY 1998 which is inconsistent with the model itself. (iv) The transverse and longitudinal dynamics keep their differ- ences up to at least Q 5 nm21, which indicates the nonlocalized nature of the vibrational states up to this Q value. We gratefully acknowledge F. Sette for many discus- sions and suggestions, P. Benassi and V. Mazzacurati for useful discussions, and U. Buchenau and S. R. Elliot for the preprints of Ref. [8] and [10], respectively. [1] S. R. Elliot, Physics of Amorphous Materials (Longman, New York, 1990). [2] U. Buchenau, H. M. Zhou, N. Nucker, K. S. Gilroy, and W. A. Phillips, Phys. Rev. Lett. 60, 1318 (1988). FIG. 4. Q dependence of the crossover frequency v [3] G. Winterling, Phys. Rev. B 12, 2432 (1975). IR result- ing from the fits of the localized modes model of Ref. [13] to [4] R. Orbach, Philos. Mag. B 65, 289 (1992). the longitudinal current spectra, calculated from NMA. The in- [5] U. Buchenau, Yu. M. Galperin, V. L. Gurevich, D. A. set shows an example of the fit quality. Parshin, M. A. Ramon, and H. R. Schober, Phys. Rev. B 46, 2798 (1992). [6] P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12 581 transverse current spectra CT Q, E at selected Q values. (1993). The transverse spectra are not experimentally accessible [7] P. Benassi, M. Krisch, C. Masciovecchio, V. Mazzacurati, and are obtained from the MD via the correlation function G. Monaco, G. Ruocco, F. Sette, and R. Verbeni, Phys. of the transverse current or from NMA substituting Q? Rev. Lett. 77, 3835 (1996). in Eq. (1) with Q3 . The transverse current spectra are [8] A. Wischnewsky, U. Buchenau, A. J. Dianoux, W. A. also reported in Fig. 2, where we observe the existence of Kamitakahara, and C. Zaretsky (to be published). a low energy excitation (E 4 meV at Q 1.47 nm21) [9] See, for example, [10], and references therein. [10] S. N. Taraskin and S. R. Elliot, Phys. Rev. B 56, 8605 which disperses with a sound velocity of 4300 m s up (1997); Philos. Mag. B (to be published). to 15 meV at Q 5 nm21. The transverse modes be- [11] F. Sette, G. Ruocco, M. Krisch, U. Bergmann, C. come almost nondispersing at E 20 meV, an energy Masciovecchio, V. Mazzacurati, G. Signorelli, and that is definitely higher than that characteristic of the bo- R. Verbeni, Phys. Rev. Lett. 75, 850 (1995). son peak. Also the propagating transverse dynamics may [12] C. Masciovecchio et al. (to be published). therefore contribute to the boson peak. It is worth noting [13] M. Foret, E. Courtens, R. Vacher, and J. Suck, Phys. Rev. that, in the Q range where the longitudinal and transverse Lett. 77, 3831 (1996). modes show a propagating character, the spectra C [14] P. Vashista, R. K. Kalia, and J. P. Rino, Phys. Rev. B 41, T Q, E and C 12 197 (1990). L Q, E show a marked difference, again indicating the nonlocalized nature of the vibrational modes. Indeed, [15] C. A. Angell, J. Chem. Phys. Solids 49, 863 (1988). for localized modes one should not expect any difference [16] K. Vollmayr, W. Kob, and K. Binder, Phys. Rev. B 54, 15 808 (1996). between "longitudinal" and "transverse" dynamics, being [17] F. H. Stillinger and T. A. Weber, Phys. Rev. A 28, 2408 Q meaningless. (1983). In conclusion, the numerical study of harmonic model [18] C. Benoit, E. Royer, and G. Poussigue, J. Phys. Condens. of y-SiO2 indicates that (i) The experimentally observed Matter 4, 3125 (1992). width of the excitations in the Q 1.4 3.5 nm21 range [19] G. Viliani, R. Dell'Anna, O. Pilla, G. Ruocco, G. is due to structural effects, i.e., to the fact that in this Signorelli, and V. Mazzacurati, Phys. Rev. B 52, 3346 range Q is no longer a good quantum number in the (1995). topologically disordered silica. Anharmonicity or other [20] J. Horbach, W. Kob, and K. Binder, J. Non-Cryst. Solids dynamical effects (interaction with two level systems, (to be published). soft potential modes, etc.) can only play a minor role. [21] U. Buchenau, M. Prager, N. Nucker, A. J. Dianoux, (ii) Both longitudinal and transverse modes are propa- N. Ahmad, and W. A. Phillips, Phys. Rev. B 34, 5665 (1986). gating in the energy range covered by the boson peak: [22] V. Mazzacurati, G. Ruocco, and M. Sampoli, Europhys. their dispersion relations saturate at much higher energy. Lett. 34, 681 (1996). (iii) The explanation of the numerical data with a [23] R. Vacher and J. Pelous, J. Non-Cryst. Solids 45, 397 crossover model gives a Q dependent crossover energy (1981). 1239