Measurement of strains in zircon inclusions by Raman spectroscopy

The discovery of coesite inclusions in garnet in the Dora-Maira massif (Chopin, 1984) was “the little mineral that changed everything” (Tremblay, 2013). The discovery clearly indicated that continental crust can experience pressures in excess of 2.8 GPa and be subsequently exhumed. The question of whether these rocks experienced lithostatic pressures of 2.8 GPa and thus depths of metamorphism of 90 km or more, or whether instead they experienced significant deviatoric stress that stabilised coesite with respect to quartz at much shallower depths, was raised almost immediately by Smith (1984) in his description of coesite inclusions in pyroxenes from the Western Gneiss Region of Norway. The origin of such over-pressures could be on the micro-scale, generated by grain–grain interactions between minerals with different thermo-elastic properties, on the scale of the outcrop as indicated by the strongly heterogeneous deformation of the Dora-Maira rocks at peak conditions (Henry et al., 1993) or be regional in scale due to tectonic overpressure (e.g., Schmalholz & Podladchikov, 2014; Gerya, 2015).

Conventional thermo-barometry is challenged in ultra-high-pressure metamorphic (UHPM) terrains because it relies on thermodynamic equilibrium being attained in the rocks, and the textures and peak phase assemblage being preserved during exhumation. A complementary method of geobarometry, based on the stress state of inclusions in high-pressure minerals, has been developed over the past two decades, based on the original ideas of Rosenfeld & Chase (1961). The basic idea is that a mineral such as coesite or zircon trapped inside a garnet that has grown at UHPM conditions will exhibit residual stresses when examined at room conditions, as a result of the contrast in thermoelastic properties (the equations of state [EoS]) of the host and the inclusion mineral. In simple terms, the inclusion mineral is constricted by the host and cannot expand freely on exhumation as would a crystal of the same mineral in the rock matrix. Knowledge of the EoS of the minerals allows possible entrapment conditions of the inclusion to be calculated and thus provides a second type of geobarometer completely independent of chemical equilibrium (Rosenfeld & Chase, 1961). Recent developments in theory have expanded the original isotropic models for calculating entrapment conditions (e.g., Gillet et al., 1984; Zhang, 1998; Angel et al., 2014) to allow for the effects of the shapes and the local environment of the inclusions (Campomenosi et al., 2018; Mazzucchelli et al., 2018). The effects of the anisotropy of the inclusion minerals trapped within a host that is elastically isotropic can also be calculated by extension of these methods. Calculations of entrapment conditions are based on the strain state of the inclusion mineral still buried in its host relative to a free crystal of the same mineral at room conditions. The challenge is that small faceted inclusions can exhibit significant strain gradients (e.g., Campomenosi et al., 2018; Murri et al., 2018) and these cannot be measured by X-ray laboratory diffraction because the large X-ray beams illuminate all or most of the inclusion and therefore return some average value of the unit-cell parameters of the inclusion and hence the strains. Such average values cannot be used to determine entrapment conditions of the inclusion. As shown in the original work on the coesite inclusions (Chopin, 1984; Smith, 1984), Raman spectra can be collected from much smaller volumes (of the order of 1 μm³) and allow local strains to be determined if the values of the phonon-mode Grüneisen tensors for the inclusion phase are known (Angel et al., 2019). The values of the components of these tensors can be determined by calculating the Raman spectra of the inclusion mineral under different strain conditions with Density Functional Theory (DFT) (Murri et al., 2018, 2019). Coesite, being monoclinic, is a challenge for this approach because it requires many more simulations of different strain states to be performed than are required for a uniaxial mineral such as quartz, and the theory for calculating the elastic relaxation of monoclinic mineral inclusions has yet to be developed. However, zircon is also ubiquitous as an inclusion mineral in garnets in UHPM rocks such as those found in Dora-Maira (Schertl et al., 1991), and has the potential to be used as an elastic geobarometer as well as providing the opportunity for radiometric dating of the time of entrapment. In this paper we present the results of DFT simulations of zircon under a range of deviatoric strains in order to determine the components of the phonon-mode Grüneisen tensors of zircon. We show that these values correctly predict the change in the Raman shifts of modes under P, but not under T, and we show that this contrast in behaviour to that of quartz is due to differences between the anisotropy of the thermoelastic properties of the two minerals. The resulting tensor components can be used to determine the strains in zircon inclusions from measurements of their Raman spectra and thus ultimately their entrapment conditions.

