# The Boson Peak in Amorphous Silica: Results from Molecular Dynamics Computer Simulations

## I Introduction

In the last few years various scattering techniques, such as neutron, Raman and X–ray scattering, have been used to investigate the so–called boson peak, a vibrational feature, which is found in the frequency spectra of many, typically strong, glass formers at a frequency of about $`1`$ THz boson\_peak . In this context various mechanisms giving rise to this peak have been proposed, such as certain localized vibrational modes or scattering of acoustic waves, and also simple models have been developed that produce an excess over the Debye behavior in the density of states boson\_peak\_th .

Especially in the case of silica molecular dynamics computer simulations have recently been used in order to gain insight into the nature of the boson peak taraskin97ab ; dellanna97 ; horbach98 . Despite the limitations of these simulations, such as the small system size (of the order of $`10^3`$$`10^4`$ particles) and high cooling rates (of the order of $`10^{12}`$ $`\mathrm{K}/\mathrm{s}`$), they are very useful because they include in principle the full microscopic information in form of the particle trajectories. Most of the recent computer simulation studies have investigated the boson peak within the harmonic approximation in that the eigenvalues and eigenvectors of the dynamical matrix have been calculated taraskin97ab ; dellanna97 . In contrast to this method we use the full microscopic information to determine quantities like the dynamic structure factor $`S(q,\nu )`$ directly from the particle coordinates. Thus, we are not restricted to the harmonic approximation and we are able to compare the dynamics of our silica model in the liquid state with the dynamics in the glass state. Moreover, by varying parameters like the size of the system and the mass of the particles we gain information on the character of the boson peak excitations.

## II Details of the simulation

The silica model we use for our simulation is the one proposed by van Beest et al. beest  which is given by

$$u(r_{ij})=\frac{q_iq_je^2}{r_{ij}}+A_{ij}\mathrm{exp}(B_{ij}r_{ij})\frac{C_{ij}}{r_{ij}^6}\mathrm{with}i,j\{\mathrm{Si},\mathrm{O}\}.$$
(1)   
The values of the partial charges $`q_i`$ and the constants $`A_{ij},B_{ij}`$ and $`C_{ij}`$ can be found in the original publication. The simulations were done at constant volume keeping the density fixed at $`2.37`$ $`\mathrm{g}/\mathrm{cm}^3`$. Our simulation box contains $`8016`$ particles with a box length of $`48.37`$ Å. We investigate the equilibrium dynamics of the liquid state as well as the glass state. The lowest temperature for which we were able to fully equilibrate our system was $`2750`$ K. At this temperature we integrated the equations of motion over $`13`$ million time steps of $`1.6`$ fs, thus over a time span of about $`21`$ ns. The glass state was produced by starting from two equilibrium configurations at $`T=2900`$ K and cooling them to the temperatures $`T=1670`$ K, $`1050`$ K and $`300`$ K with a cooling rate of $`1.810^{12}`$ $`\mathrm{K}/\mathrm{s}`$. The details of how we calculated the time Fourier transformations can be found elsewhere horbach99 .

## III Results

We investigate the high frequency dynamics of silica by means of the dynamic structure factor

$$S(q,\nu )=N^1_{\mathrm{}}^{\mathrm{}}𝑑t\mathrm{exp}(i2\pi \nu t)\underset{kl}{}\mathrm{exp}(i𝐪[𝐫_k(t)𝐫_l(0)]),$$
(2)   
and its self part $`S_\mathrm{s}(q,\nu )`$ which can be extracted from Eq. (2) by taking into account only the terms with $`k=l`$. In the following we will consider only $`S(q,\nu )`$ for the oxygen–oxygen correlations because the oxygen–silicon and the silicon–silicon correlations behave similarly with respect to the features which are discussed below.

As we have reported elsewhere for the liquid state at the temperature $`T=2900`$ K horbach98 , apart from optical modes with frequencies $`\nu >20`$ THZ, two types of excitations are visible in the dynamic structure factor for $`q>0.23`$ Å<sup>-1</sup>. The first one corresponds to the boson peak which is located, essentially independent of $`q`$, around $`1.8`$ THz. The second one corresponds to dispersive longitudinal acoustic modes. Note that the latter are not like longitudinal acoustic excitations in harmonic crystals because, due to the disorder, they cannot be described as plane waves.

Having found the aforementioned two features at $`T=2900`$ K we will now look at the temperature dependence of $`S(q,\nu )`$, which is shown in Fig. 1, by plotting $`S(q,\nu )/T`$ versus $`\nu `$ at $`q=1.7`$ Å<sup>-1</sup> in the frequency range below $`20`$ THz. From this figure we can conclude that the dynamic structure factor scales roughly with temperature which is expected if the harmonic approximation is valid. Moreover we can clearly identify for all temperatures two peaks: the boson peak located at $`1.75`$ THz (vertical line) and a peak which corresponds to the longitudinal acoustic modes located around $`17`$ THz. Even at $`T=3760`$ K the excitations giving rise to a boson peak at lower temperatures are at least partially present in that a shoulder can be recognized in the frequency region of the boson peak.

