Introduction

Mineral inclusions entrapped in ultrahigh-pressure metamorphic rocks can provide fundamental information about geological processes such as subduction and continental collision. When a host-inclusion pair is exhumed from depth to the Earth’s surface non-lithostatic stresses are developed in the inclusion because of the contrast in their elastic properties (Angel et al. 2015 and references therein). The inclusion is under compressive stress if it is more compressible than the host. If correctly interpreted, these stresses on the inclusion allow the stress conditions at entrapment to be estimated. However, the theory for elastic geobarometry for host-inclusion pairs is only developed for crystals immersed in a hydrostatic stress field. This is valid for an isotropic and spherical inclusion entrapped in an isotropic host as it is subject to isotropic strains imposed by the host and will, therefore, exhibit isotropic stresses; that means the inclusion will be under hydrostatic pressure. However, if inclusions of elastically anisotropic minerals such as quartz, coesite, or zircon are entrapped in a cubic host such as garnet the same isotropic strain imposed by the garnet will result in anisotropic (i.e., deviatoric) stresses in the inclusion (see Fig. 1). Thus, the stress state of the inclusion will be different from the hydrostatic case, and it cannot be characterized by a single “pressure” value (Anzolini et al. 2018). Furthermore, the shifts of Raman mode frequencies under deviatoric stress are unknown in general. The limited experimental evidence (e.g., Briggs and Ramdas 1977) is that the Raman modes change differently from those measured under hydrostatic pressure. Therefore, the key question remains, posed by Korsakov et al. (2010), can shifts of Raman modes measured in hydrostatic experiments be used to interpret the Raman shifts from an inclusion under non-hydrostatic stress? And what errors does this introduce into estimates of inclusion stress and entrapment conditions?

Figure 1

Deviatoric stress |σ₁ – σ₃| in a quartz inclusion in a pyrope garnet (where σ₁ and σ₃ are parallel to the a and c axes in quartz, respectively) at room conditions as a function of P (GPa) and T (°C) of entrapment. The EoS of quartz and pyrope are from Angel et al. (2017) and Milani et al. (2015). The anisotropic relaxation was calculated using the elastic tensors for quartz and pyrope (Sinogeikin and Bass 2002; Lakshtanov et al. 2007).

Since it is challenging to perform experiments under controlled deviatoric stress conditions, we have performed ab initio Hartree-Fock/Density Functional Theory (HF/DFT) calculations to determine how the Raman modes of quartz change under both hydrostatic pressure and deviatoric stress conditions. We show how the HF/DFT simulations can be used to determine the strains (and by inference the stresses) within a crystal inclusion via mode Grüneisen tensors and measurement of its Raman shifts.

Grüneisen tensor and strains

It is well known that when stresses are applied to a crystal, for example by changing the pressure, the phonon wavenumbers are shifted. This is most easily seen in the change in Raman peak positions measured under hydrostatic pressure. The shifts of certain modes of materials (e.g., the 464 cm^–1 mode of quartz) have frequently been used as secondary pressure standards. However, there is a common misconception that the observed wavenumber shifts of the Raman-active modes are directly related to the magnitude of the applied pressure or stress. If that were true, then Raman modes would not show a change in wavenumber when the temperature of a free crystal is changed at ambient pressure. Instead, Raman modes generally exhibit a decrease in wavenumber as the temperature is increased. As Figure 2 shows, the shift of the 464 cm^–1 mode of quartz with both pressure and temperature shows the same dependence upon volume. These observations point to the correct interpretation that the shift Δω of the wavenumber ω of a vibrational mode is primarily due to the strains of the crystal induced by the applied temperature or pressure. Thus, the Raman peak position of a crystal under a strain ɛ (i.e., the full strain tensor) is determined by the second-rank symmetric tensor: the mode Grüneisen tensor γ^m (Ziman 1960; Key 1967; Cantrell 1980) characteristic of each phonon mode m, which can be written in Voigt (1910) notation as:

Figure 2

