# The specific heat of amorphous silica within the harmonic approximation

## I Introduction

It is well known that many low-temperature properties of crystals, such as the specific heat or the vibrational dynamics of the atoms, can be calculated if the wave-vector dependence of the density of states (DOS) is known. At higher temperatures the system will in general become anharmonic and thus the validity of the harmonic approximation breaks down. For amorphous systems at low temperatures the situation is similar to the one of crystals. However, if the temperature is increased not only anharmonic effects have to be taken into account but also the relaxation dynamics of the system since the latter will in general even take place at temperatures significantly below the glass transition temperature $`T_g`$ and thus make the harmonic approximation invalid. That the low temperature specific heat $`C_V`$ of silica can indeed be calculated reliably from the DOS was demonstrated by Buchenau et al. who used the DOS, as determined from neutron scattering, to calculate $`C_V`$ between 5 and 20 Kelvin . Very recently Taraskin and Elliott presented the results of a computer simulation in which they had calculated the specific heat of silica in a similar temperature range . They found that the calculated specific heat is smaller than the experimental one and conjectured that the discrepancy might be due to the relatively small size of their system.

What so far has not been investigated is to what extent the harmonic approximation can be used to calculate the specific heat also at higher temperatures, i.e. in the range 100 K $`TT_g`$, where $`T_g=1450`$ K is the glass transition temperature. In the present work we therefore present the results of a calculation of the specific heat in this temperature range. For this we use the DOS obtained from a large scale molecular dynamics computer simulation and compare the so obtained results with the experimental values.

The rest of the paper is organized as follows: In the next section we give the details of the model used and of the molecular dynamics simulation. In Sec. III we present the results and summarize these in the last section.

## II Model and Details of the Simulation

In this section we give the details of the simulations and describe how the glass samples have been generated.

The silica model we are using in this simulation is the one proposed by van Beest et al. . In this model the interactions $`\varphi (r_{ij})`$ between two ions $`i,j`$ that are separated by a distance $`r_{ij}`$ is given by

$$\varphi (r_{ij})=\frac{q_iq_je^2}{r_{ij}}+A_{ij}e^{B_{ij}r_{ij}}\frac{C_{ij}}{r_{ij}^6}.$$
(1)   
The values of the partial charges $`q_i`$, measured in units of $`e`$, as well as the values of the constants $`A_{ij}`$, $`B_{ij}`$, and $`C_{ij}`$ can be found in Refs. . One of the remarkable features of this potential is that it contains only two-body interactions, thus making it very appealing from a computational point of view, since the evaluation of computationally demanding three-body forces is avoided. Despite its relative simplicity, previous investigations have shown that this model gives a realistic description of the static properties of silica glass  (temperature dependence of the density, structure factor, bond-angle distribution functions), and also of the dynamical properties of silica melts  (diffusion constants, viscosity, intermediate scattering functions).

For the present simulation we use 8016 ions, in a rigid cubic box of size 48.37Å. This corresponds to a density of 2.37 g/cm<sup>3</sup> which is close to the experimental value for amorphous silica and a pressure which is a bit higher than the normal pressure. The reason for choosing a system size that, in the field of simulations of supercooled liquids and glasses, is relatively large, is that the dynamics of strong glass formers  and also the DOS (see below) shows strong finite size effects which have to be avoided.

The equations of motion have been integrated in the microcanonical ensemble with the velocity form of the Verlet algorithm using a time step of 1.6 fs. To generate the glass we proceeded as follows: First we equilibrated the system for 4 million time steps at a relatively high temperature (2900 K). At this temperature the melt is already quite viscous (120 Poise)  and the typical relaxation times are of the order of 3 ns . Subsequently the system was coupled to an external heat bath whose temperature was decreased linearly in time for 1.0 million time steps to zero Kelvin. This corresponds to a cooling rate of about 1.8$`10^{12}`$ K/s. With this cooling rate the system falls out of equilibrium at around 2850 K and thus this temperature corresponds to about the value of the fictive temperature of the glass. During the cooling procedure we stored periodically copies of the sample. These copies were then used as starting configurations for a run in the canonical ensemble with 500,000 time steps (0.8 ns) in order to anneal the configurations. Afterwards (microcanonical) runs were started in which the time Fourier transform of the velocity autocorrelation function was measured which was subsequently used to calculate the density of states. In order to improve the statistics of the results we averaged over four independent runs.

## III Results

