Quartz under stress: Raman calibration and applications of metamorphic inclusions to geobarometry

An experimental calibration of the shifts of three major Raman peaks of quartz with hydrostatic pressure and uniaxial differential stress is presented, and implications for their use in geobarometry based on Raman spectroscopy of quartz inclusions are discussed. The position of the 206 cm⁻¹ peak depends only on hydrostatic pressure P, and its pressure dependence is recalibrated with a peak-fitting procedure that is more adequate for Raman barometry than previous calibrations. The position of the 128 and 464 cm⁻¹ peaks depends on P and also on differential stress σ, which can be determined from the position of these two peaks knowing hydrostatic pressure from the position of the 206 cm⁻¹ peak. The results obtained here are different from those inferred previously from first-principles calculations. The present calibration provides direct relationships between Raman shifts and stress, with a simple formulation of residual pressure and differential stress assuming uniaxial stress along the c axis of quartz inclusions. It is tested on data from experimental and natural inclusions. Residual pressures from the present calibration are similar within uncertainties to those obtained with previous experimental calibrations. Residual differential stresses obtained from the 128 and 464 cm⁻¹ peaks are very sensitive to the precision of Raman measurements. Experimental inclusions yield residual pressures consistent with synthesis pressure. Differential stresses obtained on some experimental inclusions are sometimes incompatible, providing a criterion for identifying inclusions under complex stress conditions that are not appropriate for geobarometry. Recent data on natural inclusions show self-consistent differential stress, consistent with the assumption of major stress along the symmetry axis of the inclusion crystals. The average pressure values from the 128 and 464 cm⁻¹ peaks are similar to the residual pressure from the 206 cm⁻¹ peak that depends only on hydrostatic pressure. It can be used to obtain pressure when the 206 cm⁻¹ peak position cannot be used due to interference with host mineral peaks. Using the 128 and 464 cm⁻¹ peaks alone, or averaging either 128 and 206 or 206 and 464 cm⁻¹ peaks, can induce systematic bias in the residual pressure determination. Applications of the present results to natural inclusions suggest that combined determination of residual pressure and differential stress may be used for both barometry and thermometry pending further calibration.

1 Introduction

Raman spectroscopy has been used for over two decades for determining residual pressures in mineral inclusions in diamonds (Izreali et al., 1999) and metamorphic minerals, principally garnets (Enami et al., 2007; Parkinson and Katayama, 1999). Entrapment pressures calculated from thermoelastic modeling provide a determination of metamorphic pressures (Angel et al., 2017b; Kohn, 2014; Zhong et al., 2019; Zhang, 1998) that is independent of those obtained from metamorphic phase equilibria through thermodynamic modeling (Connolly, 1990; Holland and Powell, 1998). The method has been validated with experimental entrapment of quartz in garnet at controlled hydrostatic pressures (Bonazzi et al., 2019; Thomas and Spear, 2018).

Quartz is often used because it is a common mineral of almost pure composition, hence the position of its Raman peaks depends on intensive parameters and not on chemical variations. Accurate determinations of residual pressures rely among others on precise determinations of the Raman peak shifts with pressure that are provided by diamond anvil cell (DAC) experiments (Schmidt and Ziemann, 2000). Differential stresses also influence the calculations of residual pressures, and determinations of strains from Raman frequency of quartz inclusions have relied on first-principles calculations of Raman peak shifts (Murri et al., 2019), a method that was also successfully used for determining differential stresses in high-pressure experiments (Reynard et al., 2019) but that is backed only by experimental data at liquid helium temperature (Briggs and Ramdas, 1977; Tekippe et al., 1973).

Figure 1 Typical spectra of quartz obtained in the back-scattering geometry at ambient conditions (0 N) and under uniaxial force of 2400 N (corresponding to a stress of ∼0.6 GPa) perpendicular to the c axis that appears vertical under the microscope. Spectra were measured with incident polarization parallel (0^∘) and perpendicular (90^∘) to the c axis to better characterize the splitting of the 128 cm⁻¹ peak under stress. No polarizer was put along the scattered light path. Asterisks mark the calibration Ne and Hg lines near 269.5 and 473.5 cm⁻¹.

New experiments were performed to directly calibrate the Raman peak shifts induced by hydrostatic and non-hydrostatic stresses. Consequences for barometric applications are discussed based on available data on synthetic inclusions from piston–cylinder experiments (Bonazzi et al., 2019) and from recent datasets on natural quartz inclusions in different suites of metamorphic rocks (Cisneros et al., 2021; Gonzalez et al., 2019; Zhong et al., 2019).

2 Experiments