The wavenumber shift (Δω) of the 464 cm^–1 mode of a free quartz crystal measured under different P, T, and stress conditions form a single trend with volume strain ɛ_V. The deviation at positive strains is caused by the transition from α to β quartz that occurs at 573 °C. Away from the transition, the experimental data are reproduced by both the HF/DFT simulations under hydrostatic pressure (red filled squares), and the prediction from the mode Grüneisen parameters (line).

This equation means that the changes in the Raman peak positions, in general, depend on all the components of the strain tensor (see Supplemental material^[1]), not just on the relative change in the volume (i.e., ɛ₁ + ɛ₂ + ɛ₃). Because the mode Grüneisen tensor is a symmetric second-rank property tensor, it is subject to the same symmetry constraints on its component values as other second-rank property tensors, such as the thermal expansion and compressibility tensors. For the trigonal symmetry of quartz γ1m=γ2m≠γ3m and γ4m=γ5m=γ6m=0. Garnet is almost elastically isotropic (Sinogeikin and Bass 2002) so the strain imposed by the garnet host does not break the symmetry of the quartz inclusion. Therefore, we are interested specifically in the cases when ɛ₁ = ɛ₂ for which the shift in the phonon wavenumbers should be given by:

To determine the values of γ1m and γ3m independently of one another ¹³it is not sufficient to measure the Raman shifts under hydrostatic pressure, because only a single series of values of ɛ₁ and ɛ₃ are measured. We, therefore, use HF/DFT simulations to calculate the Raman modes at different imposed strains on the crystal to determine its structure and properties at those strains.

Raman shifts and strains from ab initio calculations

Ab initio HF/DFT simulations have been performed by means of CRYSTAL 14 code (Dovesi et al. 2014) employing the hybrid Hamiltonian WC1LYP (Wu and Cohen 2006) that has been demonstrated to be particularly suitable for the accurate reproduction of the elastic and vibrational properties (Zicovich-Wilson et al. 2004; Prencipe et al. 2014). Further computational details are reported in the Supplemental material^[1].

We will only discuss in detail the Raman-active modes near 464 cm^–1 (non-degenerate A₁ mode) and 696 cm^–1 (doubly degenerate E mode) because they give peaks in the Raman spectra that are easily resolved from the peaks of the host garnets. Note that all E modes in quartz are polar, and therefore have longitudinal optical (LO) and transverse optical (TO) components that generate two separate Raman peaks, whose intensity ratio varies depending on the scattering geometry. Because of the polarization mixing (Loudon 1964), the wavenumber of the LO component depends on the angle between the triad axis of quartz and the phonon-propagation direction that in the case of backscattering geometry coincides with the direction of the laser beam. For the E mode of interest in this study, the LO-TO splitting is small. At ambient conditions w_ETO = 696 cm^–1, while the maximum w_ELO ~ 697.5 cm^–1 and it corresponds to the case when the c axis is perpendicular to the laser beam. The two components are thus very close to each other and in specific experimental geometries may generate peaks of similar intensity. To avoid possible contribution from the LO component that may lead to a subtle artificial shift of the corresponding Raman peak toward higher wavenumbers, one should either rotate the sample about the direction of the laser beam to verify that the peak position does not change, or find the orientation at which the wavenumber of the Raman peak is lowest (which should be when the c axis is perpendicular to the polarization of the incident light).

