We have carried out *ab initio* hybrid Hartree-Fock/Density Functional Theory simulations to determine the structure and vibrational modes of zircon, ZrSiO_{4}, as a function of different applied strains. The changes in phonon-mode wavenumbers are approximately linear in the unit-cell strains, and have been fitted to determine the components of the phonon-mode Grüneisen tensors of zircon which reproduce the change in measured Raman shifts with pressure. They can therefore be used to convert Raman shifts measured from zircon inclusions in metamorphic rocks into strains that in turn can be used to determine the metamorphic conditions at the time that the inclusion was trapped. Due to the strong anisotropy in the thermal pressure of zircon, the phonon-mode Grüneisen tensor is not able to reproduce the temperature-induced changes in Raman shifts. Because zircon inclusions are normally measured at room conditions this does not prevent the calculation of their entrapment conditions.

## Introduction

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^{3}) 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.

## Computational details

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 *I*4_{1}/*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, ε_{1} = ∆*a*/*a*, ε_{2} = ∆*b*/*b*, and ε_{3} = ∆*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^{−1} (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/.

## Results

### Structure

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 *I*4_{1}/*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_{4}, 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_{4} 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 ε_{1} and ε_{3}. Under non-hydrostatic conditions the length of the O–O shared edge between the SiO_{4} tetrahedra and the ZrO_{8} dodecahedra, which is perpendicular to the **c**-axis (Fig. 1), is almost unchanged under compression of the **c**-axis, but it responds to strains ε_{1} in the **a–b** plane (Fig. 2a). Thus, when the **c**-axis is compressed (ε_{3} < 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_{4} 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 ε_{1} (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 ε_{3} (Fig. 2c). In contrast, because the shorter Zr–O bonds to polyhedral edges shared between adjacent ZrO_{8} polyhedra (Fig. 1) are sub-parallel to (001), their length is dependent mostly on the strain ε_{1} and is insensitive to ε_{3} (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_{4} tetrahedra and the ZrO_{8} polyhedra (Fig. 2e, f) following the isochors of the unit-cell volume (lines of constant 2ε_{1} + ε_{3}) while their distortions (Fig. 2g) are almost entirely dependent on the deviatoric strain (ε_{3}−ε_{1}). 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.

### Vibrational modes

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 $\omega 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).

*E*

_{g}mode near 356 cm

^{−1}clearly responds more to ε

_{1}than ε

_{3}, whereas the opposite is true for the

*E*

_{g}mode near 224 cm

^{−1}(Fig. 3b, c). In general, the isoshift lines are not parallel to isochors, which are lines of constant 2ε

_{1}+ ε

_{3}(

*e.g.*, Fig. 3b). The variation of slopes of the isoshift lines between modes also clearly means that they are not, in general, parallel to isobars which are lines of equal average stress (2σ

_{1}+ σ

_{3})/3. This is a general behaviour for crystals, which we have previously described for quartz (Murri

*et al.*, 2019). It indicates that the phonon-mode wavenumbers change approximately linearly with applied strains and can therefore be described by the phonon-mode Grüneisen tensor γ

^{m}(

*e.g.*, Ziman, 1960; Barron

*et al.*, 1980; Angel

*et al.*, 2019) which is a linear model. For uniaxial crystals in which the symmetry is not broken by the strain (as we have for our DFT simulations), ε

_{1}= ε

_{2}and the wavenumber change ∆ω

^{m}is given by:

*et al.*, 2019). The values of $\gamma 1m$ and $\gamma 3m$ (Table 2, and

*cif*in Supplementary Material) were determined for each mode by a least-squares fit of Equation (1) to the values of ∆ω

^{m}determined from the 47 DFT simulations at different strain states, using the unstrained simulation as the reference state for the calculation of both the strains and ∆ω

^{m}. Equation (1) shows that the values of $\gamma 1m$ and $\gamma 3m$ calculated in this way depend on both the observed ∆ω

^{m}and the $\omega 0m$ of the reference state (Angel

*et al.*, 2019; Murri

*et al.*, 2019). Because $\gamma 1m$ and $\gamma 3m$ will be used to calculate changes in phonon-mode wavenumbers in zircon relative to room conditions, and not to 0 K, their values (Table 2) have been calculated not with $\omega 0,DFTm$ but with the experimentally observed wavenumber at room conditions, $\omega 0,298Km$.

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 $\gamma 1m$ and $\gamma 3m$. The majority of modes have positive values for both $\gamma 1m$ and $\gamma 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 $\gamma 1m$ and $\gamma 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_{4} tetrahedra.

In contrast, the *B*_{1u} mode near 131 cm^{−1} (which is neither Raman nor infra-red active) corresponds to a rigid rotation of the SiO_{4} tetrahedra around the **c**-axis, and has very large negative values ~ −6 for $\gamma 1131$ and $\gamma 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^{−1} and the *E*_{g} mode near 201 cm^{−1} 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_{4} tetrahedra around the **c**-axis as does the B_{1u} soft mode. Similarly, the *E*_{g} mode at 224 cm^{−1} also involves rigid rotations of the SiO_{4} tetrahedra (Table 2) and has $\gamma 1224<0$ and $\gamma 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 $\gamma im$ from soft modes that drive displacive phase transitions to low-symmetry phases at high pressure (large negative $\gamma im$), and the modes that are affected by, or support, the same phase transitions (small negative $\gamma 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 $\gamma 1m$ and $\gamma 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^{−1} and 224 cm^{−1} 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^{−1} 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 $\gamma 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^{−1}.

## Discussion

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.

^{1}If the ratio of the linear thermal expansion coefficient $\alpha i=1li\xb7dlidT$ to the compressibility is the same for all directions

*i*in a crystal, then in uniaxial crystals such as zircon and quartz:

*et al.*, 2019) which also follows from Equation (1):

Because we have eliminated the phonon-mode Grüneisen tensor components $\gamma 1m$ and $\gamma 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.

## Conclusions

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 $\gamma 1m$ and $\gamma 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^{−1}*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^{−1} and the peak near 439 cm^{−1} 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 $\omega 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^{−1} are stronger (Fig. 6b) and exhibit larger shifts. The difference in the ratios $\gamma 1m/\gamma 3m$ for these three modes means that their isoshift lines have different slopes when plotted against the strains ε_{1} and ε_{3} (*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 $\u2206\omega expm$= $\omega expm-\omega 0,expm$. Strains were then calculated from the values of $\u2206\omega 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 $\gamma 1m$ and $\gamma 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 ε_{1} and ε_{3}. 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*.

## Acknowledgements

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.