Raman spectra under stress were obtained in the back-scattering geometry using a Horiba™ HR Evolution (532.31 nm excitation from a diode-pumped solid-state laser) and 1800 g mm⁻¹ at LGL in ENS de Lyon. Only the three most intense Raman peaks at 128, 206, and 464 cm⁻¹ are analyzed and will be referred to by these numbers in the remainder of the text. A series of experiments were first performed with calibration against silicon standard (crystals #1 to 4). In order to check for instrumental drift, sharp Hg 18312.55 or Ne 18516.59 cm⁻¹ lines (269.46 and 473.5 cm⁻¹ shift with respect to the excitation laser at 18786.05 cm⁻¹, Fig. 1) were recorded along with all spectra in a second series of experiments (crystal #5 and higher). The maximum drift of the reference lines' position during one experimental session is 0.5 cm⁻¹ and was systematically corrected. The measured separation between the two reference lines is 204.03 cm⁻¹ on average with a maximum deviation of 0.15 cm⁻¹. The splitting of the two components of the 128 cm⁻¹ peak under stress perpendicular to the c axis was determined by fitting two peaks with identical width. The incident polarization was rotated parallel and perpendicular to the c axis of the crystal face to enhance one component with respect to the other (Fig. 1). Quartz inclusions in eclogitic garnets from the Bergen Arcs (Zhong et al., 2019) were studied at the Freie Universität Berlin. Raman spectra were obtained with a Horiba LabRAM confocal Raman spectrometer using the 532 nm line of a Nd:YAG, 1800 g mm⁻¹ gratings, and Olympus 50× and 100× objectives. The laser energy was set at ca. 10 mW. The typical duration for each measurement was 1 to 2 min with >5 times repetition to obtain good signal-to-noise ratio. Immediately after each inclusion measurement, a gem-quality quartz crystal was measured as stress-free reference material to obtain the wavenumber shift of quartz inclusions due to residual stress.

Peak positions were fitted using symmetric Voigt profiles, although some bands, in particular the 206 cm⁻¹ peak, are slightly asymmetrical in shape. This ensures that the present calibration is comparable with data obtained on natural inclusions where a Voigt profile must be used because interference with bands of the host garnet does not allow us to use more complex functions to fit peak positions. Nominal frequencies of the three main quartz peaks of 127.1(2), 206.1(2), and 464.2(2) are obtained from internally calibrated spectra (Fig. 1). For the sake of comparison with previous studies (Bonazzi et al., 2019; Schmidt and Ziemann, 2000; Cisneros et al., 2021; Gonzalez et al., 2019; Zhong et al., 2019), all data were normalized to match the ambient condition values of 128, 206, and 464 cm⁻¹.

Raman spectra of quartz were measured under hydrostatic pressure at ambient temperature on randomly oriented small crystals (<10 µm) in a DAC in a hydrostatic methanol–ethanol–water mixture (Reynard et al., 2015, 2019) with small pressure steps up to 2 GPa. Pressure was measured with ruby fluorescence (Piermarini et al., 1975) for comparison with previous data (Schmidt and Ziemann, 2000). A Deben™ 5 kN compression cell was used to compress parallelepipeds of quartz that were cut along or perpendicular to the c axis at ambient temperature under uniaxial force, in a similar fashion to former experiments at 4 K (Tekippe et al., 1973). Crystals were placed between the two flat anvils of the system. A buffer made of an annealed 200 µm thick 301 stainless-steel foil was placed between the crystal faces and the anvils to lower the shear stress on the crystal ends through viscous flow of the metal. Apiezon grease was applied between the anvils, metal foils, and crystals to reduce shear stress. The viscous grease helped to maintain the crystal stuck to one anvil while the system was being closed with a force <10 N (pressure <2 MPa), sufficient to hold it in place. Cutting orientation was performed optically with respect to crystalline faces of the initial crystal. Orientation of final rod faces with respect to simple crystallographic orientations was determined with conoscopic observations. The deviation from crystallographic axes was determined on unused rods with Raman intensities (Zhong et al., 2021b) and found to lie within the 2^∘ range (Fig. S1 in the Supplement). The applied force is converted to stress by dividing with the cross section of the crystal rod, which was measured with an estimated uncertainty of 2 %. A supplementary data sheet provides results from successful experiments (4 out of 20) and one example of an unsuccessful experiment with early cracking of the crystal.

3 Elastic modeling

Interpretation of the Raman data on quartz inclusions relies on modeling of the residual pressure and differential stress for inclusions entrapped into a given host, generally garnets. The residual stresses are calculated based on an analytical model (Zhong et al., 2021a) derived assuming a pure almandine garnet host. The model treats anisotropic inclusion entrapped in an infinite and isotropic host using the classical Eshelby solution and the equivalent eigenstrain method (Eshelby, 1957; Mura, 1987). The unit cell parameters of quartz are fitted based on the X-ray data measured at ambient pressure and high temperature (Carpenter et al., 1998) and at ambient temperature and high pressure (Angel et al., 1997) with an equation of state (EoS) taking into account a curved alpha–beta transition (Angel et al., 2017a). The fitted results represent the axial EoS applicable to high P–T (pressure–temperature) conditions considering the alpha–beta transition. For the almandine host, we used the EoS from Milani et al. (2015). Entrapment stress is assumed to be hydrostatic. Due to the thermoelastic anisotropy of the inclusion, the residual stress becomes non-hydrostatic after cooling and decompression. Here, compressive stress is defined as negative. A similar approach has been used (Alvaro et al., 2019) to deal with different axial stress (strain) for quartz inclusion entrapped into an isotropic host.

Figure 2 Model residual pressure and residual differential stress on quartz included in almandine. Expected residual pressure and stress corresponding to conditions of entrapment of experimental and natural inclusions are shown for comparison with measured values (Figs. 5 and 6). Actual pressures for natural inclusions may be shifted by less than 0.05 GPa due to compositional effects on the garnet equation of state.

Results show that expected residual pressures are generally positive for most crustal metamorphic conditions and can be negative if the entrapment conditions are low pressures close and beyond the α–β quartz transition (Fig. 2). The residual differential stress is uniaxial with the symmetry axis of the stress parallel to the c axis of the quartz inclusion. Stress on the inclusion is hydrostatic (${\mathit{\sigma }}_c-{\mathit{\sigma }}_a=\mathrm{0}$) for a thermal gradient of ∼220 ^∘C GPa⁻¹ or ∼6 ^∘C km⁻¹. Stress along the c axis is higher than along the a axis for higher temperature and lower pressure conditions, and lower for lower temperatures and higher pressures that are uncommon in metamorphic rocks. Residual pressure is more sensitive to entrapment pressure and residual stress to entrapment temperature (Fig. 2). A combined determination of the two has potential geothermobarometric application (Alvaro et al., 2019).

Figure 3 Dependence of the 206 cm⁻¹ peak on (a) hydrostatic pressure at ambient temperature and (b) uniaxial stress. The small discrepancy between present hydrostatic compression data and former experimental data (Schmidt and Ziemann, 2000) is due to fitting with a different peak shape (symmetrical Voigt used here instead of asymmetrical Pearson IV). Uniaxial stress dependence obtained here is higher than that measured at 4 K (Briggs and Ramdas, 1977; Tekippe et al., 1973) due to temperature effects. The variations in peak position depend little on compression direction for the 206 cm⁻¹ peak in both studies, indicating that its position depends essentially on hydrostatic pressure. (c) Grüneisen plot for hydrostatic compression showing the constant Grüneisen parameter (γ) for the 464 cm⁻¹ peak and strong curvature for other peaks. Curves were fitted to the data with a volume dependence of γ expressed as the q parameter (see text and Table 1). For the sake of easy comparison, relative frequency shifts are multiplied by 2 and 4 for the 128 and 464 cm⁻¹ peaks, respectively. Pressure was converted to volume using the quartz equation of state (Angel et al., 1997).

4 Results

4.1 Hydrostatic compression

The hydrostatic pressure effects were calibrated in the range of 0–2 GPa that covers applications of Raman piezospectroscopy to metamorphic quartz inclusions (Fig. 3). For the 128 and 464 cm⁻¹ peaks, the present results are consistent with previous calibrations (Schmidt and Ziemann, 2000), with a maximum deviation of 0.4 cm⁻¹ and an average deviation of 0.1 cm⁻¹, i.e., within the inferred accuracy of the Raman spectroscopy when internal standards provided here by Ne or Hg emission lines are used.

Table 1 Pressure and uniaxial stress dependence of the three major Raman peaks of quartz.

^a First-principles results for trigonal strains from Murri et al. (2019) are shown in italics; ^b compression perpendicular to the c axis, the 128 cm⁻¹ E mode is split due to symmetry reduction to a monoclinic structure; ^c compression parallel to c axis; ^d q is fixed at zero for the 464 cm⁻¹ peak.

For the sake of comparison with former studies, data were fitted with the equation

$$\begin{array}{}\text{(1)}& P=A\mathrm{\Delta }{\mathit{\nu }}_h+B\mathrm{\Delta }{{\mathit{\nu }}_{\mathrm{h}}}^{\mathrm{2}},\end{array}$$

where Δν_h is the shift under hydrostatic stress with respect to the ambient conditions wavenumber (ν) of the Raman peak. The values of A and B are reported in Table 1 for the three main peaks of quartz. The Raman shift and its derivative are

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{h}}=[-A+(A^{\mathrm{2}}+\mathrm{4}BP{)}^{\mathrm{1}/\mathrm{2}}]/\mathrm{2}B,\text{(3)}& {\displaystyle }{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial P}}\right)}_{\mathit{\sigma }}=(A^{\mathrm{2}}+\mathrm{4}BP{)}^{-\mathrm{1}/\mathrm{2}}.\end{array}$$

