We have carried out ab initio hybrid Hartree-Fock/Density Functional Theory simulations to determine the structure and vibrational modes of zircon, ZrSiO4, 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 μm3) 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 I41/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 I41/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 SiO4, 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 SiO4 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 SiO4 tetrahedra and the ZrO8 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 SiO4 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 ZrO8 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 SiO4 tetrahedra and the ZrO8 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 ω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).

Figure 3 shows contour plots of the calculated changes in the wavenumbers of five phonon modes that have been selected to illustrate the range of behaviour found in zircon. The contours are lines of equal change in the wavenumbers of the phonon modes from the unstrained reference state, and in Fig. 3 they are simply interpolated between the calculated wavenumbers of the DFT simulations, without fitting a model. They join strain states that have the same wavenumber for the phonon mode and thus, if they are Raman-active modes, these strained crystals will exhibit the same Raman shift. As a consequence, it is not possible to determine the strain, or the stress, in an inclusion by measuring the position of a single Raman peak. We call the contours “isoshift lines”. All of the phonon modes exhibit isoshift lines that are approximately linear, parallel and equally spaced (Fig. 3). Different modes have different slopes of the isoshift lines with respect to the strain axes; for example the Eg mode near 356 cm−1 clearly responds more to ε1 than ε3, whereas the opposite is true for the Eg 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: 
-ωmω0m=2γ1mε1+γ3mε3,
(1)
in which ω0mis the wavenumber of the phonon mode in the unstrained reference state, and γ1m and γ2m=γ3m are the only two non-zero independent components of the phonon-mode Grüneisen tensor (Angel et al., 2019). The values of γ1m and γ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 γ1m and γ3m calculated in this way depend on both the observed ∆ωm and the ω0m of the reference state (Angel et al., 2019; Murri et al., 2019). Because γ1m and γ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 ω0,DFTm but with the experimentally observed wavenumber at room conditions, ω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 γ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 SiO4 tetrahedra.

In contrast, the B1u mode near 131 cm−1 (which is neither Raman nor infra-red active) corresponds to a rigid rotation of the SiO4 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 B2g mode near 265 cm−1 and the Eg 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 SiO4 tetrahedra around the c-axis as does the B1u soft mode. Similarly, the Eg mode at 224 cm−1 also involves rigid rotations of the SiO4 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 Eg modes at 201 cm−1 and 224 cm−1 that are coupled to the soft B1u 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 γ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 Eg 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.

When the strains in a crystal are induced by a change in hydrostatic pressure at constant temperature, we can rewrite Equation (1) in terms of the axial compressibilities βi=-1li·dlidP:  
ωmPT=2γ1mβ1+γ3mβ3ω0m.
(2)
Hazen & Finger (1982) established that, in general, directions in minerals that are stiff and therefore have large linear moduli and small values of linear compressibility, tend also to be the directions of relatively low linear thermal expansion. This makes intuitive sense in a simplistic way; stiffer directions in a mineral structure are often associated with chains of strong inter-atomic bonds, and the strong bonding might be expected to lead to low values of linear thermal expansion.1 If the ratio of the linear thermal expansion coefficient αi=1li·dlidT to the compressibility is the same for all directions i in a crystal, then in uniaxial crystals such as zircon and quartz: 
α1β1=α3β3=αVβV.
(3)
Within 1% this is true for quartz at room conditions (Table 3). Substitution from (3) into (2) gives: 
ωmPT=2γ1mα1+γ3mα3βVαVω0m.
(4)
Comparison to the expression for the effect of isobaric temperature change on the phonon wavenumbers (Angel et al., 2019) which also follows from Equation (1): 
ωmTP=-2γ1mα1+γ3mα3ω0m,
(5)
shows that, in this particular case of isotropy, the relationship between the isobaric and isothermal changes in phonon wavenumbers is: 
ωmPT=-βVαVωmTP.
(6)

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 αii (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 γ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−1Eg 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 ω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 γ1m/γ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 ωexpm= ωexpm-ω0,expm. Strains were then calculated from the values of ωexpmby 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 ε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.

1
Of course this is an overly simplistic view. Stiffness in a given direction in a crystal structure is due to the mutual repulsion of adjacent atoms, not to their attraction, and it is a property of the static structure. Whereas the thermal expansivity is related to the anharmonicity of the inter-atomic interactions and the strength of the bonding and is therefore an intrinsically dynamic property!
This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Commercial use right is not granted.

Supplementary data