Phonon mechanism for the negative thermal expansion of zirconium tungstate, ZrW 2 O 8

Negative thermal expansion (NTE) in ZrW 2 O 8 was investigated using a flexibility analysis of ab initio phonons. It was shown that no previously proposed mechanism adequately describes the atomic-scale origin of NTE in this material. Instead it was found that the NTE in ZrW 2 O 8 is driven, not by a single mechanism, but by wide bands of phonons that resemble vibrations of near-rigid WO 4 units and Zr–O bonds at low frequency, with deformation of O–W–O and O–Zr–O bond angles steadily increasing with increasing NTE-phonon frequency. It is asserted that this phenomenon is likely to provide a more accurate explanation for NTE in many complex systems not yet studied.

In this paper we offer an advance in our understanding of NTE in ZrW 2 O 8 using an approach we have successfully applied to a number of other NTE materials, including ScF 3 , 37 Zn(CN) 2 , 38 Cu 2 O, 39  40 and Y 2 W 3 O 12 . 41 The idea is to model the dynamics of ZrW 2 O 8 in terms of different flexibility models, and then to map the vibrational modes of these models onto the vibrational spectrum of the real material computed using accurate ab initio methods. These are then compared with the contribution of each mode to the NTE, and from this comparison we are able to determine how NTE arises from specific flexibilities in the material. In some cases the NTE can be mapped against very specific phonon modes and specific types of flexibility, but in the case of ZrW 2 O 8 discussed here we are able to show that its NTE is caused not by a single mechanism but, rather, by a much wider range of phonons.
Thermal expansion of materials is generally understood in terms of the Grüneisen model. [42][43][44] Normally, the anharmonicity of atomic bonding causes phonon frequencies to decrease with an increase in volume, due to a corresponding decrease in the effective interatomic force constants as neighbouring atoms are pulled slightly apart. This leads to the normal positive thermal expansion, as the increase in energy due to increasing the crystal volume is offset by an increase in entropy through the resultant decrease in phonon frequencies. However, in the case of NTE materials, we expect to find a significant number of phonons whose frequencies increase with an increase in volume. Given that the thermodynamic weighting of any vibration typically scales as the inverse of the square of its frequency, we expect to find that the most influential phonons for NTE will be those that exist at lower frequencies.
Quite why a particular vibration might have a frequency that decreases on decreasing volume is often understood in terms of the 'tension effect'. 2,3 In the case of ZrW 2 O 8 , for which this was first articulated, 6 the crystal structure -illustrated in Fig. 1 -shows nearly linear Zr-O-W linkages. Since it is envisaged that the Zr-O and W-O bonds are relatively rigid, it follows that transverse vibrations of an oxygen atom in a Zr-O-W linkage might pull the structure inwards, thereby giving rise to NTE. 6,28,34 The question that this idea poses is whether there are specific phonons generating this motion with sufficient thermodynamic weighting. In this regard, the discussion of the tension effect in the simpler material ScF 3 is interesting, where the specific phonons can easily be identified, 45 but the issue of why there is sufficient thermodynamic weighting hinges on the bending flexibility of the ScF 6 octahedra. 37,46 In the case of ZrW 2 O 8 , experimental and computational evidence has led to several different proposals for the realisation of the tension effect and the mechanism of NTE mechanism. The Rigid Unit Mode (RUM) model 29,30,[47][48][49][50] suggests that the low frequencies can be accomplished by normal modes that show minimal deformation of coordination polyhedra, namely the ZrO 6 octahedra and WO 4 tetrahedra. [18][19][20]29,30 An alternative model proposed that the ZrÁ Á ÁW distance remains unchanged, but significant deformation in one or more of the ZrO 6 and WO 4 tetrahedra may be allowed. [21][22][23]31 A further possibility was that adjacent WO 4 tetrahedra may have some additional interaction between each other 33 which would, in turn, mean that both tetrahedra may move as a single unit for low frequency NTE phonons; this model is illustrated in Fig. 2. However, thus far, no proposed mechanism has been comprehensively tested against the overall phonon spectrum in this material-an omission which we now rectify.
In the next section we describe our methodology in detail, namely the methods of calculation, the analysis of NTE, and the analysis of structure flexibility. In Section 3 we examine the phonon spectrum to identify which modes contribute to NTE. Then in Section 2.3 we interpret the phonon spectrum in terms of 16 different models of network flexibility. The analysis is then assessed in the Discussion (5) and Conclusions (6) sections. Finally, the analysis of the flexibility models is presented in an Appendix.