Figures 3a and 3b are contour maps that display the HF/DFT wavenumber shifts of the A1 mode near 464 cm^–1 and ETO mode near 696 cm^–1 as a function of the two independent strain components. For small strains (close to the origin) the iso-shift lines are parallel to one another and equally spaced, confirming that the values of γ₁^m and γ₃^m are constants over these strain ranges. The contour lines for the two modes have different slopes (e.g., Figs. 3a, 3b, and 3c) indicating that the values of their Grüneisen components are different. Furthermore, Figure 3a shows that the iso-shift lines are not parallel to isochors, which are represented by lines of equal volume strain [ɛV = 2ɛ₁ + ɛ₃]. If the wavenumber shifts are plotted against the stress components σ₁ and σ₃ (Fig. S4, Supplemental^[1] material) then the iso-shift lines are not parallel to isobars [lines of constant pressure (2σ₁ + σ₃)/3 = –P]. Therefore, in general, Raman shifts do not measure either volume or pressures; for uniaxial crystals like quartz, they indicate the principal normal strain components ɛ₁ and ɛ₃. The Grüneisen tensor components for all modes were then determined by fitting Equation 2 by least-squares to the wavenumber shifts at different strain states simulated by HF/DFT calculations. For the 464 mode, the maximum misfit was 1.76 cm^–1 (see Supplemental material^[1] for residual plots and full details).

Figure 3

(a, b, and c) Wave number shifts Δω (cm^–1) of quartz calculated by HF/DFT as a function of the two independent strain components (ɛ₁ = ɛ₂ ≠ ɛ₃). The reference value at zero strain is the ab initio datum at 0 K and 0 GPa (static pressure). The symbols indicate strains at which HF/DFT simulations were performed. The colored bands are the iso-shift lines (lines with the same frequency shift for a specific Raman mode). The iso-shift lines are approximately parallel to one another and equally spaced. The iso-shift lines are not parallel to isochors defined by constant volume strain ɛ_V = 2ɛ₁ + ɛ₃.

HF/DFT calculations are performed under static stress at absolute-zero temperature. Inclusions are measured at room temperature. To use the mode Grüneisen tensors to determine strains in inclusions, one has to also demonstrate that their values are independent of P and T. This is equivalent to the solid behaving according to the quasiharmonic approximation (QHA). We achieve this by calculating the strains of a free quartz crystal relative to room conditions from the known unit-cell parameter variation of quartz (Angel et al. 2017) with P and T. We then compare the wavenumber shifts from room-condition values using the Grüneisen components (γ1m and γ3m) for the two modes at 464 and 696 cm^–1 (i.e., 0.60 and 1.19; 0.50 and 0.36) determined by HF/DFT and Equation 2 with experimental data. The line in Figure 2 shows that for the mode near 464 cm^–1 the experimentally measured wavenumber shifts under high pressure at room T and those up to ~400 °C at room P are reproduced by the Grüneisen components determined by HF/DFT. At higher temperatures, the predicted shifts differ from the experimental data because of pre-transition effects associated with the α-β quartz transition that cannot be accounted using QHA. Furthermore, the fundamental soft mode near 206 cm^–1, which is heavily involved in the temperature α- to β-quartz phase transition (Scott 1968), clearly violates the QHA, and HF/DFT simulations cannot be used to determine the mode Grüneisen components. To use this band to determine strains in inclusions, the variation of its position with strain must be determined experimentally.

Validation from Raman scattering and X‑ray diffraction