It is worth noting that if the calibration is performed including data at pressures above 2 GPa, the initial slope is underestimated with respect to the present value, especially for the 206 cm⁻¹ peak. This is related to the complex pressure dependence of quartz Raman spectrum beyond the 2 GPa pressure range (Morana et al., 2020; Hemley, 1987). Using our own data above 2 GPa (Reynard et al., 2019), we noticed a similar systematic deviation. The present calibration with its small pressure steps (15 points within 2 GPa) is particularly adapted for applications to quartz inclusions that are under residual pressures below 2 GPa, a condition that applies even to those that underwent maximum pressure conditions near the stability field of coesite with maximum residual pressure around 1.3 GPa (Alvaro et al., 2019; Korsakov et al., 2009).

4.2 Uniaxial compression

The effects of uniaxial stress, applied either parallel or perpendicular to the c axis of quartz, differ strongly from peak to peak (Table 1). Stress effects are linear within uncertainties within the investigated range of stress up to ∼0.6 GPa beyond which crystals break. No symmetry breaking could be detected for compression along the c axis in spite of the 2^∘ deviation from the exact crystallographic direction (Fig. S1), hence analysis in the trigonal symmetry is valid. Symmetry breaking was observed for compression perpendicular to the c axis, marked by the splitting of the two components of the 128 cm⁻¹ E mode, in line with former experiments (Briggs and Ramdas, 1977; Tekippe et al., 1973).