Simulation method
In this work we obtained the phonon spectrum of ZrW 2 O 8 using plane-wave Density Functional Theory calculations using the CASTEP software. 51,52 The Perdew-Burke-Ernzehof (PBE) functional 53,54 was used with optimised norm-conserving pseudo-potentials. 55,56 A plane-wave cutoff energy of 750 eV was used, and a Monkhorst-Pack grid 57 of size 3 Â 3 Â 3 was used to carry out reciprocal-space integration of electronic states.
The unit cell was relaxed at a set pressure of 0 GPa to give a lattice parameter of a = 9.26705 Å. Full details of the crystal structure are given in Table 1, comparing the optimised and experimental 6 (room temperature) crystal structures. The optimised structure has a lattice parameter that is only 1.2% larger than the experimental structure, which is well within the usual tolerance of DFT with the PBE functional. The nearestneighbour W-O and Zr-O distances, given in Table 2, agree to around 0.5% in almost all cases, with one Zr-O distance differing by 2.3%. This too is within the usual tolerance.
Phonon frequencies were calculated using Density Functional Perturbation Theory. 58 A Monkhorst-Pack grid of phonon wave vectors of size 3 Â 3 Â 3 was used to compute the dynamical matrix, and Fourier interpolation of the resulting Fig. 1 The crystal structure of the ordered phase of ZrW 2 O 8 . 6 (a) The structure in ball-and-stick representation. Green spheres are Zr atoms, grey spheres are W atoms, and red spheres are O atoms. (b) The same structure represented as green ZrO 6 octahedra and grey WO 4 tetrahedra. The space group of this material is P2 1 3. Fig. 2 Two interpretations of the interaction between adjacent WO 4 tetrahedra in a pair. The blue atom is the under-bonded O that exists between the two W. (a) Both W are tetrahedrally coordinated and there is only weak interaction between the two units (b) the under-bonded O that exists between the two tetrahedra is part of a stronger bond with a second tetrahedron. One W has standard tetrahedral coordination and the other has '4+1' coordination. Table 1 Comparison between the crystal structures of ZrW 2 O 8 obtained in this study by DFT calculations and from an experimental study using neutron powder diffraction at 293 K. 6 The difference between the two values of the lattice parameter are within expectations for a DFT calculation using the PBE functional. Note that the crystal symmetry, space group P2 1 3, fixes the constraints x = y = z for some atoms  dynamical matrices was used to calculate frequencies for any point within the Brillouin zone. Calculated phonon dispersion curves for the frequency range 0-12 THz are shown in Fig. 3, together with the corresponding phonon density of states. The phonon spectrum shows a gap that extends to 21.5 THz, above which it shows a set of sharp bands of frequencies corresponding to W-O and Zr-O stretch motions up to around 31 THz (Fig. 4). Phonon frequencies were calculated for wave vectors along the main symmetry directions in reciprocal space for plotting dispersion curves. Frequencies were also calculated for 757 wave vectors distributed randomly in reciprocal space for the calculation of the phonon density of states, giving nearly 100 000 independent frequency values (Fig. 4).

Methodology for analysis of phonon contributions to thermal expansion
Here we recall that the Grüneisen theory of thermal expansion [42][43][44] leads to a simple expression for the volumetric thermal expansivity, a V = C V g/BV, where C V is the constant-volume heat capacity, and B is the bulk modulus, V is the volume. The fourth quantity g is the weighted average of the individual mode Grüneisen parameters g i,k , Here i labels the phonon branch and k labels the wave vector, and o i,k is the corresponding angular frequency. The D operators indicate a finite difference obtained from performing calculations at two slightly different volumes. Thus a second set of frequencies was calculated for the ZrW 2 O 8 structure relaxed at a fixed off-equilibrium cell parameter of a = 9.2905 Å. It can be seen from Fig. 3 that for any wave vector there are many modes of similar frequency. The exact matching of the same modes to compute the set of values of Do i,k was performed using the matching of mode eigenvectors with custom software described previously. 41