The geometry of the zircon structure was optimised by ab initio HF/DFT simulations that were performed with the CRYSTAL17 code (Dovesi et al., 2018) using the same hybrid WC1LYP Hamiltonian (Wu & Cohen, 2006), basis sets and parameters that were used for previous simulations of zircon under hydrostatic pressure (Stangarone et al., 2019). One full optimization of both the unit-cell parameters and the atomic positions was performed at zero pressure without any constraints except those imposed by the I4₁/amd space-group symmetry. This provided the unstrained reference structure. Atomic fractional coordinates of zircon were then optimized at 47 different strained configurations with fixed unit-cell parameters and the a and b unit-cell parameters kept equal to maintain the tetragonal symmetry. The unit-cell parameters were chosen to cover the range of strains, ε₁ = ∆a/a, ε₂ = ∆b/b, and ε₃ = ∆c/c, expected to be found in zircon inclusions in garnets. All of the computed geometries are at static equilibrium (i.e., at 0 K, with no vibrational effects included). For each optimized geometry the normal modes of vibration and their wavenumbers were calculated within the limit of the harmonic approximation as described in Stangarone et al. (2019). The average mismatch between phonon-mode wavenumbers calculated in this way and experimental data for a wide variety of minerals is 5 cm⁻¹ (Prencipe, 2019). The calculated phonon-mode frequencies and the structural data (atomic coordinates and unit-cell parameters) for the 48 optimisations of zircon are reported in the crystallographic information file (CIF) deposited as Supplementary Material linked to this article and freely available at https://pubs.geoscienceworld.org/eurjmin/.

The structure of zircon contains Si in tetrahedral coordination by oxygen, and Zr in 8-fold coordination by oxygen in the form of triangular-faced dodecahedra (Fig. 1). The Si and Zr atoms occupy positions that are fixed by the I4₁/amd space-group symmetry. The Si and Zr are separated along [001] by a distance of exactly c/2, and their coordination polyhedra are joined by a shared O–O edge which form continuous columns of edge-sharing polyhedra which are parallel to the c-axis. The Zr–O bonds to the O atoms in the shared edge are 0.13 Å longer than those to O not involved in edge-sharing with the SiO₄, and the O–Si–O angle to the same shared edge is only 96.5° (e.g., Kolesov et al., 2001; Finch & Hanchar, 2003). These structural features indicate that there is a strong repulsive interaction between the neighbouring Zr and Si atoms along [001] (e.g., Kolesov et al., 2001), which makes this direction more than twice as stiff as the perpendicular directions under high pressure (e.g., Hazen & Finger, 1979; Pina-Binvignat et al., 2018).

In the absence of experimental data for the structure and unit-cell parameters of zircon under deviatoric stress and strain, the only test of the validity of our DFT results is to compare the simulations under hydrostatic stress with experimental data. In a previous paper (Stangarone et al., 2019) we have shown that the simulations under hydrostatic stress slightly over-estimate the unit-cell volume by 0.6%, and the Si–O and the shorter Zr–O bond distances by 0.7–0.8%. This arises from the small error in calculating electron self-interactions in the core regions of the atoms (Cremer, 2001) in the simulations. Nonetheless, the important distortions of the SiO₄ tetrahedra are well-reproduced, the predicted structural evolution with pressure is consistent with the limited experimental data (Hazen & Finger, 1979), and the elastic anisotropy of the structure is correct in the DFT simulations. Similarly, both the Raman shifts calculated at zero pressure (Table 1) and their change with pressure are in good agreement with experimental data (Syme et al., 1977; Pina-Binvignat et al., 2018).