In order to calculate the specific heat within the harmonic approximation we need the density of states (DOS), $`g(\nu )`$. To obtain $`g(\nu )`$ one can, e.g., proceed as follows: For any given configuration of particles one uses a steepest descent method to determine the location of the nearest (metastable) minimum of the potential energy. The DOS can then be calculated from the eigenvalues of the (mass-weighted) Hessian matrix at this local minimum. This approach was, e.g., used in Ref.  to determine how the DOS depends on the cooling rate with which the sample was cooled from the high temperature phase to temperatures below $`T_g`$. Below we will come back to these results.

A different method to obtain the DOS is to determine the time Fourier transform of the velocity autocorrelation function and to make use of the relation 

$$g(\nu )=\frac{1}{Nk_BT}\underset{j}{}m_j_{\mathrm{}}^{\mathrm{}}𝑑t\mathrm{exp}(i2\pi \nu t)\stackrel{}{𝐯}_j(t)\stackrel{}{𝐯}_j(0).$$
(2)   
It should be noted that this approach gives an effective DOS which, however, coincides at low temperatures, i.e. when the harmonic approximation is valid, with the real one.

In Fig. 1 we show this effective DOS for three different temperatures, $`T=30`$ K, $`T=300`$ K, and $`T=1050`$ K. Note that all of these temperatures are well below the glass transition temperature ($`T_g`$=1450K) and thus no relaxational processes take place. From the figure we recognize that the main feature of the DOS is a double peak at high frequencies and a relatively flat region at intermediate frequencies, as already discussed in Refs. . It has been shown that the two peaks at high frequencies are due to intra-tetrahedral motions of the ions whereas the frequencies of the inter-tetrahedral motions are in the broad region of $`g(\nu )`$ at lower frequencies . From the temperature dependence of the DOS we see that the main change in $`g(\nu )`$ occurs at high frequencies in that, with decreasing temperature, the height of the two peaks increases and that their location shifts to higher frequencies, i.e. that the corresponding vibrations become faster. (We also note that the curve for $`T=1050`$K shows a tail which extends beyond 40THz. We have checked that this feature is not an effect of the sampling procedure and thus conclude that it originates from the anharmonicity of the system.) At intermediate and low frequencies the temperature dependence of the DOS is much weaker. All these observations are in qualitative agreement with the results of Vollmayr et al. where a similar dependence was found when the glass transition temperature $`T_g`$ was varied . It should be noted, however, that in that work the dependence of $`g(\nu )`$ on $`T_g`$ was related to the fact that the structure of the system changes when the cooling rate is varied. In contrast to this we investigate here only one cooling rate and the temperature dependence of the effective DOS stems only from anharmonic effects and not from a structural change.

We also note that at very small frequencies the DOS shows a gap. As already discussed in Ref.  this is a finite size effect in that excitations which have a spatial extension larger than the size of the simulation box are not present in this system. Some of these excitations are, e.g., acoustic phonons with a large wave length. In the Debye theory the density of states of these phonons is given by

$$g_D(\nu )=3\nu ^2/\nu _D^3$$
(3)   
where the Debye frequency $`\nu _D`$ is related to $`c_t`$ and $`c_l`$, the transverse and longitudinal speed of sound, by

$$\nu _D=\left(\frac{9N}{4\pi V}\right)^{1/3}\left(\frac{2}{c_t^3}+\frac{1}{c_l^3}\right)^{1/3}.$$
(4)   
Here $`N`$ is the number of particles in the volume $`V`$. By measuring the dispersion relation of the transversal and longitudinal acoustic waves at small wave-vectors, we have determined the values of $`c_t`$ and $`c_l`$ . We have found that these values, $`c_t=3772`$m/s and $`c_l=5936`$m/s, agree very well with the experimental values $`c_t=3767`$m/s  and $`c_l=5970`$m/s , thus demonstrating that with respect to these quantities our model is quite realistic<sup>*</sup><sup>*</sup>*Note that these figures are valid for 300 K. Also it has to be taken into account, that the density of our system was fixed at 2.37g/cm<sup>3</sup>, which is slightly higher than the experimental value 2.2g/cm<sup>3</sup>. Hence our sound velocities have to be multiplied by $`\sqrt{2.37/2.2}`$ before they are compared with the experimental values. If this is done one obtains $`c_t=3915`$m/s and $`c_l=6161`$m/s.. For $`\nu _D`$ we obtain 10.65THz which compares well with the experimental value which is 10.40THz . Therefore we have closed the mentioned gap in our DOS by adding to our measured $`g(\nu )`$ the Debye-law given by Eq. (3) in the frequency range $`0\nu 0.73`$THz. We will show below that this modification leads to a significant improvement of the specific heat at low temperatures.