Flexibility modelling
The existence of a tension effect acting as the mechanism for NTE hinges on the existence of some degree of rigidity in the structure. This may be as simple as single bonds, as advocated in some other works, 34,59 but, as we argued in the case of Cu 2 O, 39 this can be too simplistic because it effectively assigns a tension effect to far too many phonons, including those that have clearly positive values of their mode Grüneisen parameters. The reason why such a mode may have positive thermal expansion is often associated with correlations of the motion of a single bond with the motions of other bonds, typically through the effects of bond-bending interactions. Thus neglecting such correlations prevents a clear understanding of how the tension effect operates to give NTE. For the tension effect to operate, the normal modes must be able to allow easy flexing of the structure whilst taking account of correlations across the crystal, with corresponding low frequency. Thus here we have set up a number of flexibility models, following a similar approach used in previous studies, 39,41 and which we now describe. The ZrO 6 octahedra and WO 4 tetrahedra were modelled in two ways. The first was to have rigid bonds (Zr-O or W-O) but flexible bond angles, and the second was to create rigid polyhedra through addition of bond-bending forces within the polyhedra. The two models for the ZrO 6 octahedra are illustrated in Fig. 5, and the models for the WO 4 tetrahedra are shown in the upper part of Fig. 6.
Models were also developed to take account of possible bonding between adjacent WO 4 units, as shown in Fig. 2. Bonding between adjacent rigid WO 4 tetrahedra was modelled in three different ways. First, trivially, was to have no bonding. Second was to have a rigid bond between adjacent rigid WO 4 tetrahedra, with a flexible W-O-W angle. Thirdly was to have a rigid bond between adjacent rigid WO 4 with a rigid W-O-W angle. These last two are shown in the lower pane of Fig. 6.
The two models for the ZrO 6 octahedra coupled with the four models for the WO 4 tetrahedra give a total of 8 models. However, It has also been suggested that there is a rigidity regarding the distances between Zr and W distances. Cao et al. 21,22 proposed the existence of rigid groups consisting of three Zr centres and a Fig. 3 Left: Low energy dispersion curves and densities of states for ZrW 2 O 8 . Right: The same data shaded according to the value of g i,k of each mode at each wave vector. The linear colour scale ranges from bright red (g i,k r À6) to white (g i,k = 0) to bright blue (g i,k Z 6). Bins that make up the density of states are shaded according to the average g i,k for each bin. Wave vector labels follow convention: G denotes wave vector (0,0,0), X denotes ( 1 2 , 0, 0) M denotes ( 1 2 , 1 2 , 0) and R denotes ( 1 2 , 1 2 , 1 2 ). single W centre, forming a tent shape. Thus to each of the 8 models we also added a Zr-W constraint, giving a total of 16 different flexibility models. The rigidity and flexibility for each of these combinations is analysed in the Appendix.
In practice here, bonds and bond angles are described as either effectively rigid or completely flexible, and the models were implemented in the lattice simulation program GULP, 60,61 which is able to calculate equilibrium structures and lattice dynamics. Rigidity within the flexibility models was enforced using harmonic potential energy functions of the form 1 2 k(r À r 0 ) 2 for bonds (r is the distance between neighbouring atoms, and values of r 0 were chosen to correspond to the DFT bond lengths) and 1 2 K(y À y 0 ) 2 for bond angles (y, with ideal values of y 0 for the two types of polyhedra), with large values used for the force constants k and K. † No additional functions, such as Coulomb interactions, were used. Given the manner in which the flexibility models were constructed, any phonon that does not require deformation of any of the designated rigid units of a specific model will be calculated to have zero frequency. Phonon spectra using these simple functions were calculated for each model at the same wave vectors as used in the ab initio calculations using GULP.
The frequencies o and eigenvectors e of the model phonons j were then mapped on to the ab initio phonons i at each wave vector k using the equation where O is a constant. A value of m i,k = 1 indicates that the flexibility model is able to perfectly recreate the atomic motions of the ab initio phonon (mode eigenvector) in question, while a value m i,k = 0 indicates that there is no relation between the two mode eigenvectors. More detailed information on the methods to perform this type of mapping are given in our earlier work. 41 3 Distribution of mode Grüneisen parameters  (1)). The dispersion curves are accompanied by the calculated phonon density of states coloured by the mean value of the mode Grüneisen parameter within each frequency bin (a wider version of the density of states encompassing the high frequency modes is shown in Fig. 4). It can be seen from Fig. 3 that the values of the mode Grüneisen parameters for different phonon branches span the entire Brillouin zone such that the density of states captures a large part of the detail of the phonon spectrum.
The strongest NTE modes exist around 1.2 THz with NTE character weakening slightly as the frequency increases to 2 THz. Bands of weaker NTE modes range from 2-6 THz and from 8-10 THz. There are also some low-frequency, positive thermal expansion (PTE) phonons at wave vectors near G and X at 2.5 THz, and near G and R at 5 THz. In all, there are 48 phonon branches in the 0-6 THz range which mostly drive NTE (with some driving PTE in the upper frequencies); 12 PTE branches in the 6.5-8 THz range; 28 weak NTE branches in the 8.5-10 THz range; 12 PTE branches in the 10-12 THz range; and (not pictured) 32 weak PTE branches above 12 THz.
A similar distribution of NTE and PTE phonon branches was found in Y 2 W 3 O 12 , 41 with the exception of the PTE branches at around 2.5 THz. This is in striking contrast to simpler materials such as ScF 3 , 37,45 Zn(CN) 2 38 and Cu 2 O 39 where the NTE can be seen to be associated with a relatively small number of phonon modes (in some cases, with a set of acoustic modes).  4 Flexibility models Fig. 7 and 8 show the phonon dispersion curves calculated by DFT (same as Fig. 3) shaded according to the corresponding values of m i,k for each of the 16 models. Fig. 7 shows results for the 8 models without the Zr-W constraint (labelled as models a-h), and Fig. 8 shows the corresponding models with the ZrÁ Á ÁW constraint (labelled as models i-p). The dispersion curves in both figures are accompanied by the corresponding phonon density of states, shaded by the average value of m i,k over all modes within each frequency bin. We noted above that the analysis of the flexibility of the different models is discussed in the Appendix, and Table 3 presented there summarises the quantitative flexibility analysis for each case.

