Lattice dynamics of α-cristobalite and the Boson peak in silica glass

Cristobalite, SiO₂, crystallizes at high temperature to a cubic structure ($\beta$-phase) and undergoes a phase transition to its $\alpha$-phase upon cooling. $\alpha$-cristobalite has space group P4₁2₁2 and is stable at room temperature and ambient pressure [1]. Natural crystals are usually twinned [2]. The lattice dynamics of $\alpha$-cristobalite has attracted considerable interest due to its auxetic behavior at molecular level [3–5] and the eventual realization of a true surface wave occurring at a single isolated point within a bulk band [6]. It furthermore exhibits a mass density (2.289(1) g cm⁻³) very comparable to ambient silica glass (2.20(1) g cm⁻³) and very similar thermodynamics properties [7]. It could in fact be shown that the excess states over the Debye law in the vibrational density of states (VDOS) in absolute numbers are very similar for $\alpha$-cristobalite and ambient silica glass. The same is true for $\alpha$-quartz and densified silica glass with matched densities. These excess states, which give rise to the so-called 'Boson peak' in glassy systems, are located at the same energies and provide the same heat capacity as in crystals with matched densities. These observations strongly suggest that, similar to crystals, the excess of vibrational states in glasses originates from the piling up of the acoustic-like branches near the boundary of the pseudo-Brillouin zone. A detailed study of the lattice dynamics of $\alpha$-cristobalite is thus an essential prerequisite to finding analogies and gaining further insight into the atomic dynamics and thermodynamics of ambient silica glass. So far, the lattice dynamics of $\alpha$-cristobalite has been investigated by inelastic neutron, Brillouin, and Raman scattering. Phonon frequencies at the $\Gamma$-point, powder averaged neutron scattering spectra, elastic constants, and thermodynamics properties could be determined [8–13]. A measurement of the full dispersion relation and the investigation of the vibration contributing most to the first peak in the VDOS is not reported. Phonon dispersion relations were calculated using a Born–von-Karman force model [14], while first principles studies are limited to the $\Gamma$-point [15].

In this study we aim at a detailed examination of the lattice dynamics of $\alpha$-cristobalite and a connection to the Boson peak in ambient silica glass. Dynamical matrices and derived properties are calculated from first principles and compared to experiment. The distribution of low energy features over reciprocal space is analyzed by x-ray diffuse scattering. The scattering intensities are confronted to calculated thermal diffuse scattering (TDS) and compared to predictions from a rigid unit model [16] and electron diffraction [17]. Selected regions in reciprocal space are further investigated using inelastic x-ray scattering (IXS), which is also used to measure the dispersion relation along selected directions and the VDOS. The vibration contributing most to the first peak in the VDOS is identified with help of the validated calculation and compared to similar features in $\alpha$-quartz and coesite. Finally, the observed displacement patterns are put into relation with atomic motions reported for vitreous silica [18] to make the connection to the Boson peak.

The natural $\alpha$-cristobalite crystals from the Ellora Caves, Hyderabad, India used in this study were kindly made available from the Harvard Mineralogical Museum collection. An octahedral single crystal with 0.2 mm along the four-fold axis was used for both diffuse and inelastic x-ray scattering studies, which were conducted at room temperature. One major twin and some mosaic spread were observed during the x-ray diffuse scattering. One domain could be isolated for the experimental study thanks to the small x-ray beam size, and only a small contribution from the other twins was observed. A synthetic polycrystalline sample (Materials Research Institute and School of Physics and Astronomy, Queen Mary University of London, UK) was used for the powder IXS measurements. The purity of the polycrystalline sample was verified by high-resolution x-ray diffraction on beamline ID28 at the ESRF—The European Synchrotron. The samples revealed pure single-phase patterns. The crystalline quality was checked by x-ray diffuse scattering. Clear patterns of diffraction rings without noticeable effects of structural disorder were observed.

The x-ray diffuse scattering experiment was conducted on beamline ID29 [19] at the ESRF. Monochromatic x-rays with 17.712 keV energy were scattered from the crystal in transmission geometry. The sample was rotated with an increment of 0.1$^{\circ}$ orthogonal to the beam direction over an angular range of 360$^{\circ}$, while diffuse scattering patterns were recorded in shutterless mode with a PILATUS 6M detector [20]. The orientation matrix and geometry of the experiment were refined using the CrysAlis software package. 2D reconstructions were prepared using locally developed software.

The single crystal IXS study was carried out at ID28. The spectrometer was operated at 17.794 keV incident energy, providing an energy resolution of 3.0 meV full-width-half-maximum and a horizontal momentum resolution of 0.3 nm⁻¹. Energy transfer scans were performed at constant momentum transfer (Q) in transmission geometry along selected directions in reciprocal space. Further details of the experimental setup and the data treatment can be found elsewhere [21].