For the 206 cm⁻¹ peak, the effects of stress are similar within uncertainties perpendicular and parallel to the c axis (Fig. 3c), with a stress dependence of 10.0(4) and 10.8(2) cm⁻¹ GPa⁻¹, respectively. Because the shift of the 206 cm⁻¹ peak is independent of the stress orientation within uncertainty, its position gives an absolute reading of the pressure P whether the stress is hydrostatic or not. Similar stress-induced shifts are obtained for crystal rods with different aspect ratios (1 and 4 along the c axis, and ∼2.5 and 4 along the a axis), ruling out potential geometrical effects (Table S1 in the Supplement). For the 128 and 464 cm⁻¹ peaks, the shift is very different when stress is applied perpendicular or parallel to the c axis. For stress applied along directions perpendicular to the c axis, the two components of the 128 cm⁻¹ E mode split due to the symmetry breaking (Tekippe et al., 1973). The splitting is maximum for uniaxial stress in the basal plane of quartz, a property that was used to measure differential stresses in non-hydrostatic experiments with DAC (Reynard et al., 2019). The shifts of the 128 and 464 cm⁻¹ peaks relative to the 206 cm⁻¹ peak are linear within the investigated stress range (Fig. 4).

Figure 4 Stress dependence of the (a) 464 and (b) 128 cm⁻¹ peaks relative to the 206 cm⁻¹ peak at ambient temperature. Black lines show the relative shift from experiments at ambient conditions (this study), and blue arrows the relative shifts for a uniaxial stress of 0.5 GPa at 4 K (Briggs and Ramdas, 1977; Tekippe et al., 1973). For the 128 cm⁻¹ peak, the two components of the E mode are split by compression perpendicular to the c axis (open circles and arrows). Stress dependence (long dashed line) is the average of that of the two components (Briggs and Ramdas, 1977; Tekippe et al., 1973). A few data points were collected at effective stresses higher than 0.6 GPa in partially broken crystals, and the stress is then not calibrated because the crystal section is not known anymore. These points plot along extrapolation of the linear trends even though they were not included in the fit.

As a check of the internal consistency of the measured shift, the sum of axial stress-induced shifts (i.e., twice the shift perpendicular to the c axis plus shift parallel to it) should compare with the slope of the bulk hydrostatic pressure dependence at zero pressure (Briggs and Ramdas, 1977). The sum of the axial stress dependences of the 206 cm⁻¹ peak is 30.8(10) cm⁻¹ GPa⁻¹, similar to the bulk pressure-induced shift of 30.7(13) cm⁻¹ GPa⁻¹ at ambient pressure, as well as for the 128 and 464 cm⁻¹ peaks, with values of 7.1 and 10.5 cm⁻¹ GPa⁻¹ for the sum of axial stress dependences, and 7.9 and 9.4 cm⁻¹ GPa⁻¹ for the hydrostatic pressure dependence, respectively (Fig. 3, Table 1). The consistency between the sum of the axial stress dependence and hydrostatic pressure experiment suggests that the effect of symmetry breaking in uniaxial compression on stress-induced shifts is of second-order magnitude.

5 Comparison with former studies

5.1 Hydrostatic compression

The present variation of the 206 cm⁻¹ peak shows an increasing discrepancy with that of Schmidt and Ziemann (2000) with increasing pressure, reaching about 1.5 cm⁻¹ at 1.5 GPa (Fig. 3a) above the accuracy of the measurement. This discrepancy is due to the use of a symmetric Voigt function in the present study while Schmidt and Ziemann (2000) used an asymmetric Pearson IV distribution. For piezospectroscopic applications, the Voigt function is widely used to fit the 206 cm⁻¹ peak position (Alvaro et al., 2019; Bonazzi et al., 2019; Cisneros et al., 2020, 2021; Enami et al., 2007; Gonzalez et al., 2019; Thomas and Spear, 2018; Zhong et al., 2019) because it avoids unrealistic fits due to interference with host mineral Raman peaks (garnets in most cases). Thus the present hydrostatic calibration of the 206 cm⁻¹ peak position is appropriate for common practice in quartz-inclusion Raman geobarometry. The pressure difference between the two calibrations is about 2 %; it does not significantly affect the pressure estimates but improves the self-consistency in the calculations of residual differential stress on inclusions that are very sensitive to minute uncertainties (see below). The Grüneisen plot (Fig. 3c) shows that the assumption of the constant Grüneisen parameters (Angel et al., 2019) is not valid in quartz, except for the 464 cm⁻¹ peak. The Grüneisen parameter (γ) dependence on volume is taken into account in fitting the data (Table 1) with two parameters (Reynard et al., 2012), the ambient pressure value γ₀, and its volume dependence q expressed as $q=\left(\partial \ln \mathit{\gamma }/\partial \ln V\right)$. Pressure was converted to volume using the equation of state (Angel et al., 1997). Non-linear Grüneisen parameters contribute to discrepancies between the present results and first-principles calculations on which constant γ were fitted (Murri et al., 2019).

5.2 Non-hydrostatic compression

The relative shifts of the 128 and 464 cm⁻¹ peaks with respect to the 206 cm⁻¹ peak are remarkably constant between present experiments at 295 K and those at 4 K (Tekippe et al., 1973). Absolute values of stress-induced shifts are higher (∼1.5 factor) in the present ambient temperature experiment than at 4 K where quartz is stiffer and has a lower volume corresponding to that of compression at about 0.25 GPa at ambient temperature. Shifts from the first-principles calculations (Murri et al., 2019) are compared to the experimental ones through the Grüneisen parameters listed in Table 1. To derive the Grüneisen tensor components from the present measurements, one needs to convert between stress and strain via the stiffness tensor. Assuming σ₁=σ₂ due to symmetry, the Raman shift is related to stress as follows:

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }\mathit{\nu }=\mathrm{2}{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}{\mathit{\sigma }}_{\mathrm{1}}+{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}{\mathit{\sigma }}_{\mathrm{3}},\end{array}$$

where ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}$ and ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}$ are directly fitted from the experiment (Table 1) and assumed constant. From linear elasticity, we have

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{1}}=\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right){\mathit{\epsilon }}_{\mathrm{1}}+C_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{3}},\text{(6)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{3}}=\mathrm{2}{C}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{1}}+C_{\mathrm{33}}{\mathit{\epsilon }}_{\mathrm{3}},\end{array}$$

where C_ij are the components of the elastic stiffness tensor. The shear components of strain and stress are not present because the crystal symmetry is preserved (i.e., σ₁=σ₂), and the off-diagonal strain components are all zero because $C_{\mathrm{41}}=-C_{\mathrm{42}}$. Substituting σ₁ and σ₃ into Eq. (4), and collecting common terms for ε₁ and ε₃, we have

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{\Delta }\mathit{\nu }& =\mathrm{2}\left[{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)+{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}{C}_{\mathrm{13}}\right]{\mathit{\epsilon }}_{\mathrm{1}}\\ & +\left[\mathrm{2}{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}{C}_{\mathrm{13}}+{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}{C}_{\mathrm{33}}\right]{\mathit{\epsilon }}_{\mathrm{3}}.\end{array}\end{array}$$

The Grüneisen tensor components are obtained by dividing the terms in the square bracket by the peak position at zero stress ν₀:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}=\left[{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)+{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}{C}_{\mathrm{13}}\right]/{\mathit{\nu }}_{\mathrm{0}},\text{(9)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{3}}=\left[\mathrm{2}{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}{C}_{\mathrm{13}}+{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}{C}_{\mathrm{33}}\right]/{\mathit{\nu }}_{\mathrm{0}}.\end{array}$$

Using the values for ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}$, ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}$ obtained here and published values of C_ij (Heyliger et al., 2003), we obtain the values reported in Table 1.

Shifts with uniaxial stress from first-principles calculations (Murri et al., 2019) are roughly consistent with experimental values for the 464 cm⁻¹ peaks and divergent for the 128 and 206 cm⁻¹ modes. These differences are possibly due to shortcomings of the first-principles calculations, such as unaccounted for anharmonic effects, or the assumption of constant Grüneisen parameters used in fitting theoretical results. Constant Grüneisen parameters are inconsistent with the pressure shifts of the 128 and 206 cm⁻¹ peaks (Fig. 3c). It shows that the discrepancy between theory and experiments, although reasonable given the absence of fitted parameters in the first-principles models, is significant when trying to apply results to quartz piezobarometry. In the following section, we will therefore reexamine the relationships between Raman peak positions and pressure and differential stress for quartz inclusions with the experimental shifts determined here.

6 Applications to piezospectroscopy of quartz inclusions

A spherical monocrystalline inclusion is subjected to triaxial normal stresses with the principal components σ₁, σ₂, and σ₃. The Raman shift ν is a function of these principal stresses and their orientation with respect to that of the crystal. Here, we set σ₃ parallel to the c axis of the inclusion. For the 206 cm⁻¹ peak, the shifts with differential stress are similar along and perpendicular to the c axis, and the Raman frequency will depend only on the residual pressure P defined as the first invariant of the stress tensor ($P=\frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{3}}}{\mathrm{3}}$). It is only for the 464 and 128 cm⁻¹ peaks that the different stresses will influence the Raman shift in response to both P and the differential stresses. Since the splitting of the 128 cm⁻¹ peak is not observed in natural quartz inclusions, it is assumed that there is no significant difference between σ₁ and σ₂.

With a uniaxial differential stress along the major symmetry axis of the quartz inclusion, the stress tensor can be written as follows:

$$\begin{array}{}\text{(10)}& \left(\begin{array}{ccc}P-\mathrm{1}/\mathrm{3}\mathit{\sigma }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& P-\mathrm{1}/\mathrm{3}\mathit{\sigma }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& P+\mathrm{2}/\mathrm{3}\mathit{\sigma }\end{array}\right),\end{array}$$

where $\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{3}}-{\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_c-{\mathit{\sigma }}_a$ is the residual differential stress in the inclusion, where a and c refer to the crystallographic axes of quartz. The shift Δν_nh of a Raman peak of wavenumber ν with the differential stress is

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }{\mathit{\nu }}_{\mathrm{nh}}=\mathrm{\Delta }\mathit{\nu }-\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{h}}=\mathrm{2}/\mathrm{3}\mathit{\sigma }\left[{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}-{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}\right],\end{array}$$

where Δν_h is the hydrostatic shift from Eq. (2), Δν is the measured shift with respect to ambient conditions, and ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_i}\right)}_{{\mathit{\sigma }}_j}$ are fitted slopes to the data under uniaxial compression. It is worth noting that while the shifts with pressure required fitting with a quadratic expression, the shifts with stress ${\left(\frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_i}\right)}_{{\mathit{\sigma }}_j}$ are assumed linear (Figs. 3 and 4, Table 1).