Model a: rigid Zr-O rods with rigid W-O rods
This model is the simplest and most flexible. It has m i,k E 1 for nearly the entire low-frequency spectrum, but with some decline for frequencies in the ranges 6-8 and 10-12 THz, indicating a degree of bond stretching for these modes.
Nevertheless, most bond stretching is clearly associated with the higher-frequency modes. In addition, there is some reduction in m i,k values for the weak PTE modes around 5 THz. The contributions to PTE at 5, 6-8 and 10-12 THz are therefore the result of part of the motions involving bond stretching with reduced tension effect. This model, by design, does not give any account of correlated rotations of bonds. However the results of this model show unambiguously that that NTE in ZrW 2 O 8 arises from a tension effect involving nearly-rigid Zr-O and W-O bonds.

Model b: rigid Zr-O rods with rigid independent WO 4 tetrahedra
This model shows values of m i,k E 1 for frequencies lower than 6 THz, but with significant reduced values of m i,k for frequencies above 6 THz compared to model a with rigid W-O rods but flexible O-W-O bond angles. This result, coupled with the result for model a, shows that much of the motion for frequencies below 6 THz involves whole-body rotations of the WO 4 tetrahedra.

Model c: rigid Zr-O rods, rigid and bonded pairs of WO 4 tetrahedra, and flexible inter-tetrahedral W-O-W bond bending
The flexibility of this model covers much of the same part of the phonon spectrum as for model b, but with a slight decrease in overall m i,k values. One interesting feature is that m i,k values fall significantly for phonons around 2.5 THz at the specific wave vectors where the same modes drive PTE rather than NTE. This suggests, firstly, that the distance between two WO 4 units in a pair does not undergo significant change at low energy; secondly, that any change in distance that does occur is not a high-energy process; and, thirdly, that modes involving nontrivial change in the W-W distance contribute to PTE.

Model d: rigid Zr-O rods, rigid and bonded pairs of WO 4 tetrahedra, and no inter-tetrahedral W-O-W bond bending
This model covers much of the same part of the phonon spectrum as models b and c, but with another slight decrease in overall m i,k values. Additionally, there is a further slight decrease in m i,k values around 4 THz, indicating that W-O-W bond-angle bending occurs to some extent around this frequency. This motion, however, does not lead to any drastic change in the NTE contribution of the phonons in that portion of the full vibrational spectrum. It appears that W-O-W bond angle bending does not occur to a significant extent in the low-frequency phonon spectrum of ZrW 2 O 8 .
Models e-h, which now follow, have the same constraints for the WO 4 tetrahedra as for models a-d, but have rigid ZrO 6 octahedra rather than rigid Zr-O bonds and flexible O-Zr-O angles.