Figure 2 is a series of contour plots of structural parameters of zircon to show how they respond to the strains ε₁ and ε₃. Under non-hydrostatic conditions the length of the O–O shared edge between the SiO₄ tetrahedra and the ZrO₈ dodecahedra, which is perpendicular to the c-axis (Fig. 1), is almost unchanged under compression of the c-axis, but it responds to strains ε₁ in the a–b plane (Fig. 2a). Thus, when the c-axis is compressed (ε₃ < 0) the distance between adjacent Si and Zr is reduced by the same amount as the unit-cell strain, the O–O edge moves towards the Si, and thus the corresponding O–Si–O angle opens and increases towards the ideal tetrahedral angle of 109.47°. As a consequence, the distortion of the SiO₄ tetrahedra as measured by their angle variance (Robinson et al., 1971) decreases significantly with compression along the c-axis, but is not changed by strains ε₁ (Fig. 2b). The length of the Zr–O long bonds to the O atoms in the shared edge is also therefore dependent solely on the strain ε₃ (Fig. 2c). In contrast, because the shorter Zr–O bonds to polyhedral edges shared between adjacent ZrO₈ polyhedra (Fig. 1) are sub-parallel to (001), their length is dependent mostly on the strain ε₁ and is insensitive to ε₃ (Fig. 2d). It is the combination of these two different responses of the long and short Zr–O bonds to applied strains, induced by the Zr–Si repulsion and the stiffness of the O–O shared edge, that leads to the volumes of both the SiO₄ tetrahedra and the ZrO₈ polyhedra (Fig. 2e, f) following the isochors of the unit-cell volume (lines of constant 2ε₁ + ε₃) while their distortions (Fig. 2g) are almost entirely dependent on the deviatoric strain (ε₃−ε₁). In summary, the anisotropy of response of this O–O shared edge, which must be due to the Si–Zr repulsion, is responsible for the anisotropy of response of the other structural parameters.

The wavenumbers of all of the phonon modes in zircon calculated for the structure at the static equilibrium (zero T and P and with no zero-point pressure included) are listed in Table 1 as $ω0,DFTm$⁠. As previously reported in Stangarone et al. (2019), they are in reasonable agreement with the values estimated by extrapolating to 0 K the experimental values measured over the range of 80–1400 K (full details in the Supplementary Material of the present paper). The level of discrepancy is similar to that obtained in previous simulations of the Raman spectra of zircon (Sheremetyeva et al., 2018) and other silicate minerals (Prencipe, 2019). For ease of understanding we will refer to the modes by their wavenumbers as measured experimentally at room conditions, and their symmetries (Table 2).

Consideration of Table 2 together with the isoshift lines of selected modes in Fig. 3 illustrates some general points about phonon-mode Grüneisen tensors. The different slopes of the isoshift lines for different modes in Fig. 3 correspond to significantly different values for $γ1m$ and $γ3m$⁠. The majority of modes have positive values for both $γ1m$ and $γ3m$ which means (Eq. (1)) that the phonon-mode wavenumbers increase with negative strains and decrease with positive strains (Fig. 3a, b). This corresponds to the normal behaviour observed in minerals, that the wavenumbers of phonon modes increase under pressure and decrease with increasing temperature. There is also a general trend for the values of $γ1m$ and $γ3m$ of the modes in zircon to become more similar in magnitude with increasing wavenumber of the mode (Table 2). This decrease in anisotropy of the mode Grüneisen tensors is due to the decreasing contribution of Zr atom motions to the higher frequency modes which become dominated by the “internal” modes of the SiO₄ tetrahedra.

In contrast, the B_1u mode near 131 cm⁻¹ (which is neither Raman nor infra-red active) corresponds to a rigid rotation of the SiO₄ tetrahedra around the c-axis, and has very large negative values ~ −6 for $γ1131$ and $γ3131$ (Fig. 3e) because it is the soft mode that drives the displacive phase transition from zircon to a high-pressure phase with lower symmetry (Stangarone et al., 2019). The large negative values indicate that the wavenumber of this mode decreases strongly with increasing pressure (Fig. 4). The B_2g mode near 265 cm⁻¹ and the E_g mode near 201 cm⁻¹ also have both phonon-mode Grüneisen tensor components negative, but they are small in magnitude (Table 2). This indicates that they are modes that soften because they either couple with the soft mode directly or with the strain. Direct coupling may be expected in this case, because both of these modes also involve rotations (or twists) of the SiO₄ tetrahedra around the c-axis as does the B_1u soft mode. Similarly, the E_g mode at 224 cm⁻¹ also involves rigid rotations of the SiO₄ tetrahedra (Table 2) and has $γ1224<0$ and $γ3224>0$⁠. The phonon-mode Grüneisen tensors therefore clearly distinguish the so-called “hard modes” (e.g., Salje, 1992; Salje & Bismayer, 1997) which have relatively small positive components $γim$ from soft modes that drive displacive phase transitions to low-symmetry phases at high pressure (large negative $γim$⁠), and the modes that are affected by, or support, the same phase transitions (small negative $γim$⁠).