The powder IXS study was also conducted at ID28. IXS spectra were taken at momentum transfers chosen away from the Debye–Scherrer rings. The Q-range covered 10–70 nm⁻¹. The data combine measurements with 1.4 meV resolution at 23.725 keV incident energy within  −25 to 25 meV and 0.2 meV energy steps and measurement with 3.0 meV resolution at 17.794 keV incident energy within  −25 to 180 meV and 0.7 meV steps. The elastic peak in the IXS spectra was subtracted using the instrumental resolution function determined from a polymethylmethacrylate (PMMA) sample at momentum transfers close to the maximum of its structure factor. The generalized x-ray weighted VDOS (X-VDOS) was obtained by summing the IXS spectra within the incoherent approximation following a previously described data treatment procedure [22].

The lattice dynamics calculations were performed using density functional perturbation theory (DFPT) [23] as implemented in the CASTEP code [24, 25]. We apply the local density approximation (LDA) with the exchange correlation functional by Perdew and Zunger [26] and use the plane-wave formalism with norm-conserving pseudopotentials of the optimized form [27]. Three reference orbitals were treated as valence states for silicon and two for oxygen. The plane wave cut-off and the sampling of the electronic grid were carefully tested by evaluating the convergence of internal forces. The electronic structure was computed on a $6\times 6\times 6$ Monkhorst–Pack grid and the plane wave cut-off was set to 800 eV. The chosen settings provided a convergence of internal forces to  <10⁻³ eV Å. The cell geometry was optimized employing the Broyden–Fletcher–Goldfarb–Shannon method [28] by varying the lattice and internal parameters. For the cell parameters of the optimized cell we find a   =   b   =   4.939 Å and c   =   6.870 Å. The obtained values agree within 1.1% to those determined by x-ray diffraction (a   =   b   =   4.97(8) Å and c   =   6.94(8) Å [2]). Phonon frequencies and eigenvectors were computed on a $7\times 7\times 7$ Monkhorst–Pack grid of the irreducible part of the Brillouin zone by perturbation calculations. Sum rules for the acoustic branches as well as the charge neutrality at the level of Born effective charges were imposed with a maximum correction of 1.2 meV at $\Gamma$. A Fourier interpolation with a grid spacing of 0.005 Å⁻¹ in the cumulant scheme including all image force constants was applied for the VDOS [29, 30]. The calculation was tested to be well converged with a maximum error in phonon energies of  <0.05 meV.

Scattering intensities for inelastic and thermal diffuse scattering were calculated in first order approximation assuming the validity of both harmonic and adiabatic approximation. The formalism can be found elsewhere [21, 31]. The momentum transfer dependence of the atomic scattering factors was taken into account using an analytic function with coefficients derived from Hartree–Fock wavefunctions [32]. IXS spectra are resolved in phonon energy and the intensity is sensitive to the eigenvectors. TDS intensities are not resolved in energies, but due to their strong energy dependence, they are sensitive to low energy features.

High symmetry reciprocal space sections of experimental diffuse scattering and calculated TDS intensity distributions are shown in figure 1. A complex distribution of diffuse scattering is noticeable. The experimental intensity distributions show signatures of the mosaic spread, which appear as arclike features and small peaks aligned on Debye–Scherrer rings. Strong diffuse scattering of acoustic phonons is noted around the Bragg peaks and streaks of diffuse intensities are visible in $<1\,0\,0>$ and $<1\,1\,0>$ directions in the HK0 plane. The H0L plane shows diffuse scattering along $<1\,0\,1>$ which is more intense for low energy transfers, in particular between ($\bar{2}$ 0 0) and (0 0 2). Its origin was proven by IXS to contain a strong elastic contribution which corresponds to the inter-twin boundary and an inelastic contribution. In the HHL plane we note diffuse features along $<1\,1\,0>$ and even stronger streaks are observed along $<1\,1\,2>$. The diffuse scattering along $<1\,1\,0>$ and $<1\,0\,1>$ are in agreement with the prediction of the rigid unit modes model [16]. In this model the lines of diffuse features correspond to the zero-frequency solution to the dynamical equations for an infinite framework of rigid corner-linked tetrahedra. Diffuse scattering along these directions were also observed by electron diffraction [17]. The TDS intensity distribution shows features arising from low energy phonons. The slight deformation of the tetrahedra results in positive frequency solutions. The overall agreement between experiment and ab initio calculation is good and shows that the eigenvectors of the low energy phonons are well reproduced across reciprocal space. The calculated TDS intensity distribution in 3D around (3 3 2) is shown in figure 1. Around this Bragg peak the low energy features are particularly pronounced in all directions. Their distribution in 3D reciprocal space partially agrees with the prediction from the rigid unit mode model but is certainly more complex and the shape differs significantly, in particular along $<1\,1\,0>$. This may be appreciated by comparing figure 1 with the distribution of rigid unit modes wavevectors (figure 12 in [16]).