Model e: rigid ZrO 6 octahedra with rigid W-O rods
This model has small, non-zero values of m i,k for the entire low frequency phonon spectrum, with m i,k values largest in the 0-2.5 THz region. This shows that O-Zr-O angle bending is a relatively low-energy process and thus occurs to some extent throughout the phonon spectrum. The energy cost associated with bending of O-Zr-O bonds means that this deformation is smaller for the phonons with lower frequencies.
4.6 Model f: rigid ZrO 6 octahedra with rigid WO 4 tetrahedra; rigid unit modes This is the traditional RUM model, which was initially explored for insights into the flexibility of the ZrW 2 O 8 structure to support the tension effect. 29 The flexibility analysis, given in Table 3 in the Appendix, shows that there is an exact balance between the numbers of degrees of freedom and constraints. As a result, the RUM analysis showed a curved surface of RUMs. This model has large m i,k for modes surrounding the G point. However it has very weak values of m i,k for the rest of the phonon spectrum, including the reciprocal space regions where RUMs of this model were previously located. Furthermore, even when specific wave vectors known to harbour RUMs were individually inspected, no significant increase in m i,k was observed. In general, the (small) m i,k values are slightly larger at the lower frequencies, indicating that there is less overall polyhedral deformation at the lowest frequency modes.
The result here is consistent with an observation about RUMs in relatively complex systems. 50 ZrW 2 O 8 has 132 normal modes per wave vector, and relatively low symmetry (compared to other cubic space groups). This means that there are many modes of the same symmetry and similar frequency, which allows for mixing of eigenvectors between such modes -this point has been discussed recently. 49 Therefore, although model b shows that WO 4 tetrahedra often move as rigid units, and model e shows that ZrO 6 octahedra deform less in the eigenvectors of modes with lower frequencies, the combination of both types of rigid-unit motion are both limited as shown by the RUM analysis, and are diluted by eigenvector mixing.

Models g, h: rigid ZrO 6 octahedra with rigid, bonded WO 4 tetrahedra without (g) and with (h) inter-tetrahedral W-O-W bond bending
Model f represents the cross-over between flexible and constrained models, as shown in the flexibility analysis of Table 3 in the Appendix. Models g and h are shown by the same analysis to be over-constrained, and this is confirmed by the extremely low of zero values of m i,k E 0 for all phonons. The results for these two final cases in this suite of models without constraints on the ZrÁ Á ÁW distances indicate that modes with rigid ZrO 6 and WO 4 units cannot exist without some relative displacement of the WO 4 units within a pair.
Models i-p which now follow use exactly the same sets of constraints used in models a-h respectively, but now with the addition of the ZrÁ Á ÁW constraint implemented as a rigid rod. Fig. 8 shows the results of m i,k mapping for the eight flexibility models described above when a rigid ZrÁ Á ÁW bond, corresponding to the proposed tent mode, 21,22 is incorporated. In cases i-m, namely the four models with rigid Zr-O bonds without O-Zr-O angular constraints, and the model with rigid ZrO 6 octahedra and rigid W-O bonds but no O-W-O angular constraints. In these cases the distributions of values of m i,k over the entire vibrational spectrum replicate what is seen in the corresponding models without ZrÁ Á ÁW bonds, but with a significant overall reduction in the values of m i,k . This outcome conflicts suggests that W-O-Zr bond bending will occur across the range of frequencies, and offers little support for the operation of the tent model.

Models n-p: over-constrained models with ZrÁ Á ÁW bonding
These three models with rigid ZrO 6 octahedra and WO 4 tetrahedra are seen from Table 3 in the Appendix to be overconstrained, and the lack of any flexibility in the vibrations is shown by the near-zero values of m i,k for all phonons in the corresponding panels in Fig. 8.

View Article Online
This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys.