To validate this approach, we performed micro-Raman spectroscopy and X‑ray diffraction measurements on quartz inclusion in garnet from a diamond-grade eclogite xenolith (TM90-1) from the Mir kimberlite pipe (Yakutiya). We selected this example as it has the highest Raman shifts reported for quartz inclusions in garnet (Korsakov et al. 2009; Zhukov and Korsakov 2015). Parallel polarized Raman spectra were collected in backscattering geometry with a Horiba Jobin-Yvon T64000 triple-monochromator spectrometer (spectral resolution of ~2 cm^–1, instrumental accuracy in peak positions of ~0.35 cm^–1 and 2 μm spot size) following the same protocol reported in Campomenosi et al. (in review). Figure 4 shows that the shift of the 464 cm^–1 Raman line changes significantly across the crystal as a consequence of edge and corner effects (Campomenosi et al. in review; Mazzucchelli et al. 2018), and is lowest at the center of the inclusion. We determined the unit-cell parameters of the inclusion by single-crystal X‑ray diffraction measurements using the 8-position centering method (SINGLE, Angel and Finger 2011) using a newly developed Huber 4-circle Eulerian cradle diffractometer equipped with point detector and microfocus source (120 μm spot size). The difference between the unit-cell parameters of the inclusion a = 4.86669(44) Å, c = 5.35408(14) Å, and those of a free crystal measured on the same instrument [a = 4.91160(7) Å, c = 5.40325(9) Å] shows that the inclusion is under strains ɛ₁ = –0.00914(9), ɛ₃ = –0.00910(26), and ɛ_V = –0.02714(30). Relative to room conditions the quartz inclusion is currently under isotropic strain within the experimental uncertainties. From these strains and the Grüneisen-tensor components for the mode near 464 cm^–1 (i.e., γ1464=0.64 and γ3464=1.19), we calculate an expected wavenumber shift of 10.12 cm^–1. This is between the minimum and maximum shifts actually measured on the inclusion (Fig. 4) because the unit-cell parameters measured by X‑ray diffraction are an average over the entire volume of the inclusion. Conversely, if we take the shifts of the Raman modes near 464 and 696 cm^–1 measured at the center of the inclusion and the calculated mode Grüneisen components we predict the following strains: ɛ₁ = –0.0093(5) and ɛ₃ = –0.0070(5). In this case, the total strain is lower than that measured by X‑ray diffraction, again because this method provides an average over the inclusion.

Figure 4

Wavenumber shifts of the 464 and 696 cm^–1 modes measured in a traverse across a quartz inclusion in pyrope from the eclogite xenolith TM90-1 (Korsakov et al. 2009). The phonon wavenumber shifts calculated from the strains determined by X‑ray diffraction and the mode Grüneisen tensor are in good agreement with the average of the measured shifts.

Implications

We have demonstrated by a combination of HF/DFT simulations and comparison to experimental Raman data that quartz has some Raman-active modes whose wavenumbers are only a function of strain (for small strains) and are not directly dependent on P and T. The HF/DFT simulations show also that in general the shift of a Raman line does not indicate either the volume change (or volume strain) or the mean normal stress but is a more complex function of the linear strains (Fig. 3). The linear strains of a quartz inclusion can be determined from the wavenumber shifts of at least two Raman peaks by using the mode Grüneisen components determined by HF/DFT. Conversely, a wavenumber shift for a specific Raman-active mode can be determined from the linear strains measured by X‑ray diffraction and by using the mode Grüneisen components determined by HF/DFT. We proved that strains and Raman shifts vary significantly across an inclusion (Fig. 4) as a result of the influence of shape combined with elastic anisotropy (Mazzucchelli et al. 2018; Campomenosi et al. in review), and should never be averaged. The Raman shifts measured in the center of inclusions are least affected by the presence of edges and should be the only ones used to infer, after shape corrections (Mazzucchelli et al. 2018), the entrapment conditions. On the other hand, X‑ray diffraction measurements give strain components that are averaged over the entire volume of the investigated crystal. For the example discussed, the interpretation of wavenumber shift of the mode near 464 cm^–1 measured at the center of the quartz inclusion with the hydrostatic calculation (Schmidt and Ziemann 2000) leads to only a small error in the estimation of the mean normal stress (i.e., pressure) of ca. 0.03 GPa. This would lead to an error in the calculated entrapment pressure of 0.06 GPa. However, this is not always the case. When the deviatoric stress (Fig. 1) in the inclusion is higher, the error is larger. For example, a quartz inclusion in garnet reset at 1.5 GPa and 1100 °C during exhumation will have a wavenumber shift of the Raman mode near 464 cm^–1 that would yield an entrapment pressure of 0.6 GPa too low when the hydrostatic calibration is used. In anisotropic hosts, the discrepancy from the hydrostatic calibration also depends on the relative orientation of the host and inclusion. For example, for quartz in zircon the error in the calculated entrapment pressure can vary from 0.2 to 0.6 GPa for the same entrapment condition (i.e., 0.3 GPa and 800 °C) depending on orientation.