At low temperatures the temperature dependence of the specific heat will be dominated by the low frequency part of the spectrum. This frequency range is also of particular interest because it includes the so-called boson peak, a dynamical feature whose nature is currently a matter of intense debate . As mentioned above, at low frequencies the DOS is expected to show the mentioned Debye behavior, i.e. $`g(\nu )\nu ^2`$. In Fig. 2 we therefore plot $`g(\nu )/\nu ^2`$ for the three different temperatures (bold lines). We see that the DOS shows a large peak at around 1.4THz, i.e. an excess DOS over the Debye value (horizontal dashed line), the mentioned boson peak . Comparing the curves for the three different temperatures, we see that in this frequency range the harmonic approximation holds up to temperatures of 300K, since only the curve for $`T=1050`$K (dotted line) differs significantly from the ones at the two lower temperatures. Note that the curves shown are the one without the mentioned Debye correction. In order to see the effect of this correction we include the DOS for $`T=30`$K when the correction is taken into account (dash-dotted curve). We see that, as expected by construction, this corrected curve goes for small frequencies to the Debye value and joins the measured DOS at 0.73THz.

Also included in the figure are two curves from the work of Vollmayr et al. . These were obtained by cooling a sample of $`N=1002`$ ions with two different cooling rates $`\gamma `$, namely $`\gamma =7.110^{13}`$K/s and $`\gamma =4.410^{12}`$K/s, to zero temperature and calculating the DOS from the dynamical matrix. From these two curves we recognize that the height of the boson-peak depends on the cooling rate, as already discussed in Ref. , and that the low frequency side of the peak of the two curves is at a higher frequency than the one of the curves for 8016 ions. One might be tempted to add the mentioned Debye correction also to the curves for the smaller system, but a look at the magnitude of these corrections shows that they are not sufficient to bring the curves for $`N=1002`$ in coincidence with the ones of $`N=8016`$. Thus we conclude that the excitations giving rise to the boson peak do have a significant spatial extension and thus it is necessary to use large system sizes if finite size effects have to be avoided.

From the DOS it is easy to calculate the temperature dependence of $`C_V`$, the specific heat at constant volume. It is given by

$$C_V=\frac{h^2}{k_BT^2}_0^{\mathrm{}}\frac{\nu ^2\mathrm{exp}(h\nu /k_BT)}{\left(\mathrm{exp}(h\nu /k_BT)1\right)^2}g(\nu )𝑑\nu .$$
(5)   
At low temperatures the specific heat is expected to show a Debye law, i.e. $`C_VT^3`$, and therefore it is useful to plot $`C_V/T^3`$ vs. $`T`$, which is done in Fig. 3. The bold dashed curve is the specific heat as obtained from Eq. (5) using the measured DOS at $`T=30`$K and the solid bold curve is $`C_V`$ from this DOS and the Debye corrections. The horizontal dashed line is the Debye law

$$C_V(T)=\frac{12\pi ^4Nk_B}{5}\left(\frac{T}{\mathrm{\Theta }_D}\right)^3,$$
(6)   
with the Debye temperature $`\mathrm{\Theta }_D=h\nu _D/k_B`$, which for our model is $`\mathrm{\Theta }_D=511`$K. We see that the difference between the uncorrected and corrected curve becomes notiecable for temperatures below 20K in that the former decreases towards zero for $`T0`$ whereas the latter approaches the Debye value given by Eq. (6).

Also included in the figure is experimental data by Sosman and by Zeller and Pohl . We see that although the corrected curve looks qualitatively similar to the experimental data, the location of the peak is overestimated by about 5 K. Furthermore also the height of the peak, as well as the other parts of the theoretical curve, underestimates the measured values by about a factor of two, a result which is in agreement with findings of Taraskin and Elliott . This discrepancy is not due to the presence of two-level systems, since these are relevant only at lower temperatures . Also shown are two curves that stem from the DOS calculated by Vollmayr et al. (see Fig. 2). A comparison of these two curves with the uncorrected one from the present simulation shows that there are at least two possible reasons for the observed discrepancy between the theoretical and experimental curves for the specific heat. First of all we see from the curves for $`N=1002`$ ions that there is a small but noticable cooling rate dependence of $`C_V(T)`$ in that the specific heat increases with decreasing cooling rate. Thus it can be expected that if our silica model would be cooled with a typical experimental cooling rate \[$`O`$(1K/s)\], the theoretical curve would be significantly higher than the one shown in the figure. Secondly we see that the theoretical curves also show a substantial system size dependence in that the one for the larger system is, at low temperatures, significantly above the one for the smaller system. (The fact that the curve for the larger system has a cooling rate which is by a factor of two smaller than the one for the smaller cooling rate of the smaller system is not sufficient to explain this discrepancy between the curves for the two system sizes.) Hence we conclude that the anomalous behavior of the specific heat is also considerably affected by finite size effects.