Discussion
Unlike simpler systems such as Zn(CN) 2 , 38 Cu 2 O 39 or ScF 3 , 37 flexibility-model eigenvectors do not neatly map onto a small number of specific modes identified in the ab initio calculations. Instead, a high degree of eigenvector mixing is present, so that a given type of flexibility ends up being spread over multiple ab initio phonons, a phenomenon made possible by the low symmetry of the structure and large number of modes. For example, as noted above, eigenvector mixing is particularly evident for the rigid ZrO 6 and WO 4 model -the RUM modelmapped in Fig. 7f Some NTE modes around 4 THz involve additional bending of the W-O-W angle between two WO 4 units in a pair. However, since there is no effect on the values of the mode Grüneisen parameters, this additional angle bending appears to make no extra contribution to thermal expansion behaviour. Meanwhile the small number of PTE modes in the 0-6 THz range appear to correlate with changes in the distance between two WO 4 units in a pair at 2.5 THz and to a small amount of bond stretching within the coordination polyhedra at 5 THz. From the dispersion curves shown in Fig. 7 it can be seen that both sets of modes (and thus the PTE that they drive) have a strong dependence on wave vector, appearing around G and X, and at G and R points respectively. The fact that the occasional separation of adjacent WO 4 units does not significantly affect mode frequency suggests that there is no strong bonding between them, but that they do move in tandem more often than not, most likely as a consequence of constraints imposed by the wider ZrW 2 O 8 framework.
The 12 PTE phonons in the range 6-8 THz involve further O-Zr-O and O-W-O bond-angle bending, as well as a significant amount of Zr-O and W-O bond stretching. This latter effect is a well-known mechanism for PTE and is what ultimately gives these phonons their PTE character. The 28 weak NTE modes in the range 8-10 THz involve negligible bond stretching but increased bending of both O-Zr-O and O-W-O bond angles. This band of modes resembles the traditional tension effect, but they exist at high frequency and thus only contribute to NTE relatively weakly due to the large amount of energetically-costly bond bending they incur. The 12 phonons in the 10-12 THz range comprise another band of weak PTE modes that involve similar distortions to those seen in the 6-8 THz range but with a much greater degree of bond stretching that gives these modes both their high frequency and their PTE behaviour. PTE modes above 12 THz consist of the remaining bond stretches.
As a result of this detailed analysis, we conclude that the NTE in ZrW 2 O 8 cannot be described by a mechanism that involves a small set of normal modes, in contrast to other NTE materials that have been studied in similar depth. [37][38][39][40] Instead, NTE is driven by a broad spectrum of phonons extending over a wide range of energies. The strongest of these exist around 1.2 THz and involve minimal deformation of the WO 4

Conclusions
Our study of ZrW 2 O 8 presented in this paper, together with the previous work on other NTE in Y 2 W 3 O 12 , 41 strongly suggests that the wide expectation that NTE might always be explained in terms of simple mechanisms will not hold for all framework materials. In the simplest of cases such mechanisms can provide an elegant description of observed behaviour. For instance, the specific modes shown to drive NTE in ScF 3 37,45 and Zn(CN) 2 38 look like RUMs, and in Cu 2 O the important modes are those in which linear O-Cu-O rods move as rigid entitites. 39 This seems logical given that deformation of a coordination polyhedron is an energetically-costly process, and a tension effect in which the important modes have minimal polyhedral deformation will exist at the necessary low frequency. However with added complexity, such as a large number of phonon branches in a low symmetry structure, we find that simple mechanisms are not sufficient to explain the NTE, even though such modes can be supported by the atomic structure of the material. Mode mixing is always present in the vibrational spectra of real materials but, in less-complex high-symmetry structures, these effects are minimal, and simpler mechanisms can be sufficient to explain NTE. It is only through the in-depth study of phonons in NTE framework materials with structural complexity that the more complete picture can be drawn out. Since many NTE materials are also have relatively complex structures, 10 the authors offer the identification of this more complete explanation of NTE in ZrW 2 O 8 as the key contribution of this paper to a general understanding of NTE across all framework materials.

Data and code availability
The main CASTEP code is available for academic use. Information is available from https://www.castep.org/CASTEP/Getting CASTEP. The GULP code is also available for academic use and can be obtained from https://gulp.curtin.edu.au/download. html. The software used for the analysis of the results presented here is available from https://osf.io/bcauq/files/osfsto rage/6425168f63995b07eae19a58.

Conflicts of interest
There are no conflicts to declare.