Combining Eqs. (11) and (2) and solving for the stress, we obtain the following:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathit{\sigma }& \approx \mathrm{3}/\mathrm{2}\left[\mathrm{\Delta }\mathit{\nu }-\left(-A+(A^{\mathrm{2}}+\mathrm{4}B{P}_{\mathrm{206}}{)}^{\mathrm{1}/\mathrm{2}})/\mathrm{2}B\right)\right]/\\ & \left[{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{1}}}-{\left({\displaystyle \frac{\partial \mathit{\nu }}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\right)}_{{\mathit{\sigma }}_{\mathrm{3}}}\right],\end{array}\end{array}$$

where P₂₀₆ is the pressure obtained from the position of the 206 cm⁻¹ peak, and A and B are the fitted coefficients from the hydrostatic experiment given in Table 1. The non-hydrostatic stress can be obtained from the positions of both the 464 and 128 cm⁻¹ peaks, and checked for self-consistency. Using the values in Table 1, where the average slopes for the two split components are used for the 128 cm⁻¹ peak, stresses (in GPa) are obtained:

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}\mathit{\sigma }\left(\mathrm{464}\right)& \approx \mathrm{1.5}\times \left[\mathrm{\Delta }\mathit{\nu }\right(\mathrm{464})\\ & -(-\mathrm{0.107}+({\mathrm{0.107}}^{\mathrm{2}}+\mathrm{0.0032}{P}_{\mathrm{206}}{)}^{\mathrm{1}/\mathrm{2}})/\\ & \mathrm{0.0016}]/(-\mathrm{2.7}),\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}\mathit{\sigma }\left(\mathrm{128}\right)& \approx \mathrm{1.5}\times \left[\mathrm{\Delta }\mathit{\nu }\right(\mathrm{128})\\ & -(-\mathrm{0.127}+({\mathrm{0.127}}^{\mathrm{2}}+\mathrm{0.0204}{P}_{\mathrm{206}}{)}^{\mathrm{1}/\mathrm{2}})/\\ & \mathrm{0.0102}]/(\mathrm{1.8}).\end{array}\end{array}$$

To summarize, discrepancies between pressures obtained from the three Raman peaks are expected with increasing residual differential stress. The position of the 206 cm⁻¹ peak is independent of differential stress σ and gives the hydrostatic pressure. If the residual differential stress is non-zero, residual pressures obtained from the 128 and 464 cm⁻¹ peaks should be shifted with respect to values obtained with the 206 cm⁻¹ peak by amounts of opposite signs. A routine for calculating pressure and stresses from the Raman peak positions is provided in the supplementary Excel spreadsheet. It is worth noting that the calibration established here directly translates the measured shifts of Raman peaks into stress, the geologically relevant variable, while approaches using the Grüneisen tensor translate peak shifts into strains that have to be converted to stresses using the elastic tensor.

As Raman measurements are subject to statistical fluctuations, so are the pressure (Eq. 1) and the two stresses calculated from Eqs. (13) and (14) likely to differ for a given inclusion (Table S1). Uncertainty in pressure will arise from uncertainty in the position of the 206 cm⁻¹ peak. Uncertainty on stress will arise from uncertainty in the position of the 128 and 464 cm⁻¹ peaks and also in the position of the 206 cm⁻¹ peak that is used to calculate pressure in Eqs. (13) and (14). A sensitivity analysis was performed to quantify the effect of small fluctuation in the Raman peak position on the calculated differential stress (Fig. S2). It is demonstrated that fluctuations as small as 0.2 to 0.4 cm⁻¹ in the Raman band position can lead to a scattering of up to 0.5 GPa in differential stress. This explains the wide scattering of natural and experimental data as shown in Figs. 5 and 6. It also illustrates the origin of the anti-correlation between σ(128) and σ(464) observed in Figs. 5 and 6 from the larger uncertainty on the 206 cm⁻¹ peak position than on other peaks.

Figure 5 Pressures and differential stresses obtained from experimental quartz inclusions (Bonazzi et al., 2019). (a) Residual pressures from thermoelastic modeling (Fig. 2) as a function of pressures from the 206 cm⁻¹ peak for selected inclusions at the indicated conditions of formation. Empty rectangles show the mean value interval at 95 % confidence level (2 standard error). Solid 1:1 line and dashed lines with 0.1 GPa offset are shown. Horizontal bar shows the uncertainty associated with 1 cm⁻¹ uncertainty on the 206 cm⁻¹ peak position. (b) Pressures from the 128 and 464 cm⁻¹ peaks as a function of pressure from the 206 cm⁻¹ peak. Residual pressures are in general agreement, except for a group of experimental inclusions that show significant departure from the 1:1 line. Ellipse shows the uncertainties associated with the 0.5 cm⁻¹ uncertainties in the 128–464 cm⁻¹ peak positions. (c) Residual differential stresses from the 128 and 464 cm⁻¹ peaks are consistent within uncertainties, except for the group of inclusions departing from the 1:1 line by more than 1 GPa (see text). Ellipse shows the uncertainties associated with 1 and 0.5 cm⁻¹ uncertainties on the 206 and 128–464 cm⁻¹ peak positions, respectively. Anti-correlated uncertainties are associated with the 1 cm⁻¹ uncertainty on the 206 cm⁻¹ peak position that has opposite effects on the differential stresses calculated for the 128 and 464 cm⁻¹ peaks using Eqs. (13) and (14), respectively (see text and Fig. S2).