Zoom In Zoom Out Reset image size

The streaks of diffuse scattering along $<1\,1\,0>$ and $<1\,1\,2>$ and some more directions were explored by IXS. A few prototypical scans are shown in figure 2 and compared to calculated spectra. Some more IXS scans are summarized in the intensity maps along certain directions in figure 3. The spectra are reasonably well reproduced by the calculation. Applying a uniform scaling factor of 1.039 to the calculated phonon energies both energies and intensity of the phonon excitations match. The scaling factor was determined from the X-VDOS, and its value is justified further below. The minima of the lowest energy phonon branch are slightly underestimated in the showcases. A small elastic line is present in all spectra and not reproduced by the calculation. It arises from crystal defects. Figure 3a and b show the IXS intensity map along $<1\,1\,0>$ and $<1\,1\,2>$, respectively. We note quite flat phonon branches with energies around 5 meV which give rise to the diffuse scattering. The calculation underestimates the energies of the lowest energy phonon bands slightly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Calculated dispersion relations along high symmetry directions together with the experimental values from the IXS measurements are shown in figure 4. The phonon energies at the $\Gamma$ point are compared to infrared [12] and Raman measurements [10, 11] as well as to values calculated ab initio in a previous work [15]. We note good agreement of our calculation with all the experimental values. On average the calculated values are less then 0.3 meV off, and in the worst case 0.9 meV. The values from the previous calculation are significantly lower. We have performed a $\Gamma$-point calculation within the general gradient approximation in Perdew–Becke–Ernzerhof parameterization, as applied in the previous work rather then LDA and used the convergence criteria of our work. The resulting phonon energies are very similar to the previous result. We can thus conclude that both calculations are well converged and can attribute the discrepancy to the exchange correlation functional.

Zoom In Zoom Out Reset image size

The X-VDOS for $\alpha$-cristobalite and the calculated real VDOS are shown in figure 5. A linear scaling of all the calculated phonon energies by 1.039 leads to an almost perfect agreement of experimental and theoretical X-VDOS over the complete energy range. The applied scaling factor is very similar to what is found for the tetrahedral coordinated silica polymorphs $\alpha$-quartz, coesite and stishovite [33–35]. We observe some very small discrepancies between experiment and theory, which are mainly due to the limited accuracy of sampling the reciprocal space with powder IXS spectra. The first peak in the VDOS is observed at a slightly higher energy as calculated (the discrepancy is 0.9 meV after the applied scaling). The partial density of states (figure 5(b)) separate the contribution of silicon and oxygen atoms. Focusing on the low energy part of the partial VDOS (figure 5(c)) we find that the first peak located at 4.56 meV is slightly dominated by the vibration of the oxygen atoms. The low energy VDOS of $\alpha$-quartz is plotted in figure 5(d) for comparison. We note that the first peak of the VDOS of $\alpha$-cristobalite is located at a significantly lower energy than in $\alpha$-quartz. Furthermore, in $\alpha$-cristobalite the partial contribution of the oxygen atoms is slightly higher than in $\alpha$-quartz.

Zoom In Zoom Out Reset image size

The localization of the critical points contributing most to the first peak of the VDOS was determined by the simultaneous application of two filters. An energy filter of $\Delta E=0.3$ meV was applied to the ab initio calculated phonon energies of the first Brillouin zone, and $1/\vert {{\nabla}_{\boldsymbol{q}}}E\left(\boldsymbol{q}\right)\vert$ was computed within this energy window. The investigation reveals two saddle points with the same phonon energy contributing to the first peak in the VDOS of $\alpha$-cristobalite. The computation of the local contribution within a cube in the reciprocal space of $\Delta q=1$ nm⁻¹ shows that the first peak arises mainly from the region around the M point (1/2 1/2 0). The local contribution is shown in figure 6(a). The phonon dispersion surface containing the critical point shows a double degeneracy along $<1\,0\,0>$, $<0\,1\,0>$ and $<0\,0\,1>$, which is split along the $<1\,1\,0>$ direction, see figure 6(b). The lower sheet forms a saddle point, whereas the upper sheet forms a minimum which is of parabolic nature close to M. The peak in the VDOS arises thus from the lower branch for topological reasons [36]. The shape of the energy surface projected on the HK0 and HHL plane is shown in figures 6(c) and (d). A second saddle point with the same phonon energy located at (0.47 0.2 0.33) was found. Its contribution to the first peak in the VDOS is about a factor of two smaller compared to the contribution from the M point and less sharp.