In order to get some insight into the properties of the vibrational modes of our silica system at frequencies around $`1`$ THz we varied the size of the system at a fixed mass density of $`2.37`$ $`\mathrm{g}/\mathrm{cm}^3`$. Fig. 2 shows the self part of the dynamic structure factor for $`N=336`$, $`1002`$ and $`8016`$ particles at the temperature $`T=3760`$ K and the three $`q`$ values $`0.37`$ Å<sup>-1</sup>, $`1.7`$ Å<sup>-1</sup> and $`4.75`$ Å<sup>-1</sup>. Whereas the curves for the different system sizes coincide for frequencies that are larger than a weakly $`N`$ dependent frequency $`\nu _{\mathrm{cut}}(N)`$, for $`\nu <\nu _{\mathrm{cut}}(N)`$ the amplitude of $`S_\mathrm{s}(q,\nu )`$ decreases with decreasing $`N`$. Note that $`\nu _{\mathrm{cut}}(N)`$ is essentially independent of the wave–vector $`q`$. We read off $`\nu _{\mathrm{cut}}1.7`$ THz for $`N=336`$ and $`\nu _{\mathrm{cut}}1.2`$ THz for $`N=1002`$. Both frequencies are marked as vertical lines in Fig. 2. $`\nu _{\mathrm{cut}}(N)`$ coincides approximately with the frequency of the transverse acoustic excitation corresponding to the lowest $`q`$ value which is determined by the size of the simulation box. To see that this is the case note that the lowest $`q`$ values for $`N=336`$ and $`N=1002`$ are $`q_{\mathrm{min}}=0.37`$ Å<sup>-1</sup> and $`q_{\mathrm{min}}=0.26`$ Å<sup>-1</sup>, respectively. In comparison to that the $`q`$ values we read off from the transverse acoustic dispersion branch for $`T=3760`$ K and $`N=8016`$, and which correspond to the frequency $`\nu _{cut}`$ for $`N=336`$ and $`N=1002`$, are $`0.32`$ Å<sup>-1</sup> and $`0.22`$ Å<sup>-1</sup>, respectively (see horbach98 ). That the latter $`q`$ values are slightly smaller than the corresponding values for $`q_{\mathrm{min}}`$ is due to the fact that the transverse dispersion branch has been determined from the peak maxima $`\nu _{\mathrm{max}}(q)`$ in the transverse current correlation function. So there is always a significant contribution of transverse acoustic modes with frequencies $`\nu <\nu _{\mathrm{max}}(q)`$ for a given $`q`$. All this can be summarized by saying that the absence of transverse acoustic modes is connected with a missing of excitations giving rise to the boson peak. Therefore, in the smaller systems only the high frequency part of the boson peak is present. Moreover, it seems that the boson peak modes are only fully present, for a given frequency, if there exist transverse acoustic excitations at the same frequency.

To learn more about the character of the boson peak excitations we varied also the masses of the silicon and oxygen atoms such that the mass density remains fixed. Fig. 3 shows $`S(q,\nu )`$ at $`T=2750`$ K and $`q=0.6`$ Å<sup>-1</sup> for the four mass pairs $`M_1=(28.086,15.999)`$, $`M_2=(14.043,23.021)`$, $`M_3=(44.085,8.000)`$, and $`M_4=(56.085,2.000)`$ where the first and the second number are the masses in atomic units for silicon and oxygen, respectively. Note that $`M_1`$ corresponds to the real masses of silicon and oxygen normally used in our simulation. From the figure we see that there is a strong dependence on the mass ratio for the two peaks visible for $`M_1`$ above $`20`$ THz which are due to localized optical modes. In contrast to that, at least within the accuracy of the statistics of our data, there is no dependence for the modes giving rise to the boson peak and the acoustic modes which means that the boson peak excitations cannot be strongly localized. This supports the aforementioned statement that the boson peak is due to the coupling to transverse acoustic modes. The fact that the boson peak is independent of the mass ratio between the silicon and oxygen mass is also a good test for theoretical models of the boson peak in silica.

In conclusion we investigated the excitations giving rise to the boson peak by means of molecular dynamics computer simulations. We find that the dynamic structure factor $`S(q,\nu )`$ scales roughly with temperature in the range $`3760`$ K $`T300`$ K, which means that our silica system is in this sense quite harmonic even for temperatures as high as $`3760`$ K. By calculating $`S_\mathrm{s}(q,\nu )`$ for different system sizes we find that the modes contributing to the boson peak are only fully present at a given frequency if there exist transverse acoustic modes at the same frequency. This is supported by the fact that the height and the width of the boson peak are independent under the variation of the masses of silicon and oxygen if the mass density is fixed. So we observe that the boson peak is due to a coupling to transverse acoustic modes. Of course, the explanation of the nature of this coupling is an interesting goal for the future.

## ACKNOWLEDGMENTS

This work was supported by BMBF Project 03 N 8008 C and by SFB 262/D1 of the Deutsche Forschungsgemeinschaft. We also thank the RUS in Stuttgart for a generous grant of computer time on the T3E.