Figure 6 Stresses obtained from natural quartz inclusions. (a) Pressures from the 128 and 464 cm⁻¹ peaks deviate from pressure from the 206 cm⁻¹ peak with roughly opposite effects of similar magnitude in Papua rocks, and differences are negligible in Syros blueschist. (b) Residual differential stresses from the 128 and 464 cm⁻¹ peaks. Mean values and confidence interval at 95 % are shown as colored ellipses. The dispersion of individual measurements is accounted for by systematic uncertainties on Raman measurements (black empty ellipse, same as in Fig. 5c). Residual stresses are negative in Papua New Guinea and Norway rocks and close to zero in Syros blueschists. Data from Cisneros et al. (2021) for inclusions in garnets (filled symbols) and epidotes (empty symbols) in blueschists from Syros (Cyclades islands, Greece), Gonzalez et al. (2019) for rocks from Papua New Guinea, and Zhong et al. (2019) for eclogitic inclusions from the Bergen Arcs (Norway). Eleven eclogitic inclusions from Norway were remeasured for this study with systematic comparison to the reference spectrum. Large uncertainties on the 128 cm⁻¹ peak position on this dataset are associated with the high cutoff of the notch filter of the spectrometer in Berlin.

If a sufficient number of inclusions is measured, and if they belong to the same population (i.e., if they formed in a single event under similar pressure and temperature conditions), the average pressure and stress values can be obtained with a greater accuracy since the standard error of the mean will be greatly reduced in spite of the fluctuations. An alternative approach would be to use the over-determination of the system (two unknowns and three independent equations or more if more peaks are used) to obtain the best-fit values of pressure and differential stress, or the strains (Bonazzi et al., 2019; Murri et al., 2019). It is our choice to compare the stresses obtained from Eqs. (13) and (14) because their sensitivity to uncertainties allows discussing the self-consistency of the Raman measurements that may otherwise be blurred if a single best-fit value is calculated for each inclusion.

6.1 Comparison with experimental inclusions

Bonazzi et al. (2019) synthesized quartz inclusions in almandine at 2.5 GPa and 800 ^∘C and 3 GPa and 775 ^∘C. The present calibration was tested on their selection of inclusions based on optical criteria (absence of cracks, sufficient distance from garnet rim, and other inclusions, …). For most of these inclusions, the residual pressures inferred from the 206 cm⁻¹ peak position are consistent with those expected from elastic modeling of quenching and decompression from the experimental equilibration conditions (Fig. 5a).

Several inclusions synthesized at 3 GPa display residual pressures lower than expected and that are inconsistent between the three peaks (Fig. 5b), indicating large differential residual differential stresses (Fig. 5c). For those, the inferred residual differential stresses are of opposite signs when estimated using the 128 or 464 cm⁻¹ peaks, with values of σ between $\sim -\mathrm{0.5}$ and 2 GPa using the 464 cm⁻¹ peak, and between ∼0.5 and −3 GPa using the 128 cm⁻¹ peak. It indicates that these inclusions are under a stress pattern that is not consistent with simple elastic deformation. The large residual differential stresses are likely due to unobserved defects around inclusions. Bonazzi et al. (2019) proposed a correction for similar effects in strains using the first-principles strain-induced shifts (Murri et al., 2019).

We propose to use the inconsistency in residual differential stress calculated here directly from Raman peak positions and Eqs. (1), (13), and (14) as a guide to eliminate those inclusions under complex stress state from the analysis of residual pressures. Thus, we discarded from the analysis inclusions with an absolute difference in residual differential stress from the 128 and 464 cm⁻¹ peak positions larger than a value of 1 GPa as defined from the sensitivity analysis (see Table S1 and Fig. S2). This changes the average residual pressure from 1.06(3) to 1.11(3) GPa. Both residual pressure averages are within uncertainties of the value of 1.08 GPa expected from the elastic model. The average value of residual differential stresses after selection of 0.00(10) and −0.03(11) is close to that of 0.0 GPa expected from elastic modeling (Fig. 2), when the one prior to selection is $\sim -\mathrm{0.5}\left(\mathrm{2}\right)$ and 0.4(2) GPa for the 128 and 464 cm⁻¹ peaks, respectively.

The inclusions synthesized at 2.5 GPa yield residual pressures of 0.82(1) after selection with the above-defined criterion from differential stresses using the 128 and 464 cm⁻¹ peak positions (Fig. 5b, Table S1). They are consistent with those of the elastic model of 0.86 GPa. Residual pressures from the 128 cm⁻¹ peak are systematically lower by ∼0.1 GPa. As a result, residual differential stress is −0.01(4) for the 464 cm⁻¹ peak, and of $\sim -\mathrm{0.4}\left(\mathrm{1}\right)$ for the 128 cm⁻¹ peak, where a value of −0.05 GPa is expected from the elastic model at 2.5 GPa and 800 ^∘C (Fig. 2). The systematic shift of residual stress values from the 1:1 line (Fig. 5c) would correspond to a systematic deviation of only 0.5 cm⁻¹ in the 206 peak position for this particular set of inclusions, due for instance to interference with garnet peaks. The method use here allows detecting and illustrating the effects of such a minute deviation.

6.2 Comparison with natural inclusions

We used recent datasets on natural inclusions from three different metamorphic suites: (1) blueschists from Syros in garnet formed at ∼1.4–1.7 GPa and 500–550 ^∘C and in retrograde epidote grown between ∼1.3–1.5 GPa and 400–500 ^∘C as well as ∼1.0 GPa and 400 ^∘C (Cisneros et al., 2021); (2) Holsnøy eclogite from the Bergen Arcs (Norway), where inclusions formed at 1.4–1.6 GPa and 680–760 ^∘C, consistent with previous estimates based on phase equilibria (Zhong et al., 2019; Bhowany et al., 2018); (3) gneiss from Papua New Guinea with formation conditions estimated from Ti-in-quartz and quartz-in-garnet inclusions as ∼1.0 GPa and 600 ^∘C, which are interpreted as the P–T conditions of garnet growth and entrapment of quartz inclusions (Gonzalez et al., 2019).