This analysis is based (Eq. (1)) on a linear relationship between the changes in wavenumbers of the phonon modes and the strains, while the DFT simulations show (Fig. 3) that there are small non-linearities in the variation of wavenumbers with strain. It therefore should be understood that the values of $γ1m$ and $γ3m$ determined from the DFT simulations are average values over the range of strains included in the analysis. Whether the values are correct, and whether they can be used for larger strains, can only be determined by comparison with experimental data. We can use the known EoS and cell-parameter variation of zircon (Zaffiro, 2019) to calculate the strains relative to room conditions induced by P and/or T, and from the calculated strains predict the change in the phonon wavenumber from Equation (1). Figure 4 shows that the Grüneisen tensor approach correctly predicts the shifts of the phonon modes with pressure, including the E_g modes at 201 cm⁻¹ and 224 cm⁻¹ that are coupled to the soft B_1u mode. A limit to the linear relationship might be indicated by the small discrepancies between predicted and observed phonon-mode wavenumber changes when ∆ω^m > 20 cm⁻¹ or P > 7 GPa, whichever occurs first.

However, Fig. 5 shows that the experimental measurements of the phonon-mode wavenumbers at both low and high temperatures at room pressure are not predicted by the mode Grüneisen tensor components determined from the DFT simulations. For the hard modes the decrease in wavenumber with increasing temperature is underestimated, by different amounts for the various modes. For the phonon modes with negative $γim$ the positive strains associated with thermal expansion mean that the Grüneisen tensor approach predicts that these modes should exhibit a positive wavenumber change with increasing temperature, whereas the opposite is observed experimentally for the E_g mode at 201 cm⁻¹.

Phonon-mode Grüneisen tensor coefficients calculated from the results of DFT simulations predict the experimentally measured mode-wavenumber changes of zircon with P at room T (Fig. 4), which shows that the results of the DFT simulations are accurate. Therefore, the failure to correctly predict the mode shifts induced by T alone (Fig. 5) must arise from the application of the concept of the phonon-mode Grüneisen tensor. This result is very different from quartz, where the phonon-mode Grüneisen tensor components correctly predict the changes in the wavenumbers of Raman modes with both P and T, at least at conditions where the effects of the alpha-beta phase transition do not dominate the behaviour of quartz (Murri et al., 2018, 2019).

This contrast in the behaviour of the Raman modes of quartz and zircon can be understood as the consequence of differences in the anisotropies of their compressibility and thermal expansion.

Because we have eliminated the phonon-mode Grüneisen tensor components $γ1m$ and $γ3m$ from the equations without making any further assumptions, we have shown that when a crystal exhibits isotropy in α_i/β_i (Eq. (3)), the same phonon-mode Grüneisen tensor can explain the changes in phonon-mode wavenumbers induced by changes in either P or T, as found for many modes in quartz (Murri et al., 2018). Similarly, it will always apply to cubic materials, for which Equation (3) is always true.

However, the c-axis of zircon is considerably stiffer than the a-axis and also has a higher thermal expansion coefficient (Table 3), as a consequence of the Zr–Si repulsion across the shared O–O polyhedral edge. Therefore Equation (3) does not hold for zircon and therefore Equation (4) and the subsequent equations cannot be derived, which means that the thermally induced changes in phonon-mode wavenumbers in zircon are independent and different from the wavenumber changes induced by hydrostatic pressure. However, the changes in phonon-mode wavenumbers induced by pressure (Fig. 4), or strain, at any fixed temperature are still correctly predicted by the phonon-mode Grüneisen tensor, meaning that we can still use the concept to determine isothermal strains in crystals from the measured Raman shifts.