Zoom In Zoom Out Reset image size

The VDOS of $\alpha$-cristobalite and the displacement pattern contributing most to its first peak is compared in figure 7 with the silica polymorphs $\alpha$-quartz and coesite [34]. We note, that the low energy part of the VDOS is different. In particular, the first peak in the VDOS is located at different energies. The critical points which are responsible for the first peak are located at the zone boundary for all three silica polymorphs under investigation. The shape of the energy surface in the vicinity of the critical points is different, but despite this, the associated displacement patterns are very similar: The largest displacement is observed for oxygen. The vibration consists mainly of a tetrahedron tilt, accompanied by a small distortion.

Zoom In Zoom Out Reset image size

The VDOS, obtained from the powder IXS spectra, probes the ensemble of vibrational states and was used to determine an overall scaling factor for the ab initio calculation. The application of a linear scaling factor deviating from unity by only 3.9% leads to a very good agreement between ab initio the calculated and measured phonon energies. The underestimation of the calculated energies can be attributed to the limited accuracy of the exchange correlation function within the local density approximation. Detailed discussions can be found elsewhere [25, 37].

By comparing all three, TDS intensities, single crystal, and powder IXS spectra over an extended region in reciprocal space, the calculation can be validated not only in terms of phonon energies but also in terms of phonon eigenvectors. The accurate determination of the complete lattice dynamics becomes possible. The full description of the lattice dynamics at arbitrary momentum transfer in turn allows for an investigation of the phonon energy landscape. The critical points responsible for Van Hove singularities can be localized by the simultaneous application of an energy and gradient filter to the calculated phonon frequencies. The eigenvectors and the shape in energy-momentum space of the vibration contributing most to the first peak in the VDOS can be assigned and compared to the tetrahedrally coordinated silica polymorphs $\alpha$-quartz and coesite. Although the first peak in the VDOS is located at different energies and the energy landscape in the vicinity of the critical points differs, a unifying picture on its origin can be drawn: The critical points are located at the zone boundary. The vibration contributing most to this first peak is mainly driven by oxygen vibrations and consists of a tetrahedron tilt, accompanied by a small distortion. Such librations of nearly undistorted SiO₄ tetrahedra were reported to cause the Boson peak in vitreous silica [18]. The similarity of the displacement patterns of $\alpha$-cristobalite and vitreous silica according to the Buchenau model of coupled rotations can be appreciated in tetrahedron representation, see figure 8.

Zoom In Zoom Out Reset image size

Comparing the VDOS of ambient silica glass and $\alpha$-cristobalite reveals that the number of excess states above the Debye levels is very similar (6.2(9)% and 6.8(9)%, respectively, [7]). For both systems, the first peak in the VDOS is located nearly at the same energy. The same holds true for densified glass and $\alpha$-quartz, if the densities are matched. Taking into account the similarity observed in the displacement patterns as well as in the VDOS of vitreous silica and crystalline polymorphs with matched densities, we can now conclude, that the Boson peak in silica glass is the counterpart of the van Hove singularities in silica polymorphs. The vibrational states causing this feature are similar in terms of absolute numbers, energies, and displacement pattern.

In summary, we have seen that the lattice dynamics of $\alpha$-cristobalite can be very well described by DFPT in LDA. The thermal diffuse scattering shows a rich structure. Some, but not all, of the strong features are in good agreement with the simplified picture of rigid unit modes. A full description of phonon frequencies and eigenvectors could be retrieved and the energy landscape of specific features in the lattice dynamics became possible. The origin of the first peak in the VDOS could be identified and linked to the Boson peak in vitreous silica.

In the framework of the silica polymorphs we are now able to provide accurate models of the lattice dynamics of $\alpha$-quartz, coesite, $\alpha$-cristobalite, and stishovite [33–35]. The calculation for these polymorphs can be extended to high pressures and allow the derivation of elastic and thermodynamical properties. Such an extension potentially builds the basis for understanding the compression mechanism and phase stability at geophysically relevant conditions. The prediction of the lattice dynamics should be possible within the employed calculation scheme for the high pressure phases of cristobalite [38] as well as for the recently discovered high-pressure silica polymorphs seifertite and others [39]. For cristobalite the knowledge of the lattice dynamics at high pressures is in particular important for understanding the pressure evolution of its auxetic behavior [40], as well as the phase stability and elasticity of its high pressure phases.

The authors would like to acknowledge the courtesy of the Harvard Mineralogical museum and personally the collection curator Dr R Alonso-Perez for providing the high quality pieces of natural cristobalite crystals (Cristobalite, HMM#97849) and Martin Dove (Materials Research Institute and School of Physics and Astronomy, Queen Mary University of London, UK) for the synthetic polycrystalline sample.