From the figure we also see that with increasing temperature the relative difference between the theoretical and the experimental curves decreases. Therefore it is interesting to compare these curves also at higher temperatures, which is done in Fig. 4. Here we show $`C_V(T)`$, as calculated from our DOS at 30K (solid line), and $`C_P(T)`$, the specific heat at constant pressure, as determined by various experiments  (lines with symbols). We recognize that for temperatures below the glass transition temperature $`T_g=1450`$K the agreement between our calculation and the experimental data is very good in that the difference is smaller than 0.05 J/gK (see inset). At $`T_g`$ the theoretical curve starts to deviate from the experimental data since at this temperature the real system starts to flow and thus the translational degrees of freedom associated with this motion contribute to the specific heat. Since the harmonic model does not include this type of motion the theoretical curve does not show any special temperature dependence at $`T_g`$ and with increasing temperature approaches the classical value of Dulong and Petit, 1.236J/gK.

Since the experimental data shown in the figure stem from measurements at constant pressure, one has to check whether the discrepancy between the experimental and theoretical curves is due to the difference between $`C_V`$ and $`C_P`$. This difference can be expressed by the thermodynamic relation

$$C_PC_V=TV\frac{\alpha ^2}{\kappa _T},$$
(7)   
where $`\alpha `$ is the thermal expansion coefficient and $`\kappa _T`$ is the isothermal compressibility. Using the experimental values $`\alpha =5.510^7`$K<sup>-1</sup>, valid in the temperature range 293K $`T`$ 593K , and $`\kappa _T=2.7910^5`$GPa<sup>-1</sup> at $`T=1173`$, one obtains an estimate for $`C_PC_V510^3`$J/gK. This is about one decade less than the observed difference between our theoretical values and the experimental data and hence we conclude that the observed (small) discrepancy is due to some inadequacy of the model. (We also mention that from the temperature dependence of the static structure factor at small wave-vectors we have shown that the compressibility of our system is in reasonable good agreement with experimental values, therefore giving support to the above estimate . Finally we mention that the observed discrepancy also does not disappear if instead of the DOS for 30 K we use the one at 300 K or at 1050 K. Hence the difference between the theory and the experiment is indeed due to the used model for silica.

## IV Conclusions

We have presented results from a molecular dynamics computer simulation of amorphous silica in which we used the velocity autocorrelation function to calculate an effective density of states $`g(\nu )`$. We find that for temperatures less than around 300 K $`g(\nu )`$ is essentially independent of $`T`$, i.e. that the harmonic approximation is valid, whereas significant deviations are observed for $`T=1050`$ K. Since this last temperature is still significantly below the experimental galss transition temperature ($`T_g`$=1450 K), and even further below the glass transition temperature of this simulation ($`T_{g,\mathrm{sim}}2850`$K) we do not expect that the system relaxes on the time scale of the simulation. Hence the temperature dependence of $`g(\nu )`$ is due to the anharmonic nature of the local potential wells in which the ions are sitting.

We demonstrate that for frequencies in the vicinity of the boson-peak, i.e. around 1-2THz, the DOS shows a quite strong dependence on the cooling rate, and more important on the system size. Hence it can be concluded that the excitations giving rise to the boson-peak have a relatively large spatial extention so that even our simulation box, which measures about 48Å, is not able to include all of them. Because of this limited system size the low-temperature specific heat we obtain from our simulation is about a factor of two smaller than the one found in experiments. For higher temperatures we find, however, that our theoretical curves agree very well with the experimental data, thus showing that this type of calculation is surprisingly reliable for temperatures below $`T_g`$.

Acknowledgements: We thank Dr. U. Fotheringham for suggesting the calculations presented in this work. This work was supported by BMBF Project 03 N 8008 C and by SFB 262/D1 of the Deutsche Forschungsgemeinschaft. We also thank the HLRZ Jülich for a generous grant of computer time on the T3E.