Our determination of the phonon modes of zircon under various strains which deviate from hydrostatic stress conditions has shown that, if the temperature is held constant, the changes in their wave numbers depends approximately linearly on the unit-cell strains, up to strains of a few % (Figs. 3, 4). This linear dependence is described by the components of the phonon-mode Grüneisen tensors of zircon, which we have reported for all of the Raman-active bands (Table 2). Because of the strong anisotropy of the ratio of linear thermal expansion to linear compressibility (the anisotropic thermal pressure tensor; Nye, 1957) in zircon, these phonon-mode Grüneisen tensors do not describe the temperature-induced shifts in the phonon-mode wavenumbers (Fig. 5). Nonetheless, they do correctly predict the changes in phonon wavenumbers with isothermal changes in hydrostatic pressure (Fig. 4), and should therefore also apply to conditions of non-hydrostatic stress. Thus, the values of $γ1m$ and $γ3m$ that we have determined for zircon (Table 2) can be used to map the residual strains in zircon inclusions trapped inside garnets, and other host minerals that are almost elastically isotropic.

For example, Fig. 6 shows the strains calculated from Raman spectra collected on two traverses across a zircon inclusion trapped about 1 cm from the rim of a 12 cm-diameter pyrope megablast from the whiteschists of the UHP Brossasco-Isasca Unit in the Dora-Maira massif originally described by Chopin (1984). The full Raman spectrum at each point was measured (Fig. 6b). The Raman shift of the 356 cm⁻¹E_g mode along the two traverses of the inclusion (Fig. 6c) is clearly different and this shows immediately that the strain in the inclusion is inhomogeneous. A total of seven Raman peaks from zircon are clearly resolved from the spectra of the garnet host (Fig. 6b). Of these, the three peaks at 200–225 cm⁻¹ and the peak near 439 cm⁻¹ exhibit small shifts with strains as the combination of relatively small coefficients of their phonon-mode Grüneisen tensors (Table 2) and small values of $ω0m$⁠. This combination of factors means that their measured shifts from the reference crystal do not significantly constrain the values of strains. On the other hand, the Raman-active modes near 356, 975 and 1009 cm⁻¹ are stronger (Fig. 6b) and exhibit larger shifts. The difference in the ratios $γ1m/γ3m$ for these three modes means that their isoshift lines have different slopes when plotted against the strains ε₁ and ε₃ (e.g., Fig. 3a, b). Therefore, the strains at any one point in the inclusion can be determined from the changes in wavenumber of these three Raman modes.

The determination of strains from the Raman spectra of inclusions is entirely dependent on the difference in peak wavenumber from that of an unstrained reference crystal. We therefore measured an unstrained reference crystal in air during the same measurement session and calculated the difference in wavenumbers as $∆ωexpm$= $ωexpm-ω0,expm$⁠. Strains were then calculated from the values of $∆ωexpm$by using Equation (1) in a least-squares procedure implemented in the stRAinMAN program (Angel et al., 2019) which is available at www.rossangel.net. The file zircon_gtensor.cif provided in the Supplementary Material for this paper and available on-line contains the values of $γ1m$ and $γ3m$ for all of the Raman active modes of zircon and can be used to perform this calculation. In this way, the strains and their gradients can be mapped across trapped inclusions by Raman spectroscopy (Fig. 6d). It is clear that the strains not only vary significantly between the different points, but the gradient in strains along the crystal is different for ε₁ and ε₃. This confirms that the strains vary inhomogeneously throughout the inclusion as a consequence of its shape (Mazzucchelli et al., 2018) and elastic anisotropy (Mazzucchelli, 2019), and consequently the inclusion is under inhomogeneous stress, with significant stress gradients.

The interpretation of the strains of the inclusion in terms of entrapment conditions can then be performed two ways. If the inclusion was approximately spherical, or ellipsoidal, then the strains would be anisotropic and homogeneous, and could be converted to stresses. The mean normal stress could then be used as the inclusion pressure within an isotropic calculation to yield the entrapment isomeke (Angel et al., 2014, 2017b) which should approximate the entrapment conditions. A more complete but more complex analysis for inclusions shaped as the one in Fig. 6 requires the anisotropic relaxation of the host–inclusion system to be calculated (Mazzucchelli, 2019) and then the individual strains at the centre of the inclusion, after correction for relaxation, can be combined to calculate the entrapment P and T.

This work was supported by ERC starting grant 714936 “True Depths” to Matteo Alvaro. MA was also supported by the F.A.R.E.–M.I.U.R. grant R164WEJAHH. We thank Mara Murri and Mattia Mazzucchelli for discussions and Gabriele Zaffiro for the EoS of zircon.