Appendix: flexibility analysis
We consider the flexibility of rigidity of the models discussion in Section 2.3 in terms of both the Thorpe-Phillips constraint counting approach 62-66 and the approach developed for analysis of polyhedral structures and rigid unit modes, 47-49 both of which have been more recently reviewed in a paper on rigid unit modes. 49 The main task is to count the number of degrees of freedom, F, and the number of constraints C, per formula unit. The difference FÀC gives the number of normal modes per formula unit in which the structure can buckle without the constraints being broken. That is, these normal modes will not stretch bonds or bend constrained bond angles. If F 4 C, the structure is floppy in some ways, and there will be a set of normal modes calculated by diagonalisation of the dynamical matrix to have zero frequency. If C 4 F, the structure is over-constrained and unlikely to have any floppy modes. In the special case F = C, symmetry can make some constraints degenerate, so there will be some floppy modes.
For atoms that are bonded but not constrained by bondbending forces, the model assigns 3 degrees of freedom to each atom, and 1 constraint per bond. Thus in the case were we consider ZrW 2 O 8 to consist of bonds with no angular forces (Section 4.1), the 11 atoms give rise to three degrees of freedom each, and hence F = 33. The Zr-O bonds give 6 constraints, and the W-O bonds from the two WO 4 tetrahedra give 8 constraints, so that in total C = 14.
When we consider rigid polyhedra, it would be possible to add bond-bending constraints (for a polyhedra with r bonds, there are 2r À 3 angular constraints), but we prefer to count each polyhedron as a rigid entity with 6 degrees of freedom (now allowing for rotations). Where one vertex is connected to the vertex of another polyhedron, we assign 3 constraints per shared vertex. Where the vertex is bonded to another cation without angular constraints, the bond itself is the constraint as before. In the case where we have rigid ZrO 6 and WO 4 polyhedra, we have F = 3 Â 6 = 18. Each WO 4 polyhedron has only 3 linked apices to consider. Avoiding double-counting, we count the number of constraints from both set of polyhedra as C = 18/2 + 18/2 = 18. Thus we have the marginal case F = C, which was discussed in the initial analysis of ZrW 2 O 8 in terms of RUMs, 29,30 and as discussed above in Section 4.6.
In Section 4 we consider two variations of the model. In the first, it is recognised that in the crystal structure of ZrW 2 O 8 a pair of neighbouring WO 4 tetrahedra have one oxygen atom on one tetrahedron that is unusually close to the W centre of the other tetrahedron. Thus it is possible to add bonding between this atom and the second tetrahedron, making a WO 5 polyhedron with a shared vertex. With no additional constraints, we can still associate 12 degrees of freedom to the W 2 O 8 unit, but we now have 3 additional constraints at the shared vertex (Fig. 6). The second variant is to add a constraint that prevents bending of the W-O-W bond, which can be achieved by adding a single constraint on the WÁ Á ÁW distance.
It has also been proposed that there are effectively rigid WZr 3 groups of atoms, in what has been called the 'tent' model. 21,22 We take account of this by adding Zr-W bond constraints, with three such connections for each Zr and W atom. In this model, we add 3 constraints per Zr atom.
For a final example of mixing different types of constraint, let us consider the case of Zr-O bonds (no angular constraints within the ZrO 6 octahedra), and rigid WO 4 tetrahedra with the additional W-O bond to form a linkage of WO 5 and WO 4 (Section 4.3). The Zr atom gives 3 degrees of freedom and 6 bond constraints. The W 2 O 8 unit gives 12 degrees of freedom and 3 constraints. Thus in total F = 15 and C = 9. If the W-O-W is made unbendable, we have C = 10. If additionally we have the Table 3 Analysis of flexibility for each of the flexibility models for ZrW 2 O 8 , as defined in Fig. 7 and 8, with the keys a-p defined by the grid in these two figures. F denotes the number of degrees of freedom per formula unit and C denotes the number of constraints. For the degrees of freedom, the subscripts ''atoms'' and ''poly'' indicate the contribution due to individual atoms (3 per atom) or structural polyhedra (ZrO 6 octahedra or WO 4 tetrahedra; 6 per polyhedron) respectively. For the constraints, the subscripts ''bonds'' and ''links'' indicate the contributions from the bonding constraints (1 per bond) and from the atoms shared by two polyhedra (3 per shared atom). Without subscripts, as in the last column, F = F atoms + F poly and C = C bonds + C links . The central horizontal space separates the models shown in Fig. 7 (above the central space, without Zr. . .W constraints) and Fig. 8 (below the central space, with ZrÁ Á ÁW constraints) constraints of the tent model, we then have C = 13. In all cases F 4 C.
Analyses for all 16 cases considered in this paper, namely models a-h as shown in Fig. 7 and models i-o as shown in Fig. 8, are given in Table 3.