The residual pressures estimated from the present calibration are within uncertainties of those published using earlier calibrations. The major difference comes from the possibility of estimating the residual differential stress here and to use it as a criterion for validating measurements (Sect. 6.1). Residual pressures from the 128 and 464 cm⁻¹ peaks are respectively higher and lower than hydrostatic pressure from the 206 cm⁻¹ peak in rocks from Papua New Guinea, indicating significant residual differential stress. They are almost indistinguishable within uncertainties for blueschists of Syros, indicating little residual differential stress (Fig. 6a, Table S1).

Residual differential stresses σ₁₂₈ and σ₄₆₄ are self-consistent and low, with average values of less than 0.3 GPa (Fig. 6b, Table S1). At these deviatoric stress levels, symmetry breaking effects are likely negligible (Murri et al., 2022). Average values are well defined, except in inclusions from Norway, where the 128 cm⁻¹ peak position is imprecise due to spectrometer configuration and because of the small number of measurable inclusions. The dispersion of individual values is accounted for by uncertainties in the Raman peak position (Figs. 6b and S2). The dispersion of residual differential stresses is smaller in natural inclusions than in experimental ones, likely in response to longer growth rates and potential anelastic relaxation times than in the relatively short experiments followed by fast quenching and pressure release. Such processes may affect the reconstructed initial P–T conditions and would deserve further investigation. With the method and calibration proposed here, all but five measured natural inclusions (see Table S1) fall within the conditions defined for the calibration of uniaxial differential stress along the major symmetry axis of the quartz inclusion at the beginning of Sect. 6. The method using strain anisotropy resulted in the rejection of 74 inclusions out of 92 (Gonzalez et al., 2019); a high rejection rate may stem from the different rejection criteria.

Average values of residual differential stress from both the 128 and 464 cm⁻¹ peaks are −0.04(4), −0.18(3), and −0.16(14) GPa for Syros, Papua New Guinea, and Norway, respectively, with a similar trend to values expected from the elastic model of −0.06(3), −0.13(3), and −0.14(4), respectively (Fig. 2). Taken at face value, residual differential stresses in Syros and Papua New Guinea translate into temperatures of 340–600 and 640–830 ^∘C, respectively, to be compared with 500–550 and 600 ^∘C estimated from independent phase equilibria and geothermometers (Cisneros et al., 2021; Gonzalez et al., 2019). For the Holsnøy eclogite (Norway), the mean value of residual stress yields ∼750 ^∘C within the range of 680–760 ^∘C from previous estimates (Zhong et al., 2019), although the temperature range is not precisely defined owing to the large uncertainty on residual stress. It suggests that independent constraint on temperature of entrapment may be obtained from residual differential stresses on inclusions. This would require reducing uncertainties by measuring residual pressures and differential stresses on a large number of inclusions, systematic calibration of Raman measurements with Ne or Hg emission lines (Sect. 2), and calibration on samples formed in well-characterized metamorphic conditions. The effect of large numbers of analyses on improving the accuracy of residual stress measurements is clearly seen when comparing mean values and uncertainties (Fig. 6b) from small (Norway, 11 inclusions), medium (Syros, 22), and large (Papua New Guinea, 92) datasets.

7 Concluding remarks

The present calibration can be used to directly estimate the residual pressures and differential stresses using simple relationships established between Raman shifts, pressure, and uniaxial stresses. A sensitivity analysis on calculated uniaxial stresses is proposed that allows detecting and evaluating systematic deviations in the Raman data. Accurate average pressure and deviatoric stress are obtained after checking that populations of inclusions belong to a single population formed at similar conditions and averaging over this population. It can be applied provided that Raman peak positions are published in studies using elastic geobarometry.

Hydrostatic pressure is obtained from the 206 cm⁻¹ peak position that is insensitive to residual stress. Self-consistency of residual differential stresses obtained from the 128 and 464 cm⁻¹ peak positions is an objective filter to detect inclusions that are under stress conditions where assumptions used for elastic modeling of entrapment conditions do not apply. Residual stress effects on the pressures determined with 128 and 464 cm⁻¹ peaks are of similar magnitude and opposite sign; hence the average of these two values is a good estimate of the hydrostatic pressure if the 206 cm⁻¹ peak cannot be measured accurately due to interference with host mineral peaks. Using the 128 and 464 cm⁻¹ peak alone, or averaging either 128 and 206 or 206 and 464 cm⁻¹ peaks, can induce systematic bias in the residual pressure determination.

Residual differential stresses observed in suites of metamorphic rock inclusions suggest that they depend on temperature at the time of entrapment. Reduction of uncertainties and calibration on rocks with various conditions of formation would be necessary for use as a quantitative geothermometer. Large numbers of inclusions should be measured and uncertainties on Raman data reduced with the systematic use of internal calibration lines. Reduction of uncertainties may also require Raman measurements at combined high pressure and high differential stress, currently an experimental challenge, and elastic constant determinations at combined high pressure and temperature to limit extrapolation in elastic model predictions.