Negative thermal expansion in isotropic crystals

Abstract

The purpose of this article is to take away some of the mystery associated with understanding the phenomenon of negative thermal expansion (NTE) in isotropic crystals. For sure, NTE is counter to our textbook intuition about the normal positive thermal expansion of chemical bonds, but there are now sufficient results from experimental and theoretical investigations to enable us to rationalise the existence of NTE in materials with cubic symmetry, and indeed to be able to predict its possible existence in any new material. In this article, we posit that there are four principles that can enable us to rationalise the origin of NTE in any cubic crystal: the existence of a network structure, the ability of the structure to support low-frequency vibrations, the spreading of the important vibrations into a sufficiently large volume of reciprocal space, and the limiting factor of dynamic disorder. We describe these with regard to a number of key examples. We point out that the same principles will operate for anisotropic crystals, but we offer the caution that NTE in anisotropic materials can also arise from factors missing in isotropic materials. We also point out that NTE arises in the face of competing mechanisms that ordinarily would give positive thermal expansion, so the practical existence of NTE arises as a balance between mechanisms pulling in opposite directions. We therefore argue that one should not point to a single simple mechanism for NTE. Instead, NTE is associated with the four principles together with chemical effects that can tilt how these principles are balanced against the drive of the chemical bonds towards positive thermal expansion. We present a practical way to understand these principles applied to the recently identified NTE material NaZr₂(PO₄)₃. There is a long-standing and well-demonstrated importance of glass ceramics for many technical and domestic applications, which practically achieve many of the goals of the NTE materials community. In the light of this, it is our feeling that the key achievement of much of NTE research over the past 30 years has not so much been in terms of finding alternative NTE materials as in what we have been able to learn from these materials, which in turn has enabled us to elucidate the general principles that underpin our understanding of how NTE materials work.

---

1 Introduction and background

Current interest in negative thermal expansion (NTE) – the phenomenon in which materials shrink rather than expand when heated – started from a study of NTE in the cubic crystal ZrW₂O₈ in 1996.^1,2 It was shown that NTE exists over the complete range of temperatures over which this material is stable and is not adversely affected by the existence of a phase transition. As these authors pointed out, NTE in ZrW₂O₈ had been identified nearly three decades earlier,³ but for some reason it did not elicit much interest at the time. Indeed, prior to 1996 it appears that interest in NTE had been minimal,^4,5 or its existence effectively dismissed.^6,7 The report of NTE in zeolites a few years earlier, predicted by simulations and confirmed by experiments,^8,9 did not appear to have the same impact as the study of ZrW₂O₈.^1,2 Perhaps in 1996 the scientific community was in a better shape than before for a more determined effort to understand NTE. Around this time there was interest in a range of other ‘negative properties’, such as negative refractive index and auxetic behaviour. Furthermore, laboratory equipment was more advanced, as was computational capability for simulation studies. What followed after 1996 was the discovery and characterisation of many new materials showing NTE.

In this article we focus on NTE materials with cubic symmetry. We have the impression that these materials are now better understood than is often appreciated, and our intention here is to present a coherent view of the various factors that can lead to NTE in cubic materials. Indeed, we believe that based on these insights and a few calculations, it ought to be possible to predict whether a given cubic material might show NTE.

It might be thought that crystals with anisotropic symmetry are not so different from those with cubic symmetry, and indeed in the early years after 1996 there was equal interest in a range of anisotropic NTE materials such as the family of materials whose crystal structures are isomorphic with that of orthorhombic Sc₂W₃O₁₂.^10–14 Certainly the mechanisms for NTE in cubic materials that we discuss here apply equally to anisotropic crystals. However, the existence of anisotropy creates two new options for NTE. One is from the effect of elastic anisotropy, which we have discussed in some detail with respect to some specific case studies.^15,16 The second is that NTE can be driven by the existence of a structural phase transition, and we have discussed this in detail for the case of Cu₂P₂O₇.¹⁷ In both cases NTE can arise quite independently of the mechanisms we discuss here for cubic materials. Some comments in this regard are given in Section 8.

There have been many reviews of NTE,^18–31 often from the perspectives of materials and experimental characterisation. In this regard, we particularly recommend the recent comprehensive review of Naike Shi,²⁸ which gives details of many families of NTE materials. We also recommend the Focus article by Coates and Goodwin,³² which discusses how to quantify NTE in isotropic materials, and in so doing provides a good overview of the range of isotropic NTE materials.

Since the publication of the primary reviews of mechanisms from a decade ago,^22,33 the general principles have been articulated for a number of key examples, with the role of simulation becoming particularly important. This paper is not intended to be a review per se, but it is offered as an attempt to understand how the various principles apply to different cases. However, it is useful to start with a brief overview of the different cubic materials known to show NTE before we move on to discuss the underlying general principles, not least because we will use many of these materials as examples of how the principles apply in practice.

Here we focus on cubic crystals whose NTE is driven by a standard phonon mechanism (summarised in Section 3.1). These materials have what is often called a network structure, in which structural polyhedra, such as SiO₄ and WO₄ tetrahedra, or ZrO₆ and ScF₆ octahedra, are connected to each other at the polyhedral vertices to form an infinite connectivity. Because we focus on the phonon mechanism, we explicitly do not consider NTE that arises from magnetic and electronic effects, except to briefly mention the case of SmB₆, because at first sight it looks deceptively like a material that could show NTE for the reasons that we have NTE in the other materials discussed here. NTE arising from magnetic and electronic effects has been reviewed by other authors, which we recommend.^{21,24,29,34,35}

2 Cubic crystals showing NTE

2.1 Notable examples of NTE in cubic crystals

Starting from the rediscovery of NTE in ZrW₂O₈ in 1996, a number of cubic NTE materials have been discovered, and we briefly review them in this section. The first two examples are older, and curiously within families of materials that mostly show positive thermal expansion (PTE) (Sections 2.1.1 and 2.1.2). Some examples show relatively simple crystal structures (Sections 2.1.6 and 2.1.8). The most recent observation of NTE in a new material, to the best of our knowledge, was in 2024 (Section 2.1.12).

The literature variously reports two types of measurement of thermal expansivity, either by volume, α_V = V⁻¹∂V/∂T, or its linear form, α_l = l⁻¹∂l/∂T, where l is typically the unit cell edge length. To a good approximation, α_V ≃ 3α_l. Typical units are 10⁻⁶ K⁻¹, which is usually written as MK⁻¹. We have typically converted from published values of α_l to α_V in this paper when publishing reported measurements of thermal expansivities.

2.1.1 Alkali halides with the rocksalt structure. These have the simplest of crystal structures of the chemical compounds – indeed, a basic textbook structure – which is not exactly a network structure, but which has some shared characteristics. NTE in RbI and RbBr, which share the rocksalt structure, was discovered to exist at temperatures below 10 K half a century ago.³⁶

2.1.2 Zinc blende structures. Many materials with the zinc blende structure show NTE, as summarised in a recent comprehensive study,³⁷ usually at low temperature (typically below ca. 100 K). The corresponding elements, silicon, germanium and tin, also show NTE, at temperatures below 50 K. CuCl is the structure with the largest NTE, with an extreme value of α_V = −25 MK⁻¹.⁵ NTE is seen in the cadmium and zinc chalcogenides and in the III–V zinc blende structures.³⁷

2.1.3 β-Cristobalite. β-Cristobalite has not been considered to be one of the important NTE materials, because in fact NTE only occurs at temperatures above 1400 K,³⁸ and for temperatures below this the thermal expansion is positive.³⁹ However, several important NTE materials have crystal structures that are based on that of β-cristobalite, as we will note below, and so it is important to discuss this material early in this brief review.

The crystal structure of β-cristobalite, shown in Fig. 1, is based on the cubic diamond structure, with the silicon atoms occupying the sites of the diamond structure, and oxygen atoms occupying, on average, the sites half-way between two silicon atoms. These sites suggest a linear Si–O–Si connection, but almost always in silicates this connection has an angle of 145°. It appears that in practice the positions of the oxygens atoms are dynamically disordered as the SiO₄ tetrahedra constantly and dynamically tilt their orientations in order to bend the Si–O–Si linkages.⁴⁰ Such motions are enabled by the spectrum of rigid unit modes (RUMs), namely normal modes of vibration in which the structural polyhedra, in this case the SiO₄ tetrahedra are able to move without distortion, as if they are moving as rigid units.^41–45

Fig. 1 Crystal structure of β-cristobalite (space group Fd3̄m). The oxygen atoms are represented by the red spheres, and the silicon atoms are represented by blue spheres within the accompanying shaded tetrahedra.

2.1.4 Zirconium tungstate, ZrW₂O₈. As we noted above, the curious thing about ZrW₂O₈ is that the existence of NTE over a wide range of temperatures (from 0 to 1030 K) was first reported in a short note, but presented merely as data without comment.³ The NTE in ZrW₂O₈ was rediscovered in 1996, together with the discovery of NTE in some isostructural compounds.^1,2 In fact ZrW₂O₈ undergoes a phase transition involving orientations of pairs of neighbouring WO₄ tetrahedra, but NTE exists in both phases with some slight difference in the magnitude.¹

The crystal structure of ZrW₂O₈ is shown in Fig. 2. It shows a network of connected ZrO₆ octahedra and WO₄ tetrahedra. Each WO₄ tetrahedron has one non-bridging bond.

Fig. 2 Crystal structure of ZrW₂O₈ (space group P2₁3). The oxygen atoms are represented by the red spheres. The zirconium and tungsten atoms are shown as green and blue spheres respectively, at the centres of their shaded coordination polyhedra.

Because of the importance of ZrW₂O₈ there have been many experimental and simulation studies. Briefly we note extensive neutron powder diffraction measurements^46,47 for the study of the average structure, and X-ray absorption spectroscopy and total scattering measurements with X-rays and neutrons^48–53 for measurements of local structure. Atomic vibrational dynamics have been studied by inelastic neutron and light scattering experiments.^54–57 The effects on its thermodynamic properties have been studied by detailed measurements of the heat capacity.^56,58 Many simulation studies focussing on the lattice dynamics have been reported,^59–65 which have been supported by recent measurements of dispersion curves using neutron inelastic scattering from a single crystal.⁶⁶

2.1.5 Zirconium vanadate, ZrV₂O₇. NTE in ZrV₂O₇ was identified in 1988,⁶⁷ and studied again^68,69 around the same time that this group was studying ZrW₂O₈.^1,2 One of the best presentations of the data for ZrV₂O₇ is in the early review by John Evans.¹⁸ ZrV₂O₇ has a similar crystal structure to that of ZrW₂O₈, Fig. 3, but the reduced number of oxygen atoms is accommodated by joining two pairs of VO₄ tetrahedra with a bridging atom. In contrast, the same two tetrahedra in ZrW₂O₈ have non-bridging bonds. But this connectivity leads to linear V–O–V bonds with shortened V–O bonds, which creates some sense of disorder similar to that seen in β-cristobalite (Section 2.1.3) at high temperature, and leads to a phase transition to a structure with a 3 × 3 × 3 supercell.⁷⁰ Replacement of vanadium by phosphorus gradually reduces and eliminates the NTE^67,68

Fig. 3 Crystal structure of Zr(V,P)₂O₇ (space group Pa3̄). The oxygen atoms are represented by the red spheres. The zirconium and vanadium/phosphorus atoms are shown as green and dark grey spheres respectively, at the centres of their shaded coordination polyhedra.

2.1.6 ReO₃, ScF₃, and double cation derivatives. These materials have simple cubic crystal structures, with the cation occupying the site at the corner of the unit cell, and the anions occupying the sites half way along each edge of the unit cell. In this crystal structure, which is shown for ScF₃ in Fig. 4, the cation lies at the centre of an octahedral arrangement of anions, with the octahedra connected at corners. This therefore is essentially a cubic perovskite crystal structure with missing 12-coordinated A-site cations.

Fig. 4 Crystal structure of ScF₃ (space group Pm3̄m). The fluorine atoms are represented by the green spheres, and the scandium atoms are represented by the pink spheres with shaded octahedra.

NTE was discovered in ReO₃ in 2008^71–73 and in ScF₃ two years later.⁷⁴ The NTE in ReO₃ is relatively small and limited to low temperatures (note that an additional higher-temperature NTE has been reported from neutron powder diffraction measurements,⁷² but this was not confirmed by other studies⁷⁵ and is likely to be associated with reduction of the oxygen content). On the other hand, NTE in ScF₃ has a large magnitude and exists up to around 1100 K. Crystal structure analysis of both ReO₃ and ScF₃ clearly shows large amplitudes of the transverse motions of the anions.^72,74

NTE has also been observed in double-cation derivatives, such as CaZrF₆.^76–78 The size of the NTE can be larger than that in ScF₃ in some cases, but with other compositions it can be reduced, and even such materials can show PTE.

There have been several studies of the local structures of these materials using X-ray absorption spectroscopy,⁷⁹ X-ray total scattering,⁷⁹ and neutron total scattering.^75,80,81 Dynamics and related aspects have been studied by inelastic X-ray scattering methods.^82–84 There are also many reported simulation studies of NTE in ReO₃,^71,85 ScF₃,^86–88 and the double cation derivatives,^89,90 although different studies frequently do not give consistent calculations of phonon dispersion curves.

Interestingly NTE is not found in the cubic phase of AlF₃⁹¹ or in any of the cubic perovskite phases.

2.1.7 Samarium hexaboride, SmB₆. NTE in SmB₆ has been observed below 125 K from diffraction studies.^92,93 The crystal structure, which is shown in Fig. 5, consists of SmB₆ octahedra connected by bridging B–B bonds, giving the appearance of a network crystal structure similar to that of Prussian blue materials (see below), with the capability to support the same sort of mechanism as found in other materials highlighted in this brief review. However, phonon calculations showed that NTE probably does not arise from a lattice dynamics origin;⁹⁴ instead X-ray diffraction and total scattering measurements suggested an electronic origin of the NTE.⁹³ For this reason, we will not discuss SmB₆ much in this article, but it is interesting as an example of a crystal structure for which we might imagine, incorrectly, that the origin of NTE is similar to that found in the other materials discussed in this article.

Fig. 5 Crystal structure of SmB₆ (space group Pm3̄m). The samarium atoms are shown as pink spheres, and the boron atoms are shown as the grey spheres. Because the octahedral arrangement of boron atoms does not have a central atom we just show the octahedral arrangement using rods for the nearest-neighbour B–B bonds.

2.1.8 Cu₂O and Ag₂O. NTE has been found to exist up to 240 K in Cu₂O^95–97 and up to 200 K in Ag₂O⁹⁷ (the high-temperature NTE of Ag₂O shows a very unusual effect of suddenly increasing above room temperature, but we suspect based on our own diffraction measurements that − as in the case of ReO₃ noted above − the reported high-temperature enhanced NTE may be associated with reduction of the sample).

The crystal structure of these phases, shown in Fig. 6, can be described as interpenetrating face-centred cubic arrangement of the metal ions and a body-centred cubic arrangement of the oxygen atoms. But more intuitively, this structure looks like two interpenetrating β-cristobalite structures of corner-linked O(Cu,Ag)₄ tetrahedra with linear O–(Cu,Ag)–O bonds. Aspects of the dynamics have been studied by EXAFS^96,97 and X-ray total scattering methods.^98,99 Several simulation studies of the phonon dynamics have shown the importance of the acoustic modes in the mechanism for NTE.^100–102

Fig. 6 Crystal structure of Cu₂O (space group Pn3̄m). The copper atoms are represented by the blue spheres, and the oxygen atoms are represented by the red spheres within shaded tetrahedra.

2.1.9 Zn(CN)₂ and related materials. Zn(CN)₂ and Cd(CN)₂ were shown to have amongst the largest NTE of all materials, spread over a wide range of temperatures.¹⁰³ A structure of the same topology, Si(NCN)₂, also shows NTE, but to a lesser extent.

The crystal structures of Zn(CN)₂ and Si(NCN)₂ are similar to that of Cu₂O, with Zn/Si occupying the tetrahedral site occupied by the oxygen site in Cu₂O, and with the linear CN⁻ and NCN²⁻ anions forming the vertices of the tetrahedra and providing bridges between connected tetrahedra. As for Cu₂O, the crystal consists of two interpenetrating β-cristobalite networks. We show the crystal structure of Si(NCN)₂ in Fig. 7. In the case of Zn(CN)₂ the head-to-tail orientations of the CN⁻ anions are disordered;¹⁰⁴ practically this gives the reason for showing Si(NCN)₂ in Fig. 7.

Fig. 7 Crystal structure of Si(NCN)₂ (space group Pn3̄m). The NCN²⁻ molecular anions are represented by the triplet of spheres (carbon is black and nitrogen is royal blue), and the silicon atoms are represented by the light blue spheres within shaded tetrahedra. The crystal structure of Zn(CN)₂ is identical except that the NCN²⁻ molecular anions are replaced by CN⁻ ions with head-to-tail disorder.

Zn(CN)₂ has been studied by X-ray total scattering experiments.¹⁰⁵ Much of the analysis has focussed on dynamics involving twisting and displacement of the molecular anions, but several simulation studies have suggested that the NTE mostly arises from the acoustic modes,^106,107 exactly as in Cu₂O.

Phillips et al.¹⁰⁸ showed that it is possible to form a single-network form of Cd(CN)₂. These authors showed that without any absorbed molecules, this phase can give a volumetric thermal expansivity as large as α_V = 100.5 MK⁻¹, but with absorbed CCl₄ molecules the NTE decreases monotonically with the occupancy, eventually turning into PTE. The volume expansivity for the double-network Cd(CN)₂ is −61.2 MK⁻¹,¹⁰³ which is itself larger than that for Zn(CN)₂, −50.7 MK⁻¹.¹⁰³

More recently there has been interest in double-cation versions of Zn(CN)₂, with measurements made for materials of general chemical formula MB(CN)₂.^109,110 Two examples, AgB(CN)₂ and NaB(CN)₂, have been shown to have NTE reaching as high as −40 MK^−1 109 and −92 MK^−1 110 respectively.

Although Si(NCN)₂ has an analogous crystal structure to Zn(CN)₂, but without disorder and perhaps more flexibility,¹¹¹ it shows much smaller NTE.¹¹² This is surprising, and we will discuss possible reasons in Section 7.

2.1.10 Prussian blue family of materials, MM′(CN)₆. These materials are to ScF₃, or more precisely to the double-cation version CaMgF₆, what Zn(CN)₂ is to Cu₂O, with the shared vertex atoms for the octahedra replaced by the CN⁻ molecular anion. Their crystal structure is presented in Fig. 8. Unlike in Zn(CN)₂, the molecular anions usually have head-to-tail order, as shown in Fig. 8. There are many materials whose crystal structures are the same as that of Prussian blue, but some of these have a significant amount of water absorbed into the crystal. The extent of NTE varies significantly for different examples, in part coupled to the existence of magnetic order, but many show NTE.^113–116

Fig. 8 Crystal structure of a generic Prussian blue material (space group Pm3̄m). The CN⁻ molecular anions are shown as the pair of atoms (nitrogen represented by grey spheres, carbon by black spheres, with a shared bond represented by the solid rod), and the two types of cations are represented by the blue spheres within shaded octahedra.

2.1.11 Cubic zeolites. A few silicates, apart from β-cristobalite, have crystal structures of cubic symmetry, but the significant exceptions are the members of the large zeolite family. These are low-density silicates with large internal pores and channels. The example of the cubic zeolite Faujasite is shown in Fig. 9.

Fig. 9 Crystal structure of the zeolite Faujasite (space group Fd3̄m). Oxygen atoms are represented by the red spheres, and silicon atoms are represented by the blue spheres within shaded tetrahedra.

As noted at the start of this article, the existence of NTE in some zeolites had been predicted from simulation and confirmed experimentally^8,9 just prior to the surge in interest in NTE following the publication of the work on ZrW₂O₈¹ in 1996. Experimental work has confirmed NTE in a number of other zeolites.^117–119

The zeolite database¹²⁰ gives 22 zeolites with cubic crystal structures. That NTE probably occurs in most cubic zeolites has been suggested by simulation studies.^121,122 Experimentally, NTE has been confirmed in the cubic zeolites LTA^123,124 and FAU.^8,125,126

2.1.12 NaZr₂(PO₄)₃ langbeinite. The langbeinite family of structures are named after the mineral with chemical composition K₂Mg₂(PO₄)₃. The crystal structure is shown in Fig. 10. It shows connected PO₄ and ZrO₆ octahedra, with the alkali cations lying within cavities in the structure.

Fig. 10 Crystal structure of the cubic Langbeinite phase of NaZr₂(PO₄)₃ (space group P2₁3). The oxygen atoms are represented by the red spheres, and the sodium atoms are represented by the blue sphere. The zirconium atoms are represented by the grey spheres within shaded octahedra, and the phosphorus atoms are represented by the green spheres within shaded tetrahedra.

Many of these materials show phase transitions, and in some cases NTE in the low-temperature anisotropic phase has been reported, presumably reflecting spontaneous strains accompanied by the displacive distortions. Data for thermal expansion of the cubic phases are sparse, but some data have been collated recently,¹²⁷ including older reports of NTE in (NH₄)₂Cd₂(SO₄)₃ and (NH₄)₂Mg₂(SO₄)₃.¹²⁸ NTE has also recently been reported in the cubic langbeinite NaZr₂(PO₄)₃.¹²⁹

2.1.13 Cubic metal–organic framework materials. Amongst the many metal–organic framework materials, cubic MOF-5 (IRMOF-1; Zn₄O(BDC)₃) is perhaps one of the better studied examples. Its crystal structure is shown in Fig. 11. The lattice sites have a central OZn₄ tetrahedron, and the Zn atoms are connected to an oxygen atom from each of the three BDC, C₆H₄(CO₂)₂, linkage organic moieties, completing the formation of linked ZnO₄ tetrahedra.

Fig. 11 Crystal structure of the metal–organic framework MOF-5, OZn₄(C₆H₄(CO₂)₂)₃ (space group Fm3̄m). The zinc atoms are represented by light-blue spheres, the oxygen atoms by red spheres, the carbon atoms by black spheres, and the hydrogen atoms by pink spheres.

MOF-5 shows a significant NTE,^130–133 although in practice relatively few MOFs show NTE. Other examples include UiO-66-type structures,^134,135 which have larger secondary building units (the central groups of atoms) than MOF-5; CUB-5 and 3DL-MOF-1,¹³⁶ which are analogues of MOF-5 with the same secondary building unit but with different linkage organic moieties; and the distinctively different Cu₃(BTC)₂ (HKUST-1).^137,138 Curiously, perhaps, the cubic ZIF-8, which has the same network topology as the zeolite sodalite, does not appear to show NTE, even with favourable doping.¹³⁹

We remark that metal–organic framework materials differ from cyanide materials in allowing the flexibility of the linker molecular ions and of the secondary building units to contribute to the mechanism for NTE, whereas the cyanide anions are effectively rigid over the relevant energy scales for NTE, and they do not form secondary building units beyond the standard tetrahedra and octahedra. Because of this we will not include metal–organic framework materials beyond the simpler cyanides in this article; they probably warrant a separate review article.

2.2 Are there new cubic NTE materials to be discovered?

As far as we are aware, the synoptic review of NTE materials with cubic crystal symmetry given above is almost comprehensive. This leads to the obvious question: are there new NTE materials of cubic symmetry yet to be discovered? For sure, substitution of the molecular anions in some materials, such as demonstrated in the cases OCu₂ → Zn(CN)₂ → Si(NCN)₂, or ScF₃ → MgZrF₆ → ScCo(CN)₆, can give possibilities of new analogue materials. Similarly there are possibilities to discover new structures with open zeolitic networks of tetrahedra. But are new families of materials possible?

Let us begin this discussion by noting that there are whole families of materials – which here we restrict to oxides for convenience – which appear not to show NTE. These include materials with the fluorite (XO₂) or related pyrochlore (X₂X₂′O₇) structures, and with the spinel (X₂X′O₄) and garnet (Y₃X₂X₃′O₁₂) structures, where we denote ions with tetrahedral or octahedral coordination by X (second ones by X′) and other cations by Y. Even cubic perovskite phases (YXO₃) do not appear to show NTE, in contrast to the analogous ReO₃ and ScF₃.

So are there other possible combinations of ionic oxides that might form as cubic structures? Let us consider the case of fully-linked polyhedra. We denote the tetrahedrally coordinated cation as T and the octahedrally coordinated cation as M. Thus we are searching for oxides of chemical formula M_mT_nO_3m+2n. For oxides, since the oxygen ion is divalent, the anion components carry a total charge of −2(3m + 2n). This will require high oxidation states on the cations. In practical terms, the tetrahedral ions must often be a high-valence p- or d-block species (e.g., Si⁴⁺, P⁵⁺, Mo⁶⁺, and W⁶⁺) so that reasonable combinations of M and T can achieve overall neutrality without interstitial compensating cations. One combination is Sc₂W₃O₁₂ (m = 2, n = 3), where Sc³⁺ and W⁶⁺ exactly balance the oxide framework to give a fully connected, interstitial-free network. However, this forms crystals of orthorhombic symmetry.

It is of course possible to create network structures in which some of the cation charge balance is achieved by some cations occupying large cavity sites. The perovskite family provides a clear example of this, as do members of the leucite (KAlSi₂O₆) family in their high-temperature phases and the garnet families. None of these families have yet been shown to have NTE in their cubic phases. However, NTE has recently been discovered in the langbeinite (K₂Mg₂(SO₄)₃) family, as discussed above.¹²⁹ The few other combinations we can imagine also tend to form non-cubic structures, such as NASICON (Na_1+xZr₂P_3−xSi_xO₁₂).

The charge balance is easier for halide anions, because the required cation charge is reduced by half. However, whilst charge balance is easier, finding tetrahedral coordination for metal cations in fluorides is difficult. Most transition metals prefer octahedral geometry in fluorides. BeF₄, BF₄ and ZnF₄ groups are exceptions, but linking them with octahedra to form a 3D network is geometrically rare compared to oxides.

It is not the purpose here to discuss potential new materials with cubic crystal structures that show NTE, but we do wonder whether the likelihood is low. Put another way, have we now discovered all the major families of cubic NTE materials? We of course hope not!

3 Understanding thermal expansion from Grüneisen theory

3.1 Review of basic theory

The key result from classical thermodynamics is that the equilibrium state is that for which the free energy of a system is at a minimum. The free energy itself is a balance between the internal energy E (or enthalpy H) and entropy S, defined as F = E − TS (or G = H − TS), where T is the temperature. A change in the structure that increases the bonding energy, but which also increases the entropy, can be favourable depending on the value of temperature.

The basis of the Grüneisen theory of thermal expansion^140–143 is that a lattice strain, , will increase the energy of the crystal, typically as the square of strain to lowest order, but will also cause a change in the frequencies of the phonons, which will affect entropy. From symmetry we expect that the change in frequency will, to the lowest order, vary as . For cubic crystals, we can define the volume strain , where V₀ is the equilibrium volume.

The total free energy is the sum of the vibrational and elastic free energies, denoted as F_ω and respectively.

The vibrational free energy depends only on the vibrational angular frequency ω, and for a set of modes for unit volume it is given aswhere the sum is over all normal modes, labelled λ, associated with the atoms within the unit volume V₀.

The elastic free energy, per unit volume, is given as

where K is the bulk modulus, defined asThe bulk modulus can be written in terms of elastic constants C_ij and for crystals of cubic symmetry, K = (C₁₁ + 2C₁₂)/3.

The equilibrium condition is . We need to form the derivatives of the component parts. For the vibrational free energy, we have

where n(ω_λ, T) is the Bose–Einstein distribution for phonons of angular frequency, ω, E_λ is the energy associated with a particular frequency given asand the mode Grüneisen parameter (for isotropic systems) is defined asThe last form is a trivial rewriting of the derivative, but is important later when we consider ω_λ² as the solution of the dynamical matrix, in preference to considering the linear frequencies ω_λ.

For the elastic free energy we have

The equilibrium condition, , gives an equation for the thermal strain:

The thermal expansion follows immediately aswhere c_λ is the contribution that the mode labelled λ makes to the overall heat capacity C_V. We can now define an overall Grüneisen parameter obtained as the weighted average of the mode Grüneisen parameters:It follows that we can write eqn (10) in familiar form:Note that this equation only applies to crystals with cubic symmetry. For all anisotropic crystals, the anisotropic form of the Grüneisen theory must be used, which is written in terms of linear rather than volumetric strain. This point is discussed in Section 8.

The mode Grüneisen parameters γ_λ and their weighting coefficients c_λ, and hence the overall Grüneisen parameter [small gamma, Greek, macron], can be readily calculated in a lattice dynamics calculation to get insights into mechanisms of NTE in specific systems. The assumption is made that the harmonic approximation applies for all values of strain and that the effects of anharmonicity are represented by changes in the force constants in the dynamical matrix caused by the strain. This approach is often called the quasiharmonic approximation (QHA). Direct phonon–phonon anharmonic interactions, whose effect increases with temperature, are neglected in the QHA.

We now want to make a comment about why positive thermal expansivity is the normal property of most materials. That is, from eqn (12), the lattice dynamics of a material usually gives positive values of the overall Grüneisen parameter, eqn (11). The theory of harmonic lattice dynamics – which we outline in a lightweight manner in the Appendix for readers who are less familiar with it – is directly dependent on the force constants between pairs of atoms. If two atoms are separated by distance r, the force constants, and hence vibrational frequencies, associated with any atomic displacements are proportional to ∂²E(r)/∂r², where E(r) is the energy between two atoms separated by distance r.

If the interactions between atoms are pairwise and dependent only on distances, let us assume that the energy between two atoms follows as E(r) ∝ ±r^−k, where k ∼ 12 with a positive coefficient for nearest-neighbour repulsive interactions, k = 6 with a negative coefficient for dispersion interactions, and k = 1 with coefficients of both signs for Coulomb interactions. The second-order derivative will scale as ϕ = ∂²E/∂r²∝ ±r^−(k+2). In many crystals, where some of the atoms are on symmetric sites, with their fractional coordinates fixed by symmetry, we can anticipate that most interaction distances scale as because the atoms do not have much freedom to rearrange their relative positions. In this case, the derivative of ϕ with respect to strain will be equal to . That is, the values of |ϕ| will decrease with increasing strain, thereby systematically lowering the phonon frequencies and increasing entropy, and hence favouring positive thermal expansion.

As we see in the Appendix, the correct treatment of lattice dynamics is to consider not r as the key variable, but the displacements of the atoms. However, we show in the Appendix that the derivatives of the energy with respect to atomic displacement end up requiring an evaluation of ϕ = ∂²E/∂r², and so the argument presented here is valid.

Where this argument gives way to NTE, however, is in the special case when the displacement is exactly perpendicular to the direction connecting two atoms. This is the situation when there is a maximum tension effect, as we will discuss in Section 3.3. In this case the relevant derivative has a zero value, and the residual force constant will be given a non-zero value by weaker angular forces. This gives the possibility that the overall force constant might increase with positive strain, and hence with the required negative values of mode Grüneisen parameters to give NTE. We will illustrate this point in Section 3.5.

3.2 NTE arising as a balance between positive and negative factors

It is tempting to look for a specific ‘explanation’ for the origin of NTE in any system, but the question is somewhat nuanced because many materials show many phonon modes with positive values of the mode Grüneisen parameters and some with negative values. The result of NTE or PTE, or NTE over some range of temperature, is more a question of balance. That is, if we have NTE, the question to ask is why there is a sufficient number of modes with negative values of the mode Grüneisen parameters to dominate over the modes with positive values?

Higher-frequency bond-stretching modes always show positive values of mode Grüneisen parameters and measurements of a few instantaneous bond lengths using total scattering or EXAFS always show that the bonds display positive thermal expansion.^81,144,145 Why this is so can be explored using the Morse potential as a representation of the energy of a chemical bond:

E(r) = D[exp(−2α(r − r₀)) − 2 exp(−α(r − r₀))] (13)

where D is equal to the energy of the bond when the atoms are separated by distance r₀, which corresponds to the minimum of the energy. α is a parameter that controls the curvature of the energy around r = r₀. The mean bond length at non-zero temperature T can be derived,¹⁴⁶ within the classical approximation and at lower temperatures, aswhich immediately gives us a positive expansion of the chemical bond. This is in keeping with the conventional interpretation of the positive thermal expansion of a chemical bond in terms of a potential energy function that rises faster with |r − r₀| than does the corresponding harmonic function for r < r₀ and slower than the corresponding harmonic function for r > r₀. We therefore will expect to find that the normal modes of higher frequencies, reflecting stretching of bonds, will nearly always contribute positive values of the mode Grüneisen parameters.

Let us illustrate the discussion about the balance between contributions to positive and negative thermal expansion by considering the case of ScF₃. Its crystal structure is shown in Fig. 12a, where we show the three-dimensional linear Sc–F–Sc connectivities. This figure represents the atoms by their thermal ellipsoids⁷⁴ and demonstrates the extent of the transverse motions of the fluorine atoms. The distribution of the values of individual mode Grüneisen parameters with frequency as calculated by recent density functional theory (DFT) calculations⁸⁸ is shown as a two-dimensional colour map in Fig. 12b. This distribution is not untypical for NTE materials, in that the positive values of the mode Grüneisen parameters are uniformly for the higher frequencies – in this case, for frequencies from 7 THz up to the maximum frequency of 19 THz – with the most negative values of the mode Grüneisen parameters for lower frequencies, beginning at 4 THz and with increasingly negative values down to 1 THz.

Fig. 12 (a) Crystal structure of cubic ScF₃ (space group Pm3̄m),⁷⁴ shown in an alternative representation to that shown previously in Fig. 4. The light blue spheres represent the scandium atoms, and the darker blue ellipsoids represent the fluorine atoms, with the ellipsoids representing the relative extent of thermal motion in three-dimensional space. This representation highlights the significant amplitudes of vibrations of the fluorine atoms in the directions perpendicular to the Sc–F–Sc linkages (transverse vibrations). (b) Distribution of the values of the mode Grüneisen parameters with frequency calculated by DFT methods,⁸⁸ shown as a colour map. Pink represents zero values, and yellow represents maxima in the distribution. (c) Accumulated overall Grüneisen parameter at high temperature, .⁸⁸

It is pertinent to this discussion to explain exactly why we would expect the most negative values of the mode Grüneisen parameters to occur at low frequencies. The intuitive explanation might follow this way. The spectrum of frequencies we see in a crystal for a given wave vector arise by combining all the force constants in different ways, some of which may be large and others small. This of course is what the theory of lattice dynamics does for us – recall that we present a lightweight version of the theory in the Appendix for readers who are less familiar with it – but let us stick to intuition for now. The high-frequency modes will arise from combinations of the force constants that mostly add constructively. On the other hand, the low-frequency modes are likely to involve combinations that give cancellations of large force constants, giving a residual value that is effectively the difference between two large numbers. Instinctively we know that the relative sensitivity of the difference between two large numbers to small changes in these numbers is much higher than that for the sum of two large numbers. Thus when the force constants change due to imposition of a small strain, the lower frequencies in the spectrum are likely to be much more sensitive than the higher frequencies. Thus we intuitively expect that the vibrations of lower frequency will have larger absolute values of their mode Grüneisen parameters.

We can say the same thing more formally from the theory of lattice dynamics. Central to this theory is the quantity called the dynamical matrix, which is described in the Appendix. It is formed from the interatomic harmonic force constants together with relevant phase factors. The important point is that the squares of the mode frequencies, ω², are obtained as the eigenvalues of the dynamical matrix. The dynamical matrix is Hermitian − its complex conjugate is its transpose – and this means that its eigenvalues are real numbers.

The eigenvalues of a Hermitian matrix that are most sensitive to small changes in the values of the matrix elements are those with the lowest values,¹⁴⁷ and this is important when we consider the link between mode Grüneisen parameters and frequencies. Small volume strains will cause small changes in the interatomic-force constants that contribute to the dynamical matrix, and it follows that the largest effects will be on the modes of lowest frequencies. This is the exact mathematical demonstration of the intuitive analysis given above.

This sensitivity of the lower-frequency eigenvalues is then enhanced by the factor of ω_λ⁻² in the equation for the mode Grüneisen parameter, eqn (7). Thus the distribution of negative values of the mode Grüneisen parameters seen in ScF₃, Fig. 12b, at the lower frequencies is typical of most NTE materials for which the same analysis has been performed.

With negative and positive values of the mode Grüneisen parameters occurring over the lower and higher ranges of frequencies respectively, the question is how finely balanced are the two positive and negative contributions. We can evaluate this by performing an accumulated integral of the form

to show how the balance changes as we continuously add higher frequencies to the overall Grüneisen parameter. This integral for ScF₃ is shown in Fig. 12c. The significant contribution to NTE is shown at the lower frequencies, but as the integral is extended to higher frequency the mean is tending towards zero. In the end, the fact that the overall Grüneisen parameter has a negative value is seen to be a result of a very fine balance between positive and negative contributions.

From eqn (11) we see that each contribution is weighted by a factor, c_λ, that depends on frequency and temperature, specifically as their ratio ħω/k_BT. At low temperature only the low-frequency modes are excited, so the values of c_λ preferentially weight the contributions of the low-frequency modes and hence typically of the NTE modes. Upon increasing temperature, the higher-frequency modes are excited, increasing the weighting factor c_λ for these modes, and hence increasing the weighting of the terms with positive values of mode Grüneisen parameters to the overall average. This explains why there are many cubic NTE materials where the NTE is only found at lower temperatures (β-cristobalite is our best-known exception,^38,39 but there are indications of the same in ZrV₂O₇ at temperatures just above the first phase transition upon cooling,¹⁸ and in the isostructural materials ThP₂O₇ and UP₂O₇⁶⁷).

The point is that this shows why asking a question regarding the origin of NTE is very nuanced. One should not point to a given vibration and assign undue significance. For example, it is tempting to point to the RUMs in ScF₃ – the normal modes in which ScF₆ rotates without distortion, which exist with wave vectors on the edges of the Brillouin zone between the special points and ^41,44,45 – and imagine that they are the origin of NTE in this material. Well, they are important, but the question as to the origin of NTE needs to consider not just why modes with negative values of their Grüneisen parameters exist, but why these modes eventually outweigh the large number of vibrations with positive values. In this regard ScF₃ offers us a challenge if we try to understand why it shows NTE: we also need to understand why many other materials with exactly the same features as ScF₃ show only PTE. We will return to this question in Section 9.

3.3 The general representation of the tension effect

The tension effect^19,22 is a simple concept, which we introduce here with reference to two model two-dimensional atomic networks shown in Fig. 13. Both examples show one type of atom, which is at the centre of a structural polyhedron. In the two-dimensional examples, the polyhedron is a square coordination. The polyhedra are connected via shared corners, which in a ceramic (Fig. 13a) is an ion, but in a hybrid material (Fig. 13c) will be a molecular ion. The closest three-dimensional representations of the two cases shown in Fig. 13 are ScF₃ (Fig. 12a) and Prussian blue, both of which are introduced above.

Fig. 13 Simple representations of two-dimensional networks and their atomic motions that can give rise to a tension effect. (a) A simple ceramic, shown as a two-dimensional ScF₃; (b) transverse displacements of the linkage anion; (c) a simple hybrid material shown as a two-dimensional Prussian blue network; (d) transverse displacement and twist motions that might lead to a tension effect. In both networks one type of atom is linked to four neighbouring atoms of the same type via an intermediate atom (a) and (b) or molecule (c) and (d).

Let us consider the ceramic network first. If we denote the polyhedral centre atoms as A, and the linkage atoms as B, the tension effect operates on the linear triplet of A–B–A atoms. The A–B bonds are considered to be relatively stiff, and likely the B–A–B right angles are also quite stiff, but the linear A–B–A angles are relatively flexible. Because of the stiffness of the A–B bonds, transverse vibrations of the B atom may actually pull the two A atoms towards each other rather than stretch the A–B bonds. A geometric argument can show that the contraction of the A–A distance in this case will be proportional to the mean-square-amplitude of the vibration,¹⁴⁸ which will vary linearly with temperature in a classical approximation. This idea actually predated the onset of the current interest in NTE by four decades,^143,149 but was revived with the later work on ZrW₂O₈.¹ It provides a natural foundation for our understanding of NTE, but as we will discuss below (Section 3.6), by itself it does not automatically explain the origin of NTE.

This discussion of the general idea of the tension effect can be extended for the case where the intermediate entity is not a single atom but is a molecular ion, as illustrated in Fig. 13c. For example, we can represent this as A–C–C–A, where C–C represents the molecular anion. In this case, not only can we find a transverse displacement of the C–C molecular anion that will give an analogous contraction of the A–A distance, but a contraction of the A–A distance can also be achieved by rotation of the C–C molecular anion. These two motions are illustrated in Fig. 13d.

We want to add a cautionary note, by way of introducing the next section, that the images in Fig. 13 typically represent optic phonon modes, but often, particularly if there are strong bending forces in the crystal, the tension effect may instead arise mostly from acoustic modes. These may of course include some aspects of the motions shown in Fig. 13b and d, but there are significant differences in how NTE can arise from a tension effect mediated by acoustic and optic modes. Readers should not consider the visualisation of Fig. 13 to represent all of the possible tension-effect mechanisms.

3.4 The tension effect from acoustic phonon modes

The common perception of the tension effect, such as in the examples shown in Fig. 13, is usually that of a representation of the optic modes of the crystal. Moreover, most visualisations are for a zero wave vector, or for a wave vector on a special point on the surface of the Brillouin zone. However, it is equally possible, at least in relatively simple crystal structures, that the acoustic modes can also give a tension effect and thereby make an important contribution to NTE. Indeed, as has been shown from calculations of the lattice dynamics, several NTE materials have their properties primarily determined by the acoustic modes.^102,106,133

At the outset we should comment on the difference between acoustic and optic modes. In the limit that the wave vector q goes to zero, the acoustic modes are the solutions to the dynamical matrix in which the centre of mass is allowed to move, whereas for the optic modes the centre of mass is fixed. Consider a crystal with 4 symmetrically related atoms in the unit cell, that is, that they have equal mass and displacements with equal magnitudes but with the possibility of positive of negative signs. There are four possible linear combinations of displacements in each unit ell that are consistent with the requirements of symmetry. Denoting the displacements by arrows, the four linear combinations are

Mode 1 is the acoustic mode, and at zero wave vector it reflects a uniform displacement of the crystal. Because this pattern of displacement will not involve stretching or bending of bonds, it incurs no cost in energy, and this mode is therefore of zero frequency at a zero wave vector. Put another way, the displacements of atoms in the acoustic mode are completely in phase. On the other hand, for Modes 2 to 4 the displacements are in some out-of-phase combination. Their displacements sum to zero, thereby maintaining a constant centre of mass. These are the optic modes at a zero wave vector. Each vector has three components in three dimensions, so that we have a total of 12 modes, with three being the acoustic modes. This effectively defines the difference between acoustic and optic modes, exactly and only defined in the limit q → 0.

In the limit q → 0, but not exactly q = 0, the transverse acoustic modes correspond to shearing of the unit cell, with displacements perpendicular to the direction of q, whereas the longitudinal acoustic modes correspond to tensile strains parallel to q. Thus for a cubic crystal with q‖[001], the frequencies of the longitudinal acoustic modes are determined by the tensile elastic constant, ω_LA²(q) ∝ C₁₁|q|², whereas the frequencies of the transverse acoustic modes are determined by the shear elastic constant, ω_TA²(q) ∝ C₄₄|q|².

The distinction between acoustic and optic modes becomes less clear upon increasing |q|, and at the Brillouin zone boundary there may be no distinction. The textbook example of the diatomic chain has a clear distinction between the in-phase and out-of-phase motions of the acoustic and optic modes as q → 0, but at the zone boundary one mode has one atom not moving at all, and the other mode has the other atom not moving. Distinctions may also be blurred upon increasing |q| in more complex crystals as the acoustic modes ‘interact’ and anti-cross with optic modes of the same symmetry, leading to a mixing of mode eigenvectors leading to the acoustic branch in the dispersion diagram showing no signature of the original acoustic mode in its eigenvectors.

We show one example of a tension effect operating due to a transverse acoustic mode in Fig. 14. This image shows the effect of a long-wavelength transverse acoustic mode with the wave vector along one of the [110] directions. It leads to a shearing of the unit cell. To minimise deformation of the two sets of angles and to preserve bond lengths, there is a small twist of the cyanide groups and a commensurate rotation of the octahedra. However, for the subtle reason that from unit cell to unit cell the displacements are progressive, meaning that the twist is actually a positive displacement of both atoms, the motions associated with this acoustic mode are formally distinct from those of an optic mode.

Fig. 14 Representation of the deformation of the structure for two-dimensional Prussian blue (shown in Fig. 13c) associated with a long-wavelength transverse acoustic mode whose wave vector is propagating along the diagonal of the diagram.

The motions associated with the transverse acoustic mode are seen in Fig. 14 to act exactly as a tension effect in lowering the volume. We will see below that acoustic modes are responsible for NTE in several cubic materials, including Cu₂O (Section 5.4.4) and Zn(CN)₂ (Section 5.4.5). We will also see below, Section 3.5, that the importance of acoustic modes for NTE arises naturally in a simple model of the tension effect in a one-dimensional crystal.

It might be asked whether we can predict whether acoustic modes are or are not important in any particular system. Our feeling is that acoustic modes are likely to be important in many systems, both isotropic and anisotropic, and the cases where they are not important, such as in ScF₃, are likely to be the exception. The other comment we can make is that in cases where there are many atoms in the unit cell, such as ZrW₂O₈, the tension effect is spread over many phonon modes, acoustic and optic, so the contributions specifically from the acoustic modes are relatively slight. On the other hand, in Cu₂O and Zn(CN)₂, the network flexibility enables the acoustic modes to have low frequencies, and thus these are markedly important compared to the relatively few optic modes that contribute to NTE.

Our point is that we should view the transverse acoustic modes as more likely than not to contribute to NTE and that their relative importance will depend on whether the overall flexibility of the network allows for smaller or larger numbers of optic modes that also contribute to NTE. It is surely the case that the role of acoustic modes has been somewhat neglected in our understanding of NTE.

3.5 Model illustration of the tension effect

We can illustrate the operation of the tension effect by considering a simple one-dimensional model crystal introduced by Fang et al.,¹⁵⁰ shown in Fig. 15. This represents a simple ceramic, perhaps a one-dimensional version of ScF₃. In fact Fang et al.¹⁵⁰ also introduced a model hybrid material, containing rod-like molecular anions, but we will merely show the results of this ceramic model here. The model only considers transverse motions of the atoms, but it does allow for a strain along the axis of the one-dimensional crystal. The model shows large and small atoms, with displacements in the unit cell labeled j represented by U_j and u_j respectively. We consider that the larger atoms are at the centres of a quasi-rigid unit composed also of two smaller atoms, and each smaller atom is the bridge between two linked rigid units.

Fig. 15 One-dimensional model showing NTE, from Fang et al.¹⁵⁰ The model has two types of atoms and three types of potential energy functions. The first, V₁(r) with force constant k, represents the stiffness of the bonds, eqn (17). The second, V₂ with force constant K, represents bending of the bond within the quasi-rigid unit, eqn (18). The third, V₃ with force constant J, represents flexing of two rigid units about their shared vertex, eqn (19).

The model has the minimum number of interactions to create realistic behaviour. The first is a strong radial interaction between nearest-neighbour atoms. Defining the equilibrium distance between atoms as r₀, and allowing for a strain , the distance between any two neighbouring atoms is given as

where i labels the unit cell, and j ∈ {0, 1}. The bond energy is thus given as

The second interaction is a term that reflects the stiffness of quasi-rigid units. In principle this could be an angular term, but it is useful to work with the linearised form in terms of the displacements,

The minimum of this function is such that u_i − U_i = −(u_i+1 − U_i).

The third term is designed to ensure that the RUMs in the model, the vibrational modes in which the quasi-rigid units rotate without distortion, have non-zero frequency, and simply give an energy cost for rotations of the quasi-rigid units. This is expressed in linear form as the transverse flexing of the linear chain, giving an energy cost for nearest-neighbours to move by different amounts,

where we assume that J ≪ k.

Here we give solutions to the dynamical matrix for two special cases. The first case is for very small wave vector q. One solution is for the acoustic mode, whose frequency is linear in the wave vector. The other solution is for the optic mode, which has zero group velocity at q = 0, namely ∂ω/∂q|_q=0 = 0. The two solutions are

The second case is for the wave vector at the boundary of the Brillouin zone, q = 1/2, where again both solutions are standing waves with zero group velocity,

The full dispersion curves for this model, for strain and for a small negative value of the strain, together with the second model representing a hybrid structure, are shown in Fig. 16. Recall that these are the modes for transverse motion, and do not include the stiff longitudinal acoustic modes or the high-frequency transverse optic modes. It is clear that all transverse phonons show a lowering of frequency under negative strain, meaning that they all have negative values of their mode Grüneisen parameters.

Fig. 16 Dispersion curves for the two one-dimensional models described in this text for a given set of values of force constants, together with the density of states (right), from Fang et al.¹⁵⁰ The top diagram is for the ceramic model described in the text, and the lower diagram is for the corresponding hybrid model described by Fang et al.¹⁵⁰ The results are given for strain as the black curves, and for a small negative value of as the grey curves.

The mode Grüneisen parameters can be calculated from . Both q → 0 and q = 1/2 give

The uniform values of the mode Grüneisen parameters across the two branches can be seen in the dispersion curves in Fig. 16. We note four points from these results. First is that the mode Grüneisen parameters are proportional to the bond-stretching stiffness k. That means that stiffer bonds will generate a larger tension effect. This holds even though k does not contribute to the mode frequencies in the strain-free limit. Second is that the mode Grüneisen parameter is lower for the lower-frequency branch, consistent with general results illustrated by the case of ScF₃ discussed above (Section 3.2). Third is that since the frequencies decrease in value with positive stress (compression), this causes the values of the mode Grüneisen parameters to become more negative, such that the NTE increases under compression. This is consistent with the general result for isotropic systems as discussed elsewhere.²² Fourth, the model also demonstrated a feature that we will meet often in this paper, namely the importance of the acoustic modes. In fact the acoustic modes in this model clearly have larger negative values of the mode Grüneisen parameters than the corresponding optic modes.

Fang et al.¹⁵⁰ used the simple model to calculate the variation of the overall Grüneisen parameters, heat capacity and thermal expansivity of both the ceramic and hybrid materials. The results for some trial parameters (those used to calculate the dispersion curves shown in Fig. 16) are shown in Fig. 17 and demonstrated that the tension effect encapsulated in the model gives NTE in the manner found experimentally in many systems (as we will discuss below). Actually the intention of this model was not merely to reproduce NTE, but to demonstrate that a system showing NTE will also show pressure-induced softening.^121,151,152 It should be noted that a simple three-dimensional model has been applied to ScF₃, where we have three interaction terms corresponding to those in eqn (17)–(19).⁸¹ This model also showed NTE and will be discussed below because of some insights it has given.

Fig. 17 Results for the two one-dimensional models, from Fang et al.¹⁵⁰ The two graphs show the variation of heat capacity, C_V/R, overall Grüneisen parameter, [small gamma, Greek, macron], and the thermal expansivity, α, for reasonable values of the model parameters.

3.6 The questions discussed in this paper

We consider that it is not an explanation of NTE to simply cite the tension effect and identify, via calculations, vibrational modes showing negative values of their Grüneisen parameters. We do consider that NTE in most inorganic and hybrid inorganic–organic materials of cubic symmetry is achieved by the tension effect (whilst noting that NTE in some special cases can arise from electronic and magnetic effects²⁴). We also do not wish to consider NTE in crystals with anisotropic symmetry, because these systems have two other possible causes of NTE that are not relevant in cubic systems, namely the role of elastic anisotropy¹⁵ and the origin of NTE from spontaneous strains that accompany structural phase transitions.¹⁷

The tension effect, as we have described it, gives a natural mechanism for the existence of negative mode Grüneisen parameters, as supported by the simple model¹⁵⁰ described in Section 3.5. But the immediate concern we may have is that the three atoms in the A–B–A linkage are never isolated, which means that we need to consider correlations between motions. Consider the simple two-dimensional model in Fig. 13a. The transverse motion of one B atom requires rotation of the corresponding A–B bonds. Each A atom has three other A–B bonds, so an isolated displacement of a single B atom will cause distortionsof two B–A–B right-angles. It is quite possible that this will cost energy, and not a small energy, which will raise the frequency of any collective vibration showing transverse vibrations unless there are correlated displacements of the other atoms within the AB₄ square coordination. What this means is that the energy cost can be avoided if the squares can rotate about the perpendicular direction as if acting as rigid entities, with neighbouring squares rotating in equal and opposite senses. This is in the manner associated with the cubic–tetragonal phase transitions in perovskites, such as SrTiO₃, and is the starting point for the RUM model.^41,45,153 In a three-dimensional cubic perovskite crystal, neighbouring layers can rotate in the same or opposite sense, because the connecting oxygen atom is not displaced in the vibration, and this means that the RUMs exist for all wave vectors along the edges of the Brillouin zone, from to and symmetrically-related directions.

The existence of RUMs can explain why some vibrations have the necessarily low frequencies to give large negative values of mode Grüneisen parameters, but not all systems can support the existence of RUMs. Therefore there may need to be a degree of bond-bending flexibility, but to what extent is not clear a priori. We also need to understand why it is possible for RUMs to give NTE in some systems, such as ScF₃ and zeolites, but not in the other related systems such as the cubic perovskites and the cubic leucites.

The related issue concerns the balance between the negative and positive mode Grüneisen parameters. This balance is illustrated when thinking about RUMs. Consider the simple two-dimensional model shown in Fig. 13a. It has a single RUM, with its wave vector at the corner of the two-dimensional Brillouin zone. The overall Grüneisen parameter is a sum over all vibrational modes, which means an integral over all wave vectors. A single point has a vanishingly small contribution to an integral over the area of the two-dimensional Brillouin zone, and in three dimensions a line of wave vectors also gives a vanishingly small contribution to the integral over the space of the three-dimensional Brillouin zone. Thus the contribution to NTE from the RUMs needs to also be expanded outwards in reciprocal space, but how is this possible?

These are amongst the questions we will now discuss, but they are not the only questions. The comparison between ceramic (atomic linkers between polyhedra) and hybrids (molecular linkers) is useful to explore. The existence of the molecular linkers may give more flexibility,¹⁵⁴ but to what extent? Perhaps greater network flexibility will lead to more modes with negative Grüneisen parameters, and this might explain why NTE in the Zn(CN)₂–Cd(CN)₂ system is relatively large,¹⁰³ but it appears difficult to create larger NTE. Why?

If some new cubic network materials could be found to show NTE, could we have predicted the NTE from the knowledge of the crystal structure, and would we be able to explain the NTE without having to do detailed calculations? Actually we suspect so. After a large body of work that has been carried out during the last few years, although there is not a single general mechanism that we can invoke – well, the tension effect underpins our understanding, but it needs too much else for it to be considered a general mechanism – we do have a set of general principles that are now sufficiently well understood to mean that NTE can be rationalised, in a detailed way, in most systems. In summary, these are

1. We need to have arrangements of atoms giving the possibility of a tension effect; this probably means we have to have a network crystal structure.

2. The vibrational modes giving the tension effect need to be of low frequency; these are the modes most affected by changes in structure. Thus we need to have low-energy flexibility of the network. This may be associated with RUMs, or else with some degree of polyhedral angle flexibility.

3. The modes giving the tension effect need to occupy a sufficient volume of reciprocal space; this may require flexibility of polyhedra in a RUM model.

4. If we have too much flexibility, the thermal motion will lead to dynamic structural disorder, which will limit the extent of NTE. This probably means it is unlikely that a new material will show a much larger extent of NTE than that in the current examples.

In the rest of this article we will articulate these general principles and explain how they are used to understand NTE in the isotropic systems introduced earlier (Section 2.1). We stress, however, that whilst these principles give the possibility of NTE and enable us to understand the mechanism of NTE, they do not automatically give the necessity of NTE occurring in a particular material. We need to always remember that there is a balance between interactions between atoms that will favour PTE, which we highlighted in Section 3.2 and will discuss again in Section 9. Accordingly, we hope that this overview will be able to guide understanding and interpretation of NTE in future studies.

4 The first principles: the existence of a suitable network structure

For a material to show NTE, it is necessary that the crystal structure has arrangements of atoms that give the possibility for the operation of the tension effect, namely a near linear connection of two stiff bonds, and with a relatively low energy cost for bending the linkage angle. Recall our discussion of the tension effect in Section 3.4.

Almost certainly for a crystal to display such an arrangement of atoms means that it will have a network structure, with the central atoms of these near-linear triplets being the linkage atom connecting two structural polyhedra, such as tetrahedra and polyhedra. Examples of network structures are shown in Fig. 1–11. Some examples of the groups of three atoms giving a possible tension effect mechanism are shown in Fig. 18. Not only is the network required to give the connectivities of atoms, but it is necessary for the inward pull of the tension-effect motions to percolate across the crystal; if the relevant groups of atoms were isolated from each other, their motions would not be correlated and therefore could cause contraction over large distances. We see something similar in low-dimensional networks, such as the one-dimensional CuCN¹⁵⁵ and the two-dimensional ZnB₂(CN)₈,¹⁶ where the thermal expansion in the perpendicular directions is positive.

Fig. 18 Examples of polyhedral connectivities in different network structures showing the opportunities for the operation of the tension effect. Some of the bonds have linear connectivities, but others are bent.

None of the examples discussed here and shown in Fig. 1–10 have the structural polyhedra sharing edges or faces. The only example of a possible NTE material where polyhedra share edges is tetragonal ZnF₂, which has been reported to have weak area NTE,¹⁵⁶ but we were unable to reproduce this in our own neutron diffraction measurements.¹⁵⁷ Not all networks need to have fully-connected polyhedra. ZrW₂O₈ is one example that has non-bridging W–O bonds,¹ but this is relatively uncommon in fact. It is a factor that we will discuss in the next section regarding the flexibility of the network. Anticipating the discussion, we would like to remark that the existence of non-bridging bonds within the structural polyhedra does not necessarily mean an increase in the flexibility of a network.

Lattice dynamics calculations on a variety of materials usually show that the bonds in network structures are stiff, giving rise to the highest frequencies.^{85–88,100,111,158} A few experiments focusing on local structure, namely total scattering^{50,51,53,79,81} and X-ray absorption spectroscopy,^49,53,96,97 show that these bonds tend to show PTE.

5 The second principle: the need for the network to support low-frequency vibrations

5.1 Why low frequencies matter

Merely identifying the existence of a group of atoms that may show the tension effect is insufficient to help us understand NTE. As we have said, the group of three atoms is a part of the larger network, and no atomic displacement can be independent of its local environment.

Consider a group of atoms with square coordination of the polyhedra in two dimensions, and with a connection between two squares that gives a tension effect, as shown in Fig. 19. The transverse displacement of the shared atom moving alone will bend the angles within the two polyhedra, which is likely to cost energy, Fig. 19a. This cost is removed if the two squares both rotate in concert, as shown in Fig. 19b. Such rotations can occur as RUMs if allowed by the network. The point is that quite probably the bending mode (Fig. 19a) will have rather higher frequency than the RUM (Fig. 19b).

Fig. 19 Two possible distortions of square coordinations with a transverse displacement of the central atom, one showing bending of the angles within the polyhedra (a), and the other showing rotations of the polyhedra (b).

Our qualitative understanding of the role of low frequencies is from the fact that the squared amplitude of atomic motion is proportional to the inverse-square of the frequency, 〈u²〉 ∝ ω⁻². That is, the lower the frequency, the larger the corresponding atomic motion. Much of the transverse atomic motion shown by the fluorine atoms in ScF₃, Fig. 12a, arises from low frequency modes.

Consider this triangle, which represents two bonds in the tension effect, each of length a₀/2, with a root-mean-square displacement of 〈u²〉^1/2:

The distance between the two horizontal ends, representing the polyhedra centres, is reduced by the tension effect from a₀ to a ≃ a₀ − 2〈u²〉/a₀ to the first order. The mean-square displacement associated with a vibration of angular frequency ω will be equal to 〈u²〉 = k_BT/mω², where T is the temperature and m is an effective mass. Thus we have a negative linear thermal expansivity equal to α = a₀⁻¹∂a/∂T = −2k_B/mω²a₀². The key point is that the NTE is large for smaller values of ω², and hence we can immediately see that the tension effect is more likely to arise from motions that cause lowest energy deformations of the structural polyhedra. This argument appears to be edging towards a “RUM theory of NTE”,¹⁴⁸ but at this point let us be prudently cautious, because in practice this need not be so.

This argument has actually sidestepped the third principle to be discussed later in Section 6, but we have at least established that requiring the network to permit low frequencies is a plausible factor in allowing it to show NTE. A firmer argument comes by considering a more-formal approach using Grüneisen theory, which we reviewed in Section 3.1. The important quantities are the individual mode Grüneisen parameters, defined in eqn (7). Considering the form of the equation in terms of ω², we see that γ_λ is the product of two quantities, ω_λ⁻² and ∂ω_λ²/∂V. Obviously, ω_λ⁻² is larger for smaller values of ω_λ. But using the argument we made earlier, that the eigenvalues of a matrix with lowest-values are the most sensitive values, it also follows that |∂ω_λ²/∂V| will be larger for lower frequencies. It follows from both factors that the mode Grüneisen parameters with the largest negative values will be those of the lower frequencies. This was clearly seen in ScF₃^86,88,90 and in most other systems for which there are comparable calculations.^{33,63,65,101,102,106,107,109,133,159}

We next ask the question, can network structures support the existence of low-frequency vibrations?

5.2 Rigidities of ceramic networks

The rigidity or flexibility of a network structure can be analysed in a first approximation using the counting of constraints and degrees of freedom, following the work of James Clerk Maxwell.¹⁶⁰ The central idea is that if the separate components of a structure have a total of F degrees of freedom, and if the connections between these components give a total number of constraints C, the structure will be rigid if C > F and flexible (or floppy) if F > C. This approach was initially introduced into condensed-matter sciences for the analysis of the rigidity and flexibility of chalcogenide glasses by Phillips and Thorpe.^161–166

In fact there are two ways to count constraints and degrees of freedom in network structures, as reviewed recently.⁴⁵ We prefer the approach of treating the polyhedra as rigid units, when appropriate, giving F = 6 degrees per polyhedron, and C = 3 constraints per shared vertex (from ensuring that the vertices of two connected polyhedra share the same coordinates x, y, z). The other approach is to treat atoms as the dynamic objects, with F = 3, to allow each bond to give one constraint, and to include angular constraints within the polyhedra where appropriate. Here we mostly follow the first approach, but incorporate bond constraints when we consider flexible polyhedra.

The important point is that if the total number of degrees of freedom exceeds the total number of constraints, F > C, the network is flexible in some way. But if the total number of constraints exceeds the total number of degrees of freedom, C > F, the network is over-constrained and hence rigid. In the marginal case, F = C, it is possible that some constraints are made redundant due to symmetry. This point has been discussed in some detail elsewhere in the context of the RUM model of phase transitions.^41,45

We start by considering the case where the polyhedra in the network are fully connected, so that the material has the general formula M_mT_nO_3m+2n and contains m MO₆ octahedra and n TO₄ tetrahedra. In the case where the polyhedra are rigid, we have F = 6m + 6n per formula unit. To count constraints, note that each vertex gives 3/2 constraints per polyhedron (3 total constraints, but shared between two polyhedra). It follows that C = 9m + 6n per formula unit, meaning that C–F = 3m. Thus the network is over-constrained unless m = 0, that is, only for tetrahedral networks might we expect to find RUMs.¹⁶⁷ In fact the cubic perovskite network is the only known exception, because the high-symmetry alignment of the octahedra allows for a significant number of constraints to become degenerate. This has been discussed elsewhere.^44,45,153

Let us now take the case where some of the atoms in the tetrahedra are non-bridging. Without making it too complicated, we can say that each tetrahedron has q bridging atoms. The chemical formula will be of the form M_mT_nO_{3m+n(4−q/2)}. Of course, we can't have half atoms, so clearly nq is constrained to be an even number. The number of degrees of freedom is the same as in the fully-connected case, but now we see a reduction in the number of constraints to C = 9m + 3nq/2, giving C–F = 3m − (6 − 3q/2)n per formula unit. This becomes flexible in the RUM-sense when C ≤ F. Let us take ZrW₂O₈ as our example: we have m = 1, n = 2 and q = 3, giving C = F = 18. In searching for new NTE materials, this simple criterion may give a possible guidance.

Now we consider the case where the tetrahedra are rigid, but the octahedra have some degree of bond-bending flexibility. We cite this example because this appears from DFT calculations to be a well-realised case.^65,168 Now we have to count slightly differently. We assign the same 6 degrees of freedom to the tetrahedra, but now we only count the degrees of freedom of the central atom in the octahedral sites, so F = 3m + 6n. Each bond in an octahedron acts as a constraint, giving C = 6m for the octahedra. Again, counting constraints for the tetrahedra will depend on how many non-bridging oxygen atoms are there. But we also need to be aware that if tetrahedra are connected to each other, we will need to add 3 constraints for each shared linkage.

Let us compare the cases of ZrW₂O₈ and ZrV₂O₇ as a practical example, in the two cases of rigid and flexible ZrO₆ octahedra, and with rigid WO₄ and VO₄ tetrahedra. For the case of rigid octahedra, we have F = 18 per formula unit in both cases. Because ZrV₂O₇ is fully connected, we have C = 21 per formula unit and hence C–F = 3. On the other hand, both WO₄ tetrahedra in ZrW₂O₈ have one non-bridging bond each, and as we discussed above we have the marginal case of F = C. The RUMs allowed in ZrW₂O₈, which arise as a result of symmetry causing some constraints to be degenerate, were found to have wave vectors lying on a surface of exotic shape in the three-dimensional reciprocal space.¹⁶⁹

When we take the case of flexible octahedra, for both cases we count F = 3 + 12 = 15 per formula unit, allowing for the rigid-body motions of two tetrahedra and the motion of the octahedral cation. For ZrW₂O₈, the only constraints we need to count are the six bond constraints within the ZrO₆ octahedra, giving F–C = 9 per formula unit. On the other hand, in ZrV₂O₇ we have to add 3 constraints for the linkage between pairs of VO₄ tetrahedra, giving F–C = 6 per formula unit. Again, ZrW₂O₈ will be more flexible than ZrV₂O₇.

These arguments could be modified if there are non-bridging bonds within the octahedra, and extended to the idea that no polyhedra are rigid, and that instead the only rigid entities are the bonds.^80,170 In such a case, for a fully connected network, and where we now consider each atom to have 3 degrees of freedom and the bonds provide the constraints, we have F = 12m + 9n and C = 6m + 4n, so F–C = 6m + 5n. However, this imputes too much flexibility to the network, and one consequence of such a situation is almost certainly that the network will be elastically soft. We do not consider that this situation is likely to be encountered, and we do not consider such an interpretation to be useful.

We should make two concluding comments in this section. The first is that if F = C, as in tetrahedral silica networks, this means that there is exactly zero flexibility unless, as mentioned above, symmetry makes some constraints redundant. In that case, we expect to see some RUMs with wave vectors at specific places in reciprocal space. In our studies of tetrahedral networks we often found that the RUMs lie on special planes in reciprocal space, which are often flat,^41,167 but to all our surprise, can also be on well-defined curved surfaces.¹⁷¹ We mentioned above the exotic surface found for ZrW₂O₈.¹⁶⁹ Also surprising was the finding that in low-density structures, like cubic zeolites of highest symmetry, we might find more than one RUM for each wave vector.^45,172,173 This gives a very high degree of flexibility. We stress that we do not consider this to be a failure of the Maxwell method¹⁶⁰ per se, but instead represents the fact that identifying constraints that are redundant is difficult. Furthermore, an inability to identify redundant constraints in the Maxwell counting will lead to us underestimating the degree of flexibility and never the other way round. The important cases where Maxwell counting overestimates the number of constraints appear to only be for the tetrahedral networks, where Maxwell gives the marginal case F = C, and in the case of corner-sharing octahedra on a cubic lattice (ScF₃ and perovskites). We have not identified other cases where there are significant numbers of redundant constraints. This statement applies also to the case of hybrid networks discussed in the next section.

The second comment is that the purpose of this analysis is to identify potential ways in which we might expect to find the necessary low frequencies in a network structure. But we cannot take this argument alone; it remains one of four principles, not a single principle. For example, in comparing the flexibility of ZrV₂O₇ and ZrW₂O₈ above, we found that the latter network is more flexible. Yet the NTE in ZrV₂O₇^18,67,68 is only slightly lower than that in ZrW₂O₈,^1,2,46 but by changing the chemical composition can significantly change the thermal expansion. The NTE in ZrV₂O₇ will turn into PTE by doping the vanadium sites with phosphorus.^67,68 In fact we will not consider the effects of chemical composition in detail here as a principle. but we will comment on some issues below (Section 9).

5.3 Rigidities of hybrid networks

In the case of ceramic networks, structural polyhedra often share oxygen or halide anions at their vertices (Cu₂O is of course a notable exception). However, it is possible for metal cations to be found in polyhedral coordination with specific atoms in linkage molecular ions. The simplest example is for the linkage molecular ion to be a simple dumbbell or chain, such as CN⁻ in Zn(CN)₂ and members of the Prussian blue family, or NCN²⁻ in Si(NCN)₂ (which has the same structure topology as Zn(CN)₂, except without the head-to-tail disorder).

For the structural topology of Zn(CN)₂, if we treat the tetrahedra and cyanide molecular anions as rigid objects, we have F = 16 per formula unit. The number of constraints, counting by vertices, is C = 12 per formula unit. There are actually two formula units in the primitive unit cell, so for the crystal as a whole we have F–C = 8 per unit cell. This means that we expect to find 8 RUMs for every wave vector. These are three modes in which the rotations of the Zn-centred tetrahedra cause transverse displacements of the cyanide molecular ions, three where the rotations give twists of the cyanide molecular ions, and also both of the transverse acoustic modes.

In principle the 8 RUMs in Zn(CN)₂ give the possibility for a tension effect, and hence for strong NTE. However, this assumes that there is a high flexibility at the vertices, as there is in silica, but it appears that there are strong directional forces. Thus the displacement mode has frequencies of around 4 THz, and the twist mode has frequencies of around 8 THz,¹⁰⁷ and their contribution to NTE appears to be relatively slight. Instead, it is the two acoustic modes that give the low frequencies that support the strong NTE in Zn(CN)₂.

In the case of the Prussian blue network with octahedra connected by shared cyanide molecular anions, the same counting gives F = 21 per M(CN)₃ formula unit. The number of constraints is C = 18. This gives F–C = 3, but there are two formula units in the primitive unit cell in the case of chemical composition MM′(CN)₆, giving a net F–C = 6. These correspond to two twist modes, two displacement modes, and two transverse acoustic modes.

The cyanide, or NCN²⁻, connectors give the simplest hybrid materials. In recent years we have seen an explosion in the number of metal–organic framework materials, a few of which crystallise with cubic symmetry and show NTE. Mostly, however, the linkages and linking groups of atoms are generally rather more complicated than we have discussed here. The isoreticular metal–organic frameworks have so-called ‘secondary building units’ instead of polyhedrally-connected atoms, such as O(ZnO₃)₄ clusters in MOF-5 (IRMOF-1), with the four ZrO₄ tetrahedra connected at the central oxygen atom, and with the outer oxygen atoms forming part of the linkage 1,4-benzenedicarboxylate molecular ions, O₂C–C₆H₄–CO₂. MOF-5 shows significant NTE.^130,131 The low-frequency flexibility in MOF-5 that leads to NTE appears to come from two sources. One is a low-frequency flexing of the 1,4-benzenedicarboxylate molecular ions,^131,133 and the other is the transverse acoustic mode.¹³³ However, there are relatively few isoreticular MOFs with cubic symmetry that show NTE.

One family of metal–organic framework materials with metal-centred tetrahedra are the zeolitic imidazolate framework structures. These exploit the fact that the angle subtended by the two bonds to metal cations from the nitrogen atoms of imidazolate molecular ions, C₃N₂H₃⁻ is around 145°, similar to the Si–O–Si angle in silicates. The ZIFs therefore are able to form crystal structures with networks that are also found within the family of zeolite structures. One of the better studied is ZIF-8, which has a cubic crystal structure with the same network topology as the sodalite zeolite structure. In this case the linkage molecular ion is methyl-imidazolate, C₃N₂CH₃H₂⁻ (mIm). The thermal expansion across the solid solution from Zn(mIm)₂ to Cd(mIm)₂ has been found to be positive.¹³⁹

5.4 Some case studies

5.4.1 Siliceous zeolites. Although NTE has been observed^8,123–126 and predicted^9,121,122 in zeolites, in truth the mechanisms have not been explored extensively. But we can make some useful remarks here.

Zeolites have been found to have surprisingly rich spectra of RUMs.^45,167,172 For example, the networks of the two cubic zeolites for which NTE have been observed, LTA^123,124 and FAU,^8,125,126 both having 4 RUMs for every wave vector,^167,172 together with a continuum of quasi-RUMs (QRUMs).⁴⁵ Such a rich spectrum of low-frequency modes will give the likely origin of a strong tension effect. However, the Si–O–Si angles in these two zeolites are typically around 145°, meaning that, geometrically, the tension effect will not be as strong as in the cases with near-linear bonds. But this is offset by the large number of phonon modes that will give NTE.

The one case where we have detailed phonon studies is for FAU,²² where many of the phonon branches below 4 THz are found to have a significant degree of RUM character and negative values of their mode Grüneisen parameters. The picture though is not as clear-cut as one might have imagined, because of the complexity of the phonon dispersion curves. With many branches at low frequencies, both RUMs and QRUMs, there will be a considerable degree of mixing of eigenvectors, as allowed by symmetry. We therefore cannot point to a single mode as being responsible for NTE.

5.4.2 Scandium fluoride. Scandium fluoride, ScF₃, is the clearest example of the relevant phonons for NTE being RUMs. Our own calculations of the phonon dispersion curves are shown in Fig. 20.⁸⁸ Aside from a few details (and taking account of issues such as the neglect of the correction for LO/TO splitting), these results are in reasonable agreement with those of other publications,^{85–87,90,174} and with X-ray inelastic scattering data.⁸² The first key point is the existence of the low-frequency band of modes extending for wave vectors between the M, , and R, , points in reciprocal space. The second key point is that in all calculations, these RUMs are found to have the largest negative values of the mode Grüneisen parameters. In fact the negative values are typically found to be very large. These modes are the well-established RUMs in this structure. Their atomic motions, namely alternating rotations of the octahedra within each layer (Fig. 19b), clearly gives the transverse motions involved in the tension effect. The facts that the motion requires no significant distortions of the polyhedra, and that there is a very low energy of flexing at the shared vertices, give the low frequency required for the large negative values of the mode Grüneisen parameters. The large negative values of the mode Grüneisen parameters are further enhanced by the relatively high stiffness of the Sc–F bonds, which shows normal positive thermal expansion^79,81 and vibrates at the highest frequency (the branch just below 20 THz in the dispersion curves shown in Fig. 20).

Fig. 20 Left shows the full set of calculated dispersion curves of ScF₃,⁸⁸ and right shows the same dispersion curves coloured according to the value of the mode Grüneisen parameters, saturating with values −4 (red) and +4 (green), with yellow representing a zero value. The same calculations were used for the generation of the data shown in Fig. 12.

Whilst we have identified a tension-effect mechanism for NTE in ScF₃, this is not sufficient to be able to say that we understand the mechanism of NTE in ScF₃. To do so, we need also to explain why NTE is not found in any oxide or halide cubic perovskite. Take for example SrTiO₃, where extensive data exist to demonstrate clearly that in its cubic phase this material shows normal positive thermal expansion.^175–177 Yet its lattice dynamics are remarkably similar to those of ScF₃.¹⁷⁸ Even more curious is that NTE is not found in AlF₃^91,179 or TiF₃.¹⁸⁰ We will identify the likely explanation when we discuss the third principle below, Section 6.

A lot of work has been carried out on NTE in the double-cation fluorides of general chemical formula MM′F₆, particularly where the one cation is divalent and the other is tetravalent. Experimental measurements have been performed on a number of examples,^76,181 as have lattice dynamics calculations^89,90,159 and molecular dynamics simulations.¹⁵⁸ The picture that emerges is similar to that discussed above regarding ScF₃, but some insights from this work will be discussed in Section 6.

5.4.3 Zirconium tungstate. Although initially it was thought that the origin of NTE in ZrW₂O₈ might be related to the RUMs in this material,¹⁶⁹ it is now believed that the mechanism is a little more complicated. It appears that the case of ZrW₂O₈ has rigid WO₄ tetrahedra but relatively flexible bending of the angles within the ZrO₆ octahedra, as discussed in some detail in Section 5.2. The phonon dispersion curves calculated by DFT methods⁶⁵ are shown in Fig. 21, where we compare two different representations with the standard representation of the dispersion curves. One highlights the modes that have negative and positive values of the mode Grüneisen parameters, and the other reflects the flexibility of the network.

Fig. 21 Left, calculated phonon dispersion curves of ZrW₂O₈ using DFT methods. Centre, the phonon branches have been coloured by the corresponding values of their mode Grüneisen parameters, with red representing negative values saturating at values of −6, blue representing positive values saturating for values above +6, and both colours fading to white as values tend towards zero. Right shows the level of matching of the phonon eigenvectors with those of a flexibility model with white to black representing matching from 0–100%.¹⁶⁸ This figure has been adapted from the published work of Rimmer et al.⁶⁵

One of the representations of the phonon dispersion curves in Fig. 21 has the different branches coloured according to the values of the mode Grüneisen parameters. It is seen that most of the branches for frequencies below 4 THz show negative values, and most of the branches for frequencies above 4 THz show positive values. The contrast with ScF₃ shown in Fig. 20 is striking. It is not just that ZrW₂O₈ has many more phonon branches (132 versus 12), but that the NTE is spread over many phonon branches with wave vectors all across the Brillouin zone.

The second comparison is with a model of the flexibility, in which the WO₄ tetrahedra are held to be rigid, as are the Zr–O bonds, but the O–Zr–O angles are completely flexible. The phonon modes from the DFT calculations have been projected onto the vibrations of the flexibility mode using an eigenvector matching algorithm.¹⁶⁸ Where the representation of the phonon branch is shaded dark, there is good matching, whereas a lighter shade of grey implies a low level of matching. It is seen that the branches with frequencies below 6 THz, more so below 4 THz, match well to this model, whereas all branches above 6 THz have very weak matching. This implies that the lower frequency modes, the ones that have negative values of their Grüneisen parameters, match this flexibility model well. It also follows that the phonon branches above 6 THz have motions that involve bending of the WO₄ tetrahedra. Not shown in the phonon dispersion diagrams of Fig. 21 are another set of modes with frequencies of 22–30 THz, which involve stretching of the Zr–O and W–O bonds. Matching to a flexibility model in which all bonds are rigid but all bond angles are flexible shows a complete match up to 12 THz,¹⁶⁸ but this fails to discriminate between rigid and flexible polyhedra.

The flexibility-matching study of the phonon dispersion curves included the case of rigid ZrO₆ and WO₄ polyhedra, representing the RUM model. This showed no signatures of RUMs for the wave vectors along the symmetry directions, except for the acoustic modes at a very small wave vector. This, however, is the only wave vector close to symmetry directions in which RUMs had been predicted.¹⁶⁹ The curved surface of RUMs does not intersect the symmetry directions in reciprocal space at any other point. We suggest, however, that for general points within the Brillouin zone, the RUM eigenvectors are able to mix with the eigenvectors of several other branches, because the phonons with general wave vectors have no symmetry other than the identity operation, and because there are many branches of relatively low frequency.

The conclusion from this study is that the low frequencies that support the negative values of mode Grüneisen parameters are enabled by some degree of flexibility of the ZrO₆ octahedra, whilst the WO₄ tetrahedra move as rigid objects over this frequency range.

5.4.4 Copper(I) oxide. We noted earlier (Section 2.1.8) that the structure of Cu₂O can be described both in terms of the positioning of the two types of atoms on simple cubic crystal lattices (interpenetrating body- and face-centred cubic lattices) and in terms of interpenetrating β-cristobalite-type frameworks of corner-linked OCu₄ tetrahedra. If we take the latter view, we might expect that the NTE can be explained in terms of a tension effect enabled by low-frequency RUMs. However, the point about lattice arrangement of atomic sites sends a red flag. In such an arrangement, the distances between Cu atoms in neighbouring tetrahedra are the same as the distances within the same tetrahedra. This immediately suggests two things, first that the electronic bonding may not lend itself to a network description of the vibrational properties, and second that there will be significant steric hindrance between neighbouring tetrahedra which will oppose large-amplitude RUM motions.

In fact it was found from lattice dynamics calculations¹⁰² that the mode eigenvectors of low-frequency modes that have some degree of RUM character also have a high degree of bond-bending flexibility. In fact the analysis of rigidity showed that there is a degree of rigidity of the linear O–Cu–O bond.¹⁰² In this case, we can treat the atoms as the objects with 3 degrees of freedom each, and the bonds as constraints, but adding 2 angular constraints for each O–Cu–O bond. Thus F = 9 and C = 8, giving one degree of flexibility per formula unit. Since there are two formula units per unit cell, we have two degrees of flexibility per wave vector. It appears that this flexibility is seen in the acoustic modes and that the NTE primarily arises from the acoustic modes.

This analysis has been quantified in the lattice dynamics results presented in Fig. 22.¹⁰² In this figure we compare the complete dispersion curves with the predictions of a flexibility model based on rigid linear O–Cu–O bonds but flexible angles, and with the calculation of mode Grüneisen parameters for each branch in the phonon dispersion curves. The role of acoustic modes giving rise to NTE is clear, as is the correlation with the description of the flexibility of the structure in which O–Cu–O moves as rigid entities but with flexible angles at the centres of the OCu₄ tetrahedra. It is important to note that assigning even more flexibility to this structure, for example in which the O–Cu–O angle is also flexible,¹⁷⁰ would predict too great a network flexibility and too soft an elastic stiffness.

Fig. 22 Left shows the calculated dispersion curves of Cu₂O. The same dispersion curves are shown in the centre panel coloured according to the size of the mode Grüneisen parameters, with red representing negative values and saturating at −8, blue representing positive values and saturating at +8, and both colours fading to white as values tend to zero. This suggests that the important modes for NTE are the transverse acoustic modes. The right panel shows a mapping of the phonon dispersion curves onto a flexibility model in which the linear O–Cu–O linkage is kept rigid, and the angular forces are zero. The motions associated with the acoustic modes reflect this rigidity model. This figure has been adapted from the work of Rimmer et al.¹⁰²

5.4.5 Zinc cyanide. Earlier, in Section 5.3, we pointed out that zinc cyanide will have 8 RUMs per wave vector. This suggests a very high degree of flexibility, but this analysis did not account for rigidity of the zinc–cyanide angle. In fact, some degree of rigidity of this angle leads to the twist and displacement modes having relatively high frequencies, particularly the twist mode. This is shown in calculations of the dispersion curves, presented in Fig. 23.¹⁰⁷ The band of modes with frequencies varying between 2 and 4 THz in the dispersion curve are the cyanide displacement RUMs, and these appear to have strongly-negative values of their mode Grüneisen parameters. On the other hand, the band of RUMs with the twist motions of the cyanide groups are at around 9 THz, and although they also have negative values of their mode Grüneisen parameters, their frequencies are too high for these to be significant. Instead, the important modes for NTE in zinc cyanide, as seen in Fig. 23, are the acoustic modes, a point that was initially noted from the lattice dynamics calculations by Joseph Zwanziger.¹⁰⁶ That is, the acoustic RUMs in zinc cyanide have the ability to give a significant tension effect, and because they are RUMs they have sufficiently-low frequency for the values of their mode Grüneisen parameters to be significant. In view of what we will say below, Section 6, it is pertinent to point out that these RUMs are spread across the whole of reciprocal space.

Fig. 23 Left shows the dispersion curves of Zn(CN)₂, calculated using a force field parameterised using ab initio calculations, performed using the virtual crystal approximation to simulate the head-to-tail orientational disorder of the cyanide anions. The centre panel shows the same dispersion curves coloured by the size of the mode Grüneisen parameters, where red represents negative values and blue represents positive values, with zero values coloured white. Right shows the mapping to a RUM flexibility model, where red represents rotational RUMs and blue represents acoustic (displacement) RUMs. On the far right is shown the corresponding atomic motions for a small wave vector. This figure has been adapted from the published work of Fang et al.¹⁰⁷

The case of single-network Cd(CN)₂ provides an interesting comparison with the double-network phase, not least because NTE is 50% stronger in the single-network case.¹⁰⁸ Zn(CN)₂ has a primitive lattice, space group Pn3̄m, with two formula units in the unit cell because of the double network. In the case of the single network, the crystal structure is face-centred cubic, space group Fd3̄m, also with two formula units per unit cell, but in this case they are part of the same network. Therefore the number of RUMs per wave vector is again 8. The enhanced NTE in the single network most likely arises from the single-network having lower frequencies of the acoustic modes, because there is no resistance from the second interpenetrating network, but this has apparently not been tested.

We noted earlier in our description of the crystal structure, Section 2.1.9, that there is a head-to-tail disorder of the orientations of the cyanide groups. Therefore most calculations have to be performed on a hypothetically ordered structure¹⁰⁶ or via parameterised force fields using the virtual crystal approximation.¹⁰⁷ In this context, two types of related compounds are interesting for study, because these have ordered orientations. The first is Si(NCN)₂, which we showed in Fig. 7 and discussed briefly in Section 2.1.9. The second has the general formula MB(CN)₄, where M is a monovalent metal cation. In this case, the boron atoms are bonded to the carbon atoms in the cyanide ion, and the M cations are bonded to the nitrogen atoms.

There are two known cubic forms of MB(CN)₄-type structures. They have the same interpenetrating β-cristobalite lattices, but with different alignments of the cations between the two substructures. In one case the M and B cations in the two sublattices are arranged in alternative (100) layers of each type, with space group P4̄3m. In this structure, one type of cation is located at the origin of the unit cell, and the other in the centre. One example is AgB(CN)₄,¹⁸² which shows an NTE that is not much smaller than that of Zn(CN)₂.¹⁰⁹ On the other hand, CuB(CN)₄ has the same structure but shows only small NTE.¹⁰⁹ Phonon calculations appear to suggest that, as in Zn(CN)₂, the acoustic modes play an important role in the NTE, with a contribution also from the lower-frequency optic RUM.¹⁰⁹ LiB(CN)₄ is another example that shows NTE.¹¹⁰ The other cubic form of this type of structure is demonstrated by NaB(CN)₄, which shows a large NTE.¹¹⁰ This has space group Fd3̄m, and in this case the neighbouring cations in the two sublattices along the [100] direction are of opposite type.¹¹⁰

In the MB(CN)₄ materials it is likely that the most rigid entity is the B(CN)₄ unit. Our own lattice dynamics calculations of the layer compound ZnB₂(CN)₄, which has a two-dimensional network of polyhedra connected by the cyanide ions, show that the B(CN)₄ units move as rigid objects for all the lower frequencies, with the modes showing significant bending of the C–B–C bond occurring for higher frequencies.¹⁶ On the other hand, the tetrahedra formed between the monovalent ions and the nitrogen atom are unlikely to be rigid. The monovalent-cation–nitrogen bond is likely to be formed by relatively weak electrostatic interactions, as in the alkali halides, with the angular rigidity formed from the repulsive interactions between neighbouring nitrogen atoms rather than from any covalent bonding. This may account for the point made by Gao et al.¹¹⁰ about some of the phonon modes showing a lack of rigidity. However, the point is not quantified. Our own unpublished calculations of LiB(CN)₄¹⁸³ indicate that in the lowest-frequency optic mode at a zero wave vector and a frequency of 65.6 cm⁻¹ (1.97 THz), which is part of a branch of optic modes with consistently negative values of the mode Grüneisen parameters, has mode eigenvectors characteristic of a RUM, with no distortions of the B(CN)₄ and LiN₄ tetrahedra. Instead these two polyhedra rotate with transverse displacements of the cyanide anions, giving rise to an obvious tension effect.

If the BC₄ and MN₄ polyhedra move as rigid units, the flexibility analysis is identical to that for Zn(CN)₂ and Si(NCN)₂, showing 2 transverse acoustic RUMs and 6 optic RUMs for every wave vector. If the B(CN)₄ unit is the rigid unit, this reduces the problem to the topology of β-cristobalite, albeit with different Bravais lattices and with the flexibility reduced to planes of RUMs in the [110] zones of wave vectors (wave vectors of the form [ζ, ±ζ, ξ]).¹⁸⁴ These include a triple-degenerate set of optic modes at a zero wave vector, and 2, 1 and 3 optic RUMs along the [100], [110] and [111] directions respectively. This degree of rigidity explains why one of the lowest-frequency optic modes double in frequency as the wave vector goes from zero (the Γ point) to (X point) in AgB(CN)₄.¹⁰⁹ On the other hand, if the MN₄ tetrahedra are flexible in terms of bending the N–M–N angle, we would get 5 RUMs per wave vector. It is quite likely that this all-or-nothing approach to rigidity is too extreme, but certainly to quantify the rigidity in the various examples of MB(CN)₄ materials would be interesting and important. It is interesting to note that we find a cyanide twisting RUM in LiB(CN)₄ at the high frequency of 335 cm⁻¹ (10 THz),¹⁸³ similar to the same RUM in Zn(CN)₄.¹⁰⁷

The same discussion will apply to other cyanides, such as the members of the Prussian blue family, Section 2.1.10. As we pointed out earlier, Section 5.3, the Prussian blue network is able to support 6 RUMs per wave vector, of which 2 are transverse acoustic modes. Although not yet clearly demonstrated, we expect that the transverse acoustic modes are most likely to be the important vibrations in giving NTE in these materials.

6 The third principle: the distribution of relevant normal modes must span a sufficient volume of reciprocal space

When we discussed ScF₃ above (Section 5.4.2) we made two points. The first concerned the low-frequency RUMs in this material, which been clearly identified through lattice dynamics calculations, Fig. 20,^{85–88,90,174} and inelastic X-ray scattering measurements.⁸² The lattice dynamics calculations have shown that the RUMs have extremely large negative values of the mode Grüneisen parameters, more so than for other vibrations in the crystal. The demonstration of the importance of these modes is given in the analysis presented in Fig. 12.⁸⁸ It is clear that the mechanism for NTE in ScF₃ involves the RUMs at its heart.

The second point is that many crystals share both the structure and features of the lattice dynamics of ScF₃ yet without showing NTE.^{91,175–177,179,180} The cubic phase of SrTiO₃ has phonon dispersion curves¹⁷⁸ that share much in common with those of ScF₃,⁸⁸ without having NTE. Therefore it is not sufficient to have identified the RUMs as giving the mechanism for NTE in ScF₃ – we also need to explain why this mechanism is important in some cases but not in others.

We noted earlier (Section 3.6) that the RUMs in ScF₃ and in all cubic perovskites lie along the edges of the Brillouin zone.† The important point to understand is that the volume of the line occupied by the wave vectors along the edge of the Brillouin zone is a vanishingly small fraction of the three-dimensional volume of the whole of the Brillouin zone. Thus, although the set of RUMs may contribute extraordinarily large negative Grüneisen parameters to the overall sum, in truth they are large irrelevant by themselves because the number of RUMs is almost infinitesimally small compared to the total number of vibrations.

In practice, however, the dispersion curves form continuous distributions of frequencies in reciprocal space, so the phonons whose wave vectors are very close to the edges of the Brillouin zone will also have low frequencies close to those of the RUMs, with correspondingly-large negative values of mode Grüneisen parameters. As the wave vector moves further from the edge of the Brillouin zone, the frequencies will increase, the modes will have less of the rotational RUM character, and increasingly the mode eigenvectors will include atomic motions that bend the bonds in the ScF₆ octahedra. Thus the values of the mode Grüneisen parameters will become less negative, possibly becoming positive. This is seen in the phonon dispersion curves in Fig. 20. The RUMs are the low-frequency curves connecting the M, , and R, , points. The R point represents the meeting of three such lines, so the RUM at R is triply degenerate. Going from R → M we see the two modes (doubly degenerate) that cease to be RUMs increase in frequency, with a commensurate lessening of the negative value of the mode Grüneisen parameter. The same is seen in the directions R → Γ, [0,0,0], and M → X, .

It appears, therefore, that there is a cylindrical volume around the edges of the Brillouin zone, within which there is a sufficient number of wave vectors with RUM-like modes – modes with the necessary low frequencies and tension-effect motions to give negative mode Grüneisen parameters – to enable the crystal to have an overall NTE. The corollary of this is that NTE is not seen in similar systems because the volume of this cylinder around the edges of the Brillouin zone may not be large enough. So what controls the size of this volume?

The size of the relevant volume around the edges of the Brillouin zone is set by how quickly the phonons that are no longer RUMs increase in frequency. In turn, this is set by the forces associated with bending the bonds within the octahedra. That is, the more rigid are the octahedra, the smaller the relevant volume around the edges of the Brillouin zone, and the less likely is the chance for the RUMs to have sufficient weight to give NTE. We tested this idea using the simple model shown in Fig. 24a.⁸¹ The model has only three forces. The first force is to control the stiffness of the Sc–F bond. A Morse potential was used for this,

E(r) = D[exp(−2α(r − r₀)) − 2 exp(−α(r − r₀))], (23)

with parameters which can be tuned to give the correct Sc–F distance and, approximately, the value of Sc–F stretching frequency in the phonon dispersion curves calculated by DFT. The second force is a bond-bending potential for the F–Sc–F right angle in the ScF₆ octahedra. We used a function of the formwhere θ is the angle and θ₀ = 90°. The third force is a bond-bending potential for the Sc–F–Sc linear connection. We used a function of the formwhere φ is the angle, with an equilibrium value of 180°. These two angles are identified in Fig. 24a.

Fig. 24 The simple model for ScF₃ used to assess the important of the size of the bond angle flexing forces.⁸¹ (a) shows the definition of the model, for which a Morse potential is used to describe the Sc–F bond length, with angular forces for the right-angle F–Sc–F bond angle θ, eqn (24), and for the linear Sc–F–Sc angle φ, eqn (25). The model parameters were tuned to best reproduce the DFT lattice dynamics calculations shown in Fig. 20. (b) Shows the effect of varying the force constant A associated with the Sc–F–Sc angle φ, with the smallest value representing the starting model. (c) Shows the effect of varying the force constant k associated with the F–Sc–F angle θ, again with the smallest value representing the starting model. This figure has been adapted from the work of Dove et al.⁸¹

The parameter A was tuned so that the model gave the value of the frequency of the RUM along the line of wave vectors from R → M, consistent with the DFT calculations of the lattice dynamics.⁸⁸ The parameter k directly controls the dispersion of the transverse acoustic mode along the Γ → X direction, so it was tuned to give the appropriate value of the elastic constant C₄₄. Whilst the model lacks electrostatic interactions, and whilst it does not contain interactions required to give a non-zero value of the elastic constant C₁₂, it nevertheless gave a full set of phonon dispersion curves in remarkable agreement with the DFT dispersion curves.⁸¹

The model was then explored using the molecular dynamics simulation method with varying values of the two force constants A and k. The results for varying the value of A are shown in Fig. 24b. As this bond becomes stiffer, the frequency of the RUM will increase, and so we see the lowering of the NTE, until at a value of A that is 40 times larger than the tuned value the thermal expansion turns positive. Qualitatively this is not surprising, given that the effect of increasing the value of A is to stiffen the structure, and through the raising of the RUM frequencies to decrease the negative values of the mode Grüneisen parameters along the RUM lines – remember the additional factor of ω⁻² in eqn (7).

The second set of results, for varying the value of k, are shown in Fig. 24c, and are particularly interesting. As the ScF₆ octahedra become more rigid, the NTE decreases, which corresponds to the narrowing of the cylinder of important wave vectors around the edges of the Brillouin zone. What is interesting, certainly in comparison with the variation of the value of A, is that the thermal expansion reaches zero with only a three-fold increase in the value of k.

It is possible to use a variation of this model in which the F–Sc–F right-angle potential, eqn (24), is replaced by a Morse potential for the F–F bond with an equilibrium distance that is times that for the Sc–F potential. This has the effect of giving a non-zero C₁₂ = C₄₄, which is actually as given by the DFT calculations,⁸⁸ because compressing a pair of opposing S–F bonds will cause compression of the F–F bonds, which is relaxed by a balanced expansion of the perpendicular Sc–F bonds. The cost is that the value of C₁₁ is also affected by the F–F force constant, thereby perhaps breaking some of the orthogonality of the model, but this might be considered a small price compared to the feeling that C₁₂ = 0 is too unrealistic. We developed such a model tuned by comparing the elastic constants given by this model with those calculated by DFT methods,⁸⁸ showing a 1% agreement, and with a calculated phonon density of states similar to that calculated by DFT. The thermal expansion of this model is also similar to that given in the DFT calculations. As with the model discussed before, with a direct F–Sc–F angle function,⁸⁸ we found from the lattice dynamics calculations that the NTE lessened and became PTE when the F–F force constant was increased by the same factor as the increase to the force constant for the F–Sc–F angular potential energy seen in Fig. 24c. Thus the conclusions regarding the importance of the flexibility of the ScF₆ octahedra are robust regarding details of the model. We will discuss the use of such models in Section 10.

We can now compare the stiffness of the octahedra in ScF₃ with that for some other systems, namely the TiO₆ octahedra in SrTiO₃ and the PbI₆ octahedra in CsPbI₃ and methylammonium lead iodide. We have performed studies on these systems^145,188,189 using neutron total scattering¹⁹⁰ analysed by the Reverse Monte Carlo method.¹⁹¹ From the realistic atomic configurations these studies have given, we have been able to extract the contributions to the atomic motions from stretching and bending of the bonds within the octahedra, and from whole-body rotations of the octahedra, using the GASP method.^43,185–187 The results for ScF₃ and the model described above,⁸¹ methylammonium lead iodide,¹⁸⁸ CsPbI₃¹⁸⁹ and SrTiO₃,¹⁴⁵ are compared in Fig. 25. It can be seen that ScF₃ has the most flexible octahedra; this correlates with the fact that only ScF₃ shows NTE, and is consistent with our contention here that the bond-bending flexibility increases the volume of reciprocal space that is required for the RUM-type motions to be significant.

Fig. 25 GASP results^43,185–187 for some different cubic networks of linked octahedra, ScF₃ and its model,⁸¹ methylammonium lead iodide (MAPbI₃),¹⁸⁸ CsPbI₃¹⁸⁹ and SrTiO₃.¹⁴⁵

It is important to note that actually ScF₃ is NTE only with a fine degree of sensitivity to the balance between the phonons that give NTE and those that give PTE. This is clearly seen in Fig. 12, particularly Fig. 12c, which shows a very large cancellation between modes with negative and positive values of the mode Grüneisen parameters. One can see that with only a slight reduction in the weighting of the modes with negative values of the mode Grüneisen parameters the thermal expansion will become positive.

This discussion might explain why the related fluorides AlF₃^91,179 and TiF₃¹⁸⁰ show positive thermal expansion. Although the necessary calculations have not yet been performed, we can learn something from the recent calculations on a range of MF₃ structures by Koiso et al.¹⁹² These calculations were primarily designed to understand why some materials such as AlF₃ show displacive phase transitions whereas ScF₃ does not. It was suggested that the RUM-type distortions in AlF₃, GaF₃ and InF₃ arise from a second-order Jahn–Teller interaction that gives distortion of the M–F–M linkage and which is strong enough to overcome the Coulomb repulsion arising from the distortion. As part of this study the authors calculated the phonon dispersion curves of the cubic phase. It is clear from the results for the acoustic modes of AlF₃ that the AlF₆ octahedra are much stiffer than the ScF₆ octahedra, consistent with the previous discussion. On the other hand, YF₃ and La₃ were predicted to show NTE, and their dispersion curves show that the acoustic modes have low frequencies, indicative of more flexible octahedra.

Confirmation of the same point also comes from the MM′F₆ systems. It follows that a cation giving softer octahedra, such as the calcium cation, might lead to an increase in the NTE. This point has been confirmed experimentally⁷⁶ and by lattice dynamics calculations.⁹⁰

It has been suggested that we could take this argument to an extreme limit in which it is assumed that there is no bond-bending interaction within the ScF₆ octahedra, but that instead the orientations of the Sc–F bonds can librate as independent Einstein oscillators.⁸⁰ ‡ This would mean that the relevant phonons for NTE would be distributed across the whole of the Brillouin zone, but detailed lattice dynamics calculations demonstrate that the relevant modes for NTE do actually lie within a cylinder around the edges of the Brillouin zone.⁸⁸ Moreover, with such a flexibility we would expect the value of the elastic constant C₄₄ to be extremely low – it isn't!⁸⁸ Not only is such a model unsatisfactory because it contradicts the elasticity seen in the experimental data for the acoustic phonons,⁸² it is also not a plausible route towards understanding NTE because it requires that crystals are unrealistically soft.

The key point is that NTE does not arise from a single phonon or a specific set of phonons, but from a sufficiently-large number of phonon modes that their overall contribution to the thermal expansion outweighs the contribution from the many phonons that would normally give PTE. When the important phonons lie on specific wave vectors, or on lines or planes in reciprocal space, it is necessary that their characteristics (low frequency, tension effect motions) bleed outwards into a sufficiently large volume of reciprocal space – our third principle.

For some systems the important phonon modes span across much of reciprocal space. This includes the zeolites, the hybrid cyanide structures, and oxides such as ZrW₂O₈ (Fig. 21) with flexible octahedra. In this case, the relevant phonon modes do not need to have large negative values of mode Grüneisen parameters, and often they don't – the large sizes seen in ScF₃^86,88 are not common. In this case, several phonon branches with higher frequencies and lower negative values of mode Grüneisen parameters are weighted by the fact that they occur for wave vectors across most of the Brillouin zone.

7 The fourth principle: disorder can be a limiting factor

Earlier we posed the question of why we do not find ever larger values of negative expansivity. If we have more flexible networks, which clearly fulfil the principles of network connectivity, low frequencies, and a distribution of NTE modes across wide volumes of reciprocal space, should we perhaps expect to find larger negative values of expansivity? But in general this is not the case.

The case of β-cristobalite presents an interesting challenge. It is surely a flexible structure. Its idealised crystal structure is shown in Fig. 1. This structure has linear Si–O–Si linkages, with large transverse displacements of the oxygen atom,¹⁹⁵ characteristic of motions associated with the tension effect (compare with ScF₃, Fig. 12a). The large-amplitude motions arise from the RUM flexibility,^41,42,196 which in turn leads to the existence of a displacive phase transition.^195,196 The presence of a significant number of RUMs in β-cristobalite has been demonstrated by inelastic neutron scattering measurements,¹⁹⁶ quantified by GASP analysis of configurations generated from neutron total scattering measurements analysed by RMC,⁴³ confirmed by the diffuse scattering calculated from the RMC configurations,⁴⁰ and measured by electron diffraction.¹⁹⁷

In spite of β-cristobalite appearing to have all the ingredients required for NTE, β-cristobalite shows PTE up to around 1400 K³⁹ and NTE only at higher temperatures,³⁸ Fig. 26. This is remarkable given that Cu₂O also crystallises with the β-cristobalite structure, but in this case there are two interpenetrating β-cristobalite lattices with OCu₄ tetrahedra, which might be thought to impede the fluctuations required to give NTE. Yet Cu₂O shows NTE at lower temperatures^95–97 when β-cristobalite does not.

Fig. 26 The variation of the lattice parameter of cubic β-cristobalite with temperature, merging data from two sources as open³⁹ and filled³⁸ black circles. The continuous black line is a guide to the eye fitted through the data. The red curve is a calculation from the simulations discussed in Section 7.¹⁹⁴

Another example is Si(NCN)₂.¹¹² This has the same crystal structure as Zn(CN)₂ (Fig. 7) with ordered NCN²⁻ molecular anions replacing the disordered CN⁻ anions. DFT calculations suggest that Si(NCN)₂ has a high degree of flexibility,¹¹¹ yet the NTE seen in Si(NCN)₂¹¹² is much smaller than that in Zn(CN)₂.¹⁰³

It is probable that the high degree of flexibility of both structures leads to the presence of a significant degree of structural disorder. In fact this has been quantified much more for β-cristobalite than that for Si(NCN)₂ in terms of structural^40,195 and simulation^198–200 studies. The total scattering data analysed using the RMC method showed that the O atoms move continuously around an annulus surrounding the Si–Si linkage, with a distribution of positions that look more like a doughnut.⁴⁰ Prior to this work a model was proposed in which the oxygen atoms occupy six sites on this annulus, but in truth the fitting of this model to the neutron powder diffraction data was little better than that for the idealised model with large anisotropic atomic displacement parameters.¹⁹⁵ The RMC results show a much more continuous, and broad, distribution of oxygen positions around the annulus. Put another way, they show a broad distribution of the orientations of the Si–O bonds and hence of the SiO₄ tetrahedra.

We have recently posed the question of whether the presence of a significant degree of dynamic structure disorder is responsible for the PTE in β-cristobalite and the relatively low NTE in Si(NCN)₂.¹⁹⁴ We have constructed some simple models, analogous to the model for ScF₃, Fig. 24, both for β-cristobalite and Si(NCN)₂, the latter using rigid NCN molecules. The model allows for stretching of the Si–O and Si–N bonds, bending of the O–Si–O and N–Si–N tetrahedral angles (an ideal value of 109.47°), and flexing of the Si–O–Si and Si–NCN angles. The assumption is made that the bond and tetrahedral angle stiffnesses are large, but the second angle is more flexible. The thermal expansion shown by the model was calculated by performing many molecular dynamics simulations across a wide range of temperatures for a range of stiffnesses of the second angle.

We consider here mostly the results for one version of the model for β-cristobalite, in which we have a potential of the form of eqn (24) for the Si–O–Si angle, which we denote as θ, with θ₀ = 145°. We have a similar function for the tetrahedral O–Si–O angle, whose force constant is held fixed at a large constant value. We also create a stiff Si–O bond using a Morse function. In this case the Morse function allows for a small positive thermal expansion of the bond.

Results for the variation of the lattice parameter with temperature for different values of the force constant k, from small to large, are shown in Fig. 27. The results show different behaviours for the two extreme values of k, large and very small values, within the results for intermediate values of k showing an interesting interpolation between the two extreme limits. This we need to unpick.

Fig. 27 Results for the simple model showing thermal expansion in β-cristobalite, in each case plotting lattice parameters obtained from molecular dynamics simulations as functions of temperature for different values of the Si–O–Si bond bending force constant (units of eV), as given in the legend. The corresponding value of the force constant for the tetrahedral angle is 3 eV.

In the case of very small values of k, the network is very floppy, and RUMs will be easily excited. However, it is important to know that the RUMs in β-cristobalite lie on planes in reciprocal space that are normal to the lattice vectors of the form 〈110〉, that is, with wave vectors of the general form {ζ,±ζ,ξ}. Thus the β-cristobalite network has a more constrained RUM flexibility than do the zeolites, and large-amplitude RUM motions will require a small distortion of the tetrahedral Si–O–Si angles. This means that for a small value of k, at low temperatures the crystal structure will align close to its state of maximum volume because that will involve minimal bending of the Si–O–Si angles. What is seen is that there is a rapid drop in the value of the lattice parameter, clearly with a large NTE, until saturating with nearly zero thermal expansion. The saturation is at the value of the lattice parameter at which the orientational disorder of the SiO₄ can no longer increase within the constraint of the network.

This result is what was also seen in simulations of the model for θ₀ = 0 and also in the simulations of the models for Si(NCN)₂ with and without a corresponding angular force. In both cases, the size of the shrinkage decreases as the angle force constant increases. This implies that the relatively low thermal expansion seen in Si(NCN)₂ may result from the saturation of the disorder. In the case of Si(NCN)₂ there is a phase transition to an unknown crystal structure, so the high NTE this model would predict at low temperatures cannot be observed.

In the other limit of large values of k, the bending term constrains θ to be close to 145°. This causes a small distortion of the SiO₄ tetrahedra above the thermal motion, expanding the standard deviation of the O–Si–O tetrahedral angle by about 0.6° at the lowest temperature, above the thermal value of 2°. Upon heating 〈θ〉 remains close to 145°, with a variance that increases with temperature. In this case there is no scope for the operation of a tension effect, because the network is locked into a constrained disordered state, and as a result the system shows only PTE. In fact the linear thermal expansion of the network is about 2/3rd of the thermal expansion of the Si–O bond, showing that fluctuations do reduce the thermal expansion to a small extent, but not by enough to give NTE.

The intermediate cases show a number of very interesting properties. First is that with a not-insignificant value of the Si–O–Si bending force constant k there is a balance between the energy costs of straightening the Si–O–Si angle and bending the O–Si–O tetrahedral angles. In practice this pulls the weight of the distribution of Si–O–Si angles towards 145°, but upon heating the increasing thermal fluctuations actually weight the distribution away from 145° and more towards 180°. This gives rise to the initial positive thermal expansion. However, upon further heating the distribution broadens, which can then only mean spreading away from 180°. However, the pull is not specifically towards 145° but towards larger angles, which then gives rise to NTE, but again saturating exactly as for the case of the weakest value of the Si–O–Si force constant k.

The interesting point from this is that in this model we see an initial PTE at lower temperatures, followed by NTE at higher temperatures. This is exactly as seen in β-cristobalite. We can in fact adjust the parameters of the model to give a reasonable agreement with the experimental data for β-cristobalite, as shown by the red curve in Fig. 26. We can also adjust the parameters to give a good agreement with the NTE seen in Si(NCN)₂.¹⁹⁴

One very curious feature in Fig. 27 is that the lattice parameters for all values of the Si–O–Si force constant k are identical at a temperature of 598 K. If we increase the value of the force constant associated with the O–Si–O tetrahedral angles there remains a point of concurrency but it moves to higher temperatures. In fact there is a linear – but not proportional – relationship between the concurrency temperature and the value of the O–Si–O force constant.

What this model has shown for both β-cristobalite – which we have discussed in more detail – and Si(NCN)₂ is that when we increase the features that can give NTE, as discussed in the first three principles, we have larger thermal fluctuations which give rise to local dynamic disorder and which then limits the size of the NTE. This leads to our fourth principle that the thermal disorder created by the network flexibility that gives NTE ends up limiting the size of the NTE that can be achieved. We might remark that even for materials that show larger NTE, it may well be the case that disorder is still a limiting factor, but this is a point that has not yet been studied.

Experimentally the saturation of the NTE with increased temperature due to disorder will be hard to differentiate from the effects of phonon–phonon anharmonicity increasing the frequencies of the modes associated with the tension effect and from the effects of exciting higher-frequency modes that contribute to PTE with positive mode Grüneisen parameters. The latter effects can be indicated by lattice dynamics calculations.

The effect due to bending of the linkage of three atoms, as in β-cristobalite, may be more common, but again difficult to identify. We can identify cases where we do not expect to see such an effect. According to the recent calculations by Koiso et al.¹⁹² it is probably the case that the Sc–F–Sc linkage in ScF₃ favours the linear arrangement. Previous work on Zn(CN)₂ also suggests that the Zn–(CN) bond also prefers to be linear, which is quite likely to be true for several other cyanide compounds. Thus the NTE in both ScF₃ and Zn(CN)₂ will not be affected by this type of disorder. On the other hand, given that it appears that AlF₃ prefers a bent Al–F–Al linkage,¹⁹² it is possible that in addition to the effects from the stiffness of the AlF₆ octahedra limiting the impact of the line of RUMs (recall the discussion regarding ScF₃ in Section 6), it is possible that the disorder of the positions of the fluorine atoms may also mitigate against NTE in the cubic phase of this material.

It might also be the case that ZrV₂O₇ shows the same behaviour as that seen in β-cristobalite. The shared oxygen in the pyrovanadate, V₂O₇, cluster is linear on average but prefers to be bent. Data published by John Evans¹⁸ show that ZrV₂O₇ shows PTE at temperatures immediately above the phase transition at around 375 K, reaching a maximum value of the lattice parameter at a temperature around 10 K higher and then showing NTE at higher temperatures. This resembles the behaviour in β-cristobalite, albeit with different ranges of temperature, and disorder of the orientation of the V–O bonds may be the explanation. The same appears to be the cases or uranium and thorium pyrovanadates.⁶⁷ The case of PTE in zirconium pyrophosphate, and also titanium pyrophosphate,⁶⁷ may also be due to orientational disorder of the P–O bonds over all temperatures measured.

8 The challenge of anisotropic materials

The ideas presented in this article will apply equally to anisotropic materials, but in this section we wish to urge caution when thinking about anisotropic thermal expansion. Whilst some anisotropic materials, such as Y₂W₃O₁₂,¹² show NTE along all three directions, many show NTE only in one (uniaxial) or two (area) directions, with PTE along the other axes.

It is commonly assumed that the tension effect applies to NTE in anisotropic materials, but this may not always be the case. In particular, there are two other causes of NTE in anisotropic materials, which frequently may be more important. One cause of uniaxial or area NTE is from elastic anisotropy. For anisotropic systems, the Grüneisen equation for thermal expansion in isotropic systems, eqn (12), does not apply. This is a point that is not always appreciated. The correct equation is

where j represents the six possible components of the strain tensor (tensile and shear) using the Voigt notation, and S_ij represents any of the components of the fourth-rank elastic compliance tensor. The overall Grüneisen parameters [small gamma, Greek, macron]_j are then defined with respect to the strains rather than volume, formed aswhere is a component of the second-rank strain tensor, and includes both tensile and shear strains.

It is tempting to assign the same significance to volume and volumetric thermal expansivity α_V to an anisotropic crystal showing NTE as to NTE in isotropic materials, but this can be misleading. The volume strain is actually a sum of the individual tensile strains,

This means that cannot be considered in isolation from separate analysis of the individual tensile strain components. For certain, this prevents the formulation of a meaningful calculation of mode Grüneisen parameters.

The point of discussing eqn (26) in this context is that it is perfectly possible for NTE to arise purely from the signs of the components of the elastic compliance tensor with positive values of the overall strain Grüneisen parameters. If the off-diagonal components of the elastic constant tensor C have positive values, it follows that the off-diagonal components of the elastic compliance tensor S = C⁻¹ are necessarily negative, and these can be large enough to give axial or area NTE even when all the overall strain Grüneisen parameters are positive. This point was discussed back in 1972⁴ and has been demonstrated in a few recent studies.^15,16,201 One of these cases, the silicate cordierite, is interesting in showing axial NTE, when at some temperatures all three overall strain Grüneisen parameters have negative values and at other temperatures all have positive values.

A second origin of NTE in anisotropic materials is from a structural phase transition, where the symmetry-breaking distortion of the crystal structure is accompanied by spontaneous strains. If these strains are positive, and if the background positive thermal expansion is not too large, the positive spontaneous strains will give the impression of NTE upon heating up to the phase transition. Examples are PbTiO₃,^202–205 Cu₂P₂O₇,^17,206,207 ammonium sulphate,^208,209 and Cd₂Re₂O₇.^210,211 Interestingly, in the last three of these cases the phase transition appears to be near-tricritical, which means that the NTE actually diverges at the phase transition. This is astonishing and gives the chance to control very large NTE, which compares with the point we have made that isotropic crystals may have a natural limit to the size of their NTE due to dynamic disorder, Section 7.

In these cases, it appears that there is a temptation to immediately perform DFT calculations of Grüneisen parameters. As we have pointed out elsewhere,¹⁷ it is highly likely that it will be possible to calculate negative values of the mode Grüneisen parameters, even when using the isotropic formulae. This is a direct consequence of the coupling between strains and order parameter. But the calculation of a negative value of a mode Grüneisen parameter when NTE is associated with a phase transition should never be taken as implying the existence of a tension effect, particularly when the thermal expansion in the high-temperature phase is positive.^207,211

The important point is that whilst there may be a tension effect operating in some anisotropic materials, sometimes this may be only for certain strain components. One example is for CuCN,¹⁵⁵ where there is a clear tension effect operating along the Cu–CN–Cu chains.²¹² Graphene is the clearest two-dimensional example, but here the tension effect operates only within a one-atom-thick layer. More interesting, perhaps, is recent work on the network layer material ZnB₂(CN)₈.^16,213 Although the biaxial NTE in this material is caused by the anisotropic elasticity, the origin of a very low overall strain Grüneisen parameter for the relevant area strain has been understood in terms of the same principles outlined in this article.¹⁶ However, in other cases, a tension effect may be slight or non-existent, with the major effects arising from the elastic anisotropy or the phase transition. We therefore caution against transferring our clear understanding of NTE in isotropic materials directly to anisotropic materials.

9 A question of balance

The principles outlined in some detail here will explain why we might find NTE in a particular material and can help also interpret the absence of NTE in other cases. But we have also argued that the features that give the phonons with negative values of their mode Grüneisen parameters are in a competitive balance with the short-range bonding interactions that favour PTE. This balance is often seen in the fact that for many materials the NTE turns to PTE at higher temperatures, as the population of modes with positive mode Grüneisen parameters but high frequencies are excited. This effect is seen particularly in Cu₂O,^95–97 but also in ScF₃,⁷⁴ for we have argued that there is a fine balance between vibrations with negative and positive values of mode Grüneisen parameters across the full spectrum of frequencies.

Particularly pertinent is the question of why, in crystal structures that are largely identical but have different chemical composition, this balance can be pushed dramatically from one side to another? For example, why do we see NTE in ScF₃ but not in AlF₃,¹⁷⁹ why in ZrV₂O₇ but not in ZrP₂O₇,⁶⁷ why strongly in AgB(CN)₄ but only weakly in CuB(CN)₄?¹⁰⁹

We have addressed some issues regarding what controls the balance. These include that of the distribution of the wave vectors of the phonons with large negative mode Grüneisen parameters and the relative volume they occupy within the Brillouin zone (Section 3.2). We have also seen from the simple model introduced in Section 3.5 that the size of the negative mode Grüneisen parameters depends on two factors. The first one, in the numerator, is the stiffness of the bond, which in turn determines the coupling of the transverse motion of the central atom in the tension effect to strain. The stiffer the bond, the higher the coupling, leading to promotion of NTE. The second, in the denominator, is the square of the mode frequency, which clearly leads to promotion of NTE if the frequency is low. All of these will be affected by the chemical composition. These factors may all be characterised in some way by the atomic size.¹⁰⁹ Larger atoms may result in softer bending potentials that lead to an increase in the relevant volume of reciprocal space and lead to lower frequencies, but at the same time larger atoms may have softer nearest-neighbour interactions, which reduce the operation of the tension effect.

It is not our intention to develop this point further. Our goal in this article has been to review the mechanisms that are required to give NTE, but in practice even with all the mechanisms in place, the issue of how the balance tips between NTE and PTE will depend to a greater or lesser extent on the detailed chemistry. These factors will affect the mechanisms in different and often competing wells, affecting both the numerator and denominator in the equation for thermal expansion of isotropic systems, eqn (12). There is scope for more quantitative research on this issue.

As a side issue, we should address the question of whether in the cases where there is a fine balance – such as in ScF₃, as evidenced in Fig. 12 – calculations by slightly different methods may give different results. For example, DFT calculations performed by different software packages, or even by different workers using the same software package, might have different flavours of the exchange–correlation functionals, different pseudopotentials, or even different cutoffs. This is a serious point because often different calculations of the phonon dispersion curves for the same system give different results, as for example for ScF₃^{85–87,90,190,214,215} and the related CaZrF₆.^76,89,90,159 It is often found that the extreme values of mode Grüneisen parameters may vary between different calculations. However, the broad results for the calculation of thermal expansion usually show reasonable qualitative agreement. Another example is for Cu₂O,^100–102 where different calculations also show reasonable qualitative agreement. It appears that in general the key factors that provide the balance, such as the spread of negative mode Grüneisen parameters across reciprocal space as given by the overall network flexibility and polyhedron flexibility, are more robust, and fine details matter less.

10 Applications of the general principles

It is interesting to apply some of the ideas discussed in this paper to a relatively new example of NTE in a cubic material, NaZr₂(PO₄)₃ langbeinite (Section 2.1.12, Fig. 10). The crystal structure, with space group P2₁3, contains m = 8 ZrO₆ octahedra and n = 12 PO₄ tetrahedra in the unit cell (we assume that the Na atoms, which occupy the large cage sites in the structure, do not have a direct effect on the thermal expansion). Each octahedron is connected to 6 tetrahedra, and each tetrahedron is connected to 4 octahedra, forming sets of rings containing 2 octahedra and 2 tetrahedra. This is an over-constrained system if both tetrahedra and octahedra are rigid. However, if we consider the system to have rigid PO₄ tetrahedra and Zr–O bonds, but flexible bond angles within the ZrO₆ octahedra, the flexibility analysis discussed in Section 5.2 indicates that F = 3m + 6n = 96 and C = 6m = 48 per unit cell, giving F–C = 48, meaning there is the chance of 48 low-frequency modes for each wave vector. This is confirmed by a simple lattice dynamics calculation with appropriate forces. By contrast, from Section 5.2, ZrV₂O₇ has F–C = 24 per unit cell and ZrW₂O₈ has F–C = 36 per unit cell (both have 4 formula units per unit cell). Thus it might be expected that the langbeinite structure has the possibility to have sufficient number of low-frequency modes to support NTE.

We have developed a simple lattice model for the Zr₂(PO₄)₃ network, setting charges to be zero, and enabling short-range interactions using the Morse function (eqn (23)) for Zr–O, P–O and O–O bonds. The parameters in this minimal model were tuned to give typical frequencies for the bond-stretching and bond-bending modes in both types of polyhedra. Because this system is over-constrained, even without adding interactions to account for bending interactions centred on the oxygen atoms, the 48 flexibility modes acquire reasonable frequencies without giving a gap between them and the other modes. This point can be seen in the phonon density of states shown in Fig. 28a.

Fig. 28 (a) Phonon density of states for our minimal model of the Zr₂(PO₄)₃ langbeinite network. The first group of modes from 0 to 8 THz contain the 48 flexible modes (rigid PO₄ tetrahedra and Zr–O bonds) and 24 modes involving motions of the Zr atoms. (b) Distribution of mode Grüneisen parameters with frequency, with the colours representing the logarithmic distribution. (c) Calculated variation of the overall Grüneisen parameter and linear thermal expansivity α with temperature.

Our starting model has tetrahedral O–P–O angles within 2° of the ideal tetrahedral angle. and octahedral O–Zr–O angles within 4° of a right angle. The variance on the P–O and Zr–O bond distances is very small, but the P–O distances show stretching from the ideal distance of the potential energy function – r₀ in eqn (23) – of 4%. This gives an indication of the internal stresses experienced by the internal constraints within the network.

The distribution of mode Grüneisen parameters with frequency for our model Zr₂(PO₄)₃ network is shown in Fig. 28b. The largest contribution to NTE clearly comes from modes with frequencies between 1 and 3 THz. These are some of the modes that are flexible in the limit of rigid PO₄ tetrahedra and Zr–O bonds and flexible O–Zr–O angles. For frequencies above 8 THz, which is the regime in which we have to account for flexible PO₄ tetrahedra and Zr–O bonds, all phonon modes have positive values of their mode Grüneisen parameters.

The combination of the mode Grüneisen parameters to form the overall Grüneisen parameter is shown in Fig. 28c, together with the corresponding linear thermal expansivity. The form of the mode Grüneisen parameters, with a deep minimum at low temperature, is fairly common for an NTE material, but it is always offset by the multiplication by the heat capacity, which falls to zero at zero temperature. The overall Grüneisen parameter remains negative to high temperature, the point at which each mode contributes equally, and so the NTE continues to high temperature, albeit to a slightly lesser extent than at lower temperatures. The experimental linear thermal expansivity also remains negative to higher temperatures,¹²⁹ but its absolute value is lower than that of our calculated value. We will comment on this below.

The point of this discussion and example is that we have used the ideas articulated in this article to construct an understanding of NTE in a new material, NaZr₂(PO₄)₃.¹²⁹ We have interpreted the network in terms of rigid PO₄ tetrahedra connected to ZrO₆ octahedra with rigid Zr–O bonds. With a minimalist force model that allows for flexibility of the PO₄ tetrahedra and more rigidity of the ZrO₆ octahedra, we have been able to compute the mode Grüneisen parameters and thereby show that the NTE arises from a tension-effect phonon mechanism. The correlation with the flexibility model, in that the low-frequency modes enabled by the flexibility of the network, demonstrates the first two principles outlined in this article. The spread of modes across reciprocal space is enabled by the fact that the flexibility model has no concentration in any particular regions of reciprocal space, in contrast, say, to ScF₃.

With a minimalist model we can ask other questions of the crystal structure that are pertinent to NTE. One example is to explore a model where we replace the O–O interactions by bond-bending interactions. As discussed in the case of ScF₃, this model does not give a reasonable value of the C₁₂ elastic constant, but for exploring the structure this is not important. One interesting point about the crystal structure of NaZr₂(PO₄)₃ concerns the distribution of the four Zr–O–P angles, which are the angles that give the tension effect. Ideally this angle should be nearly linear, but in the langbeinite structure the angles are more varied, with one as low as around 140°. The lower angles are less effective in transmitting the tension effect than are linear bonds, but in the langbeinite structure the other three angles are around 20° away from linear. It can be seen that changing the relative size of the Zr–O distance compared with the P–O distance can increase all four angles, which should give a stronger tension effect.

The use of minimalist models is probably unexploited in NTE research, yet is an ideal mechanism for understanding the key principles concerning the origin of NTE, as was demonstrated in the case of ScF₃.⁸⁸ Minimalist models actually allow for more immediate contact with the key principles than do more complex models. In fact it is more common to go immediately to DFT calculations. We would like to encourage the use of minimal models first – in fact we can ask the question of whether DFT calculations are able to give much information beyond minimalist models, other than providing a closer quantitative link with experimental results? In the calculations on the Zr₂(PO₄)₃ network we have presented here, we suspect that the minimalist model is sufficient to demonstrate the origin of NTE. We look forward to seeing whether DFT calculations performed on NaZr₂(PO₄)₃ in the near future will add more information beyond what we have presented here.

We commented that the minimalist model over-estimates the size of the NTE compared to experimental value, although it does reproduce the trend towards a constant NTE at higher temperatures with larger NTE at lower temperatures. The over-estimate is probably due to the absence of the sodium cations, which will act to stiffen the lattice and hence increase the value of the bulk modulus K, which appears in the denominator of eqn (12). Does this matter? In terms of understanding, we think that probably not.

11 Conclusion

We have briefly reviewed most of the crystals of cubic symmetry that show NTE, but our purpose has not been so much to review as to draw out the key factors that underpin the thermal expansion of these materials. We have identified four principles – discussed in detail in Sections 4–7 – that allow us to understand why NTE exists in cubic crystals, but also which might limit NTE. We have illustrated these principles with pertinent examples. We have sought to go deeper into the high-level intuition about the importance of concepts such as ‘tension effect’ or ‘network flexibility’ to see how these apply in practice, and particularly with regard to the otherwise unarticulated constraints.

It is our contention that due to the efforts of the NTE community over the past three decades, we can say that we now have a detailed and general understanding of NTE in cubic crystals, such that most of the mystery has now been taken away. That is not to say that there are no more open questions. We have articulated above (Section 9) that there is still quantitative work to be done to better understand how chemical composition can change the size of, or even existence of, NTE within the same structure family. What we have provided here is the quantitative theoretical framework into which such an analysis can be woven.

Whilst one aspect of understanding the existence of a property is to explain the factors behind its existence in the known case, two other aspects are predicting behaviour in new materials and also in explaining its non-existence in other cases. Sometimes the latter case is most challenging. For example, we have highlighted the question of why NTE is found in ScF₃ but not in AlF₃ or in SrTiO₃. Our principles have allowed us to understand these systems, and they apply also to other oxides such as spinels and garnets that do not show NTE.

We can comment on the case we mentioned at the outset, SmB₆ (Section 2.1.7). Whilst this shows NTE, other members of the family, such as the alkali-earth hexaborides, do not. At first sight these look as if they can show a tension effect, but in practice the network is too rigid and the thermal expansion is dominated by the expansion of the bonds.⁹³ Although it took a calculation to demonstrate this, once it is clear that a certain network is too rigid, it follows that we can understand the absence of NTE in the whole family. The NTE in SmB₆ then arises from electronic effects.

We have also raised the challenge of whether there are more families of cubic materials showing NTE to be discovered. We have expressed some pessimism in the case of ceramic oxides, but we would be delighted to be proven wrong, especially if proven wrong several times!

There is also the possibility for new hybrid network materials with cubic symmetry. Many cyanide networks are based on the β-cristobalite or ScF₃ structures, but there are many other tetrahedral networks that could possibly be reproduced with cyanide linkers instead of oxygen. The same may be true with other ligands, following the example of the ZIF structures. The four principles we have discussed in this article should apply in these cases, but the network flexibility may be constrained by angular forces, as identified in the cyanide networks.

A decade and a bit after Cora Lind²⁰ posed the question, “[After] two decades of negative thermal expansion research: where do we stand?”, we contend that our understanding has matured considerably. In part this is due to findings associated with new materials, or more often new variants of the existing topologies, but perhaps more importantly from the maturing of supporting calculations and quantitative analysis. We hope this is reflected in this article.

It is often said by authors that one goal of NTE research is to find new materials. Perhaps we might be bold to suggest that it may be unlikely that many new materials will replace the traditional glass ceramics, which have had wide-ranging market applications for many decades. After all, many of the materials we have reviewed are toxic or unstable. But it is our conviction that our understanding of the principles associated with NTE in isotropic crystals will be invaluable in assisting the development of improved high-performance glass ceramics.

Author contributions

TL: formal analysis, investigation, visualization, writing – review and editing; MTD: conceptualization, investigation, methodology, visualization, writing – original draft.

Conflicts of interest

There are no conflicts to declare.

Data availability

No primary research results, software or code have been included and no new data were generated or analysed as part of this article.

Appendix

Appendix: Lightweight guide to the theory of lattice dynamics

The theory of lattice dynamics is well-established and articulated, to various levels of detail, in many books and articles. Mathematically it is not difficult, but the mathematics looks cumbersome because of the large number of subscripts that are required, which is further multiplied in three dimensions! Accordingly we give the theory in one dimension, written in a way such that it is relatively easy to generalise to three dimensions. Most text books will discuss one-dimensional examples, often with only 1 and 2 atoms in a unit cell, and in a way that avoids recasting the theory into professional terms, but that misses the point. So here we will imagine a one-dimensional crystal containing several atoms, such as in Fig. 29, and we will write the equations immediately in professional form. We have articulated elsewhere in pedagogical form the link between the simplified model-specific treatments and the professional way to write the equations.^216,217

Fig. 29 One dimensional crystal with four atoms in the unit cell. The horizontal line indicates the span of the unit cell.

Without worrying about issues connected with long-range electrostatic interactions – which are not a problem but make it a little more complicated that is helpful here – we can imagine that the atoms interact through pairwise forces that depend only on the distance between them. Because the crystal is periodic and can be assumed to stretch forever, we will only need to consider atoms in the unit cell at our chosen origin interacting with the atoms in neighbouring unit cells.

Let us assume that we express the energy between two atoms in terms of their separation r, where their equilibrium separation is r₀, as a Taylor expansion:

The harmonic approximation includes only the quadratic term. This can be justified on the basis that in most cases the high-order terms only give a small perturbation to the results obtained using the harmonic approximation, without significantly changing much. The first term, E(r₀), is irrelevant for the theory of lattice dynamics and so can be dropped.

We can write r − r₀ in terms of the displacements of the two atoms. Let us denote the displacement of atom j in unit cell ℓ as u_j,n, in which case we can write

r − r₀ = u_j′,n′ − u_j,n (30)

The generalisation of the harmonic energy from eqn (29) iswhereandThe latter identify can be obtained from the chore of expanding the first part of eqn (31). The factor of 1/4 in eqn (31) is the original factor of 1/2 in the harmonic expansion multiplied by another factor of 1/2 because we count every term twice in the summation.

We can quickly connect the derivatives of eqn (29) and (32). From our definition of r − r₀ above, it follows that ∂r/∂u_j′,n′ = +1 and ∂r/∂u_j,n = −1. Thus ∂E/∂u_j′,n′ = +∂E/∂r and ∂E/∂u_j,n = −∂E/∂r. For the case when the two atoms are different, that is when (j′, n′) ≠ (j, n), it follows that

However, when we consider a single atom, (j′, n′) = (j, n), we haveThese two equations can be generalised for three dimensions for radial pairwise interactions.

We note in passing that the different signs of the two derivatives can be rationalised when we consider the force constant for two atoms. The harmonic energy is

As we would expect, the energy is zero when u_j',n' = u_j,n since there is no stretching of the bond and is positive when u_j',n' = −u_j,n, that is, when the bond is stretched.

The Born–von Kármán theory of lattice dynamics assumes that we can associate the displacements of atoms at time t with travelling waves, so we can write

where A is a general amplitude for the wave, which will depend on temperature and angular frequency ω, m_j is the mass of atom j, is the relative displacement of atom j, and q is the wave vector, with magnitude given by wavelength λ as |q| = 2π/λ, where a is the lattice parameter so that the position of the origin of the unit cell is equal to R_n = na, and e_j is the relative displacement of atom j in the unit cell, such that . The approach in the Born–von Kármán theory is to write Newton's equation, force = mass × acceleration, as

When we do the derivatives we get Aexp(i(qna − ωt)) on both sides of the resultant equations; this will cancel, leaving us, after a little shuffling, with

The complex exponential function now gives a spatial phase factor for two atoms in different unit cells.

This is now looking like a matrix equation. So let us write a matrix element as

to form the matrix D(q). The reason for including the factor of in the definition of e_j should now be obvious.

We also pack the variables e_j into a column matrix

We can then rewrite eqn (39) as

ω²e = D(q)·e ⇒ ω² = e^T·D(q)·e (42)

The matrix D(q) is called the Dynamical matrix.

There is, in one dimension, as many solutions as there are numbers of atoms, N. So let us put everything into one master equation, such that we define two square matrices

andwhere the superscripts in square brackets denote the particular solution; e(q) is merely a left-to-right packing of the column arrays of eqn (41).

Eqn (42), with N different solutions, can thus be collated as

Ω(q) = e^T(q)·D(q)·e(q) (45)

There are a few points to make about eqn (45). The first is that it is a classic eigenvalue equation, such that the diagonal matrix Ω(q) is the array of eigenvalues of D(q), and e(q) is the corresponding array of eigenvectors. Note that we have the normalisation that e^T·e is the unit matrix.

The second point to make is that the dynamical matrix D(q) has the property that its complex conjugate is equal to its transpose: D(−q) = D^T(q). This is the defining characteristic of what is called a Hermitian matrix; the key property of Hermitian matrices, in the context of the theory of lattice dynamics, is that its eigenvalues, namely the components of the diagonal matrix Ω(q), are purely real, even though the Hermitian matrix itself is complex.

The point is that the individual solutions ω² for each vibration is a real number, although it can be positive or negative. If all values of ω² are positive it means that the crystal is stable in the configuration for which the calculation has been performed: it is in its lowest energy state with regard to any small deformation. On the other hand, if one value of ω² is negative, it means that the crystal is unstable against a deformation that will be exactly characterised by the corresponding column of the eigenvector matrix e. This point is very relevant for our understanding of displacive phase transitions: in the high temperature phase, a positive value of ω² for one vibration is only sustained by the perturbation from anharmonic interactions,²¹⁸ whose effect diminishes upon cooling, and the phase transition occurs at the temperature at which the value of ω² falls to zero. The crystal then undergoes a spontaneous symmetry-breaking distortion that is characterised by the corresponding form of the eigenvector e.

The mathematics can easily be generalised to three dimensions, but at the cost of adding two sets of indices to everything. The wave vector q becomes the three-dimensional vector q, and the displacements u_j,n then become the vectors u_j,n, where the 3 × 1 array n defines the neighbouring unit cell in three spatial dimensions. The distance from the origin to the lattice point R_n becomes the vector R_n. The individual elements of the force constant matrices, and , themselves become 3 × 3 sub-matrices and . Thus the dynamical matrix in three dimensions, D(q), becomes a 3N × 3N matrix, with 3N eigenvalues. Writing this out with the new pairs of subscripts merely makes the equations look more opaque, which is why we present the formalism in one dimension, and leave it to the readers' imaginations to appreciate the generalisation to three dimensions.

Acknowledgements

We thank Profs Anthony Phillips (London), John Evans (Durham) and Paul Attfield (Edinburgh) for invaluable advice in the finalisation of this manuscript. We acknowledge financial support from the National Natural Science Foundation of China, grant numbers 12174274 and 12350710177, from the Joint Fund of Bijie City and Guizhou University of Engineering Science, Bijie Science and Technology Union Contract [2026] number 2, and from Guizhou Provincial Innovation Talent Workstation (KXJZ (2025)044 – Qiankehe Platform).

References

1. T. A. Mary, J. S. O. Evans, T. Vogt and A. W. Sleight, Negative thermal expansion from 0.3 to 1050 Kelvin in ZrW₂O₈, Science, 1996, 272, 90–92,  DOI:10.1126/science.272.5258.90.
2. J. S. O. Evans, T. A. Mary, T. Vogt, M. A. Subramanian and A. W. Sleight, Negative Thermal Expansion in ZrW₂O₈ and HfW₂O₈, Chem. Mater., 1996, 8, 2809–2823,  DOI:10.1021/cm9602959.
3. C. Martinek and F. A. Hummel, Linear thermal expansion of three tungstates, J. Am. Ceram. Soc., 1968, 51, 227–228,  DOI:10.1111/j.1151-2916.1968.tb11881.x.
4. R. W. Munn, Role of the elastic constants in negative thermal expansion of axial solids, J. Phys. C-Solid State Phys., 1972, 5, 535,  DOI:10.1088/0022-3719/5/5/005.
5. T. H. K. Barron, J. A. Birch and G. K. White, Thermal expansion and heat capacity of cuprous chloride at low temperatures, J. Phys. C-Solid State Phys., 1977, 10, 1617,  DOI:10.1088/0022-3719/10/10/006.
6. J. Alamo and R. Roy, Ultralow-expansion ceramics in the system Na₂O–ZrO₂–P₂O₅–SiO₂, J. Am. Ceram. Soc., 1984, 67, c78–c80,  DOI:10.1111/j.1151-2916.1984.tb19516.x.
7. R. Roy and D. Agrawal, Thermal-expansion materials not so new, Nature, 1997, 388, 433–434,  DOI:10.1038/41243.
8. J. W. Couves, R. H. Jones, S. C. Parker, P. Tschaufeser and C. R. A. Catlow, Experimental verification of a predicted negative thermal expansivity of crystalline zeolites, J. Phys.: Condens. Matter, 1993, 5, L329–L332,  DOI:10.1088/0953-8984/5/27/001.
9. P. Tschaufeser and S. C. Parker, Thermal Expansion Behavior of Zeolites and AlPO₄s, J. Phys. Chem., 1995, 99, 10609–10615,  DOI:10.1021/j100026a026.
10. J. S. O. Evans, T. A. Mary and A. W. Sleight, Negative thermal expansion in Sc₂(WO₄)₃, J. Solid State Chem., 1998, 137, 148–160,  DOI:10.1107/S0108768198000962.
11. P. M. Forster, A. Yokochi and A. W. Sleight, Enhanced negative thermal expansion in Lu₂W₃O₁₂, J. Solid State Chem., 1998, 140, 157–158,  DOI:10.1006/jssc.1998.7967.
12. P. M. Forster and A. W. Sleight, Negative thermal expansion in Y₂W₃O₁₂, Int. J. Inorg. Mater., 1999, 1, 123–127,  DOI:10.1016/S1466-6049(99)00021-5.
13. J. S. O. Evans and T. A. Mary, Structural phase transitions and negative thermal expansion in Sc₂(MoO₄)₃, Int. J. Inorg. Mater., 2000, 2, 143–151,  DOI:10.1016/S1466-6049(00)00012-X.
14. B. A. Marinkovic, P. M. Jardim, R. R. de Avillez and F. Rizzo, Negative thermal expansion in Y₂Mo₃O₁₂, Solid State Sci., 2005, 7, 1377–1383,  DOI:10.1016/j.solidstatesciences.2005.08.012.
15. M. T. Dove and L. Li, Negative thermal expansion in cordierite, Mg₂Al₄Si₅O₁₈, Matter, 2025, 8, 101943,  DOI:10.1016/j.matt.2024.101943.
16. J. Guo, G. Cai, K. Zhao, B. Durham, K. Refson and M. T. Dove, Planar negative thermal expansion in the layered hybrid network material ZnB₂(CN)₈, Phys. Rev. B: Condens. Matter Mater. Phys., 2026, 113, 064108,  DOI:10.1103/gc3b-h3c2.
17. M. T. Dove and N. Shi, The phase transition in copper pyrophosphate, Cu₂P₂O₇: evidence from diffraction data and implications for the interpretation of negative thermal expansion, Matter, 2025, 8, 102149,  DOI:10.1016/j.matt.2025.102149.
18. J. S. O. Evans, Negative thermal expansion materials, J. Chem. Soc., Dalton Trans., 1999, 3317–3326,  10.1039/a904297k.
19. G. D. Barrera, J. A. O. Bruno, T. H. K. Barron and N. L. Allan, Negative thermal expansion, J. Phys.: Condens. Matter, 2005, 17, R217–R252,  DOI:10.1088/0953-8984/17/4/R03.
20. C. Lind, Two Decades of Negative Thermal Expansion Research: Where Do We Stand?, Materials, 2012, 5, 1125–1154,  DOI:10.3390/ma5061125.
21. J. Chen, L. Hu, J. Deng and X. Xing, Negative thermal expansion in functional materials: controllable thermal expansion by chemical modifications, Chem. Soc. Rev., 2015, 44, 3522–3567,  10.1039/c4cs00461b.
22. M. T. Dove and H. Fang, Negative thermal expansion and associated anomalous physical properties: review of the lattice dynamics theoretical foundation, Rep. Prog. Phys., 2016, 79, 066503,  DOI:10.1088/0034-4885/79/6/066503.
23. K. Takenaka, Progress of Research in Negative Thermal Expansion Materials: Paradigm Shift in the Control of Thermal Expansion, Front. Chem., 2018, 6, 267,  DOI:10.3389/fchem.2018.00267.
24. J. P. Attfield, Mechanisms and Materials for NTE, Front. Chem., 2018, 6, 371,  DOI:10.3389/fchem.2018.00371.
25. F.-X. Coudert and J. D. Evans, Nanoscale metamaterials: meta-MOFs and framework materials with anomalous behavior, Coord. Chem. Rev., 2019, 388, 48–62,  DOI:10.1016/j.ccr.2019.02.023.
26. Q. Li, H. Zhu, L. Hu, J. Chen and X. Xing, Negative Thermal Expansion in Nanosolids, Acc. Chem. Res., 2019, 52, 2694–2702,  DOI:10.1021/acs.accounts.9b00260.
27. B. Wei, Q. Sun, C. Li and J. Hong, Phonon anharmonicity: a pertinent review of recent progress and perspective, Sci. China: Phys., Mech. Astron., 2021, 64, 117001,  DOI:10.1007/s11433-021-1748-7.
28. N. Shi, Y. Song, X. Xing and J. Chen, Negative thermal expansion in framework structure materials, Coord. Chem. Rev., 2021, 449, 214204,  DOI:10.1016/j.ccr.2021.214204.
29. E. Liang, Q. Sun, H. Yuan, J. Wang, G. Zeng and Q. Gao, Negative thermal expansion: mechanisms and materials, Front. Phys., 2021, 16, 53302,  DOI:10.1007/s11467-021-1070-0.
30. J. O. Cardoso, J. P. Borges and A. Velhinho, Structural metamaterials with negative mechanical/thermomechanical indices: a review, Prog. Nat. Sci.: Mater. Int., 2021, 31, 801–808,  DOI:10.1016/j.pnsc.2021.10.015.
31. Q. Li, K. Lin, Z. Liu, L. Hu, Y. Cao, J. Chen and X. Xing, Chemical Diversity for Tailoring Negative Thermal Expansion, Chem. Rev., 2022, 122, 8438–8486,  DOI:10.1021/acs.chemrev.1c00756.
32. C. S. Coates and A. L. Goodwin, How to quantify isotropic negative thermal expansion: magnitude, range, or both?, Mater. Horiz., 2019, 6, 211–218,  10.1039/C8MH01065J.
33. R. Mittal, M. K. Gupta and S. L. Chaplot, Phonons and anomalous thermal expansion behaviour in crystalline solids, Prog. Mater. Sci., 2018, 92, 360–445,  DOI:10.1016/j.pmatsci.2017.10.002.
34. F. Hu, F. Shen, J. Hao, Y. Liu, J. Wang, J. Sun and B. Shen, Negative Thermal Expansion in the Materials With Giant Magnetocaloric Effect, Front. Chem., 2018, 6, 438,  DOI:10.3389/fchem.2018.00438.
35. Y. Song, N. Shi, S. Deng, X. Xing and J. Chen, Negative thermal expansion in magnetic materials, Prog. Mater. Sci., 2021, 121, 100835,  DOI:10.1016/j.pmatsci.2021.100835.
36. G. K. White, J. G. Collins and K. A. G. Mendelssohn, The thermal expansion of alkali halides at low temperatures - II. Sodium, rubidium and caesium halides, Proc. R. Soc. London, Ser. A, 1973, 333, 237–259,  DOI:10.1098/rspa.1973.0060.
37. K. Yang, J. Xiao, J.-W. Luo, S.-S. Li, S.-H. Wei and H.-X. Deng, A systematic study of the negative thermal expansion in zinc blende and diamond-like semiconductors, New J. Phys., 2019, 21, 123015,  DOI:10.1088/1367-2630/ab5cb3.
38. E. Bourova and P. Richet, Quartz and Cristobalite: high-temperature cell parameters and volumes of fusion, Geophys. Res. Lett., 1998, 25, 2333–2336,  DOI:10.1029/98GL01581.
39. I. P. Swainson and M. T. Dove, On the thermal expansion of β-cristobalite, Phys. Chem. Miner., 1995, 22, 61–65,  DOI:10.1007/bf00202681.
40. M. G. Tucker, M. P. Squires, M. T. Dove and D. A. Keen, Dynamic structural disorder in cristobalite: neutron total scattering measurement and reverse Monte Carlo modelling, J. Phys.: Condens. Matter, 2001, 13, 403–423,  DOI:10.1088/0953-8984/13/3/304.
41. A. P. Giddy, M. T. Dove, G. S. Pawley and V. Heine, The determination of rigid-unit modes as potential soft modes for displacive phase transitions in framework crystal structures, Acta Crystallogr., Sect. A: Found. Crystallogr., 1993, 49, 697–703,  DOI:10.1107/s0108767393002545.
42. K. D. Hammonds, M. T. Dove, A. P. Giddy, V. Heine and B. Winkler, Rigid-unit phonon modes and structural phase transitions in framework silicates, Am. Mineral., 1996, 81, 1057–1079,  DOI:10.2138/am-1996-9-1003.
43. S. A. Wells, M. T. Dove, M. G. Tucker and K. Trachenko, Real-space rigid-unit-mode analysis of dynamic disorder in quartz, cristobalite and amorphous silica, J. Phys.: Condens. Matter, 2002, 14, 4645–4657,  DOI:10.1088/0953-8984/14/18/302.
44. M. T. Dove, Flexibility of network materials and the Rigid Unit Mode model: a personal perspective, Philos. Trans. R. Soc., A, 2019, 377, 20180222,  DOI:10.1098/rsta.2018.0222.
45. L. Tan, V. Heine, G. Li and M. T. Dove, The Rigid Unit Mode model: review of ideas and applications, Rep. Prog. Phys., 2024, 876, 126501,  DOI:10.1088/1361-6633/acc7b7.
46. J. S. O. Evans, W. I. F. David and A. W. Sleight, Structural investigation of the negative-thermal-expansion material ZrW₂O₈, Acta Crystallogr., Sect. B, 1999, 55, 333–340,  DOI:10.1107/S0108768198016966.
47. W. I. F. David, J. S. O. Evans and A. W. Sleight, Direct evidence for a low-frequency phonon mode mechanism in the negative thermal expansion compound ZrW₂O₈, Europhys. Lett., 1999, 46, 661,  DOI:10.1209/epl/i1999-00316-7.
48. D. Cao, F. Bridges, G. Kowach and A. Ramirez, Frustrated soft modes and negative thermal expansion in ZrW₂O₈, Phys. Rev. Lett., 2002, 89, 215902,  DOI:10.1103/PhysRevLett.89.215902.
49. D. Cao, F. Bridges, G. R. Kowach and A. P. Ramirez, Correlated atomic motions in the negative thermal expansion material ZrW₂O₈: a local structure study, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 014303,  DOI:10.1103/PhysRevB.68.014303.
50. M. G. Tucker, A. L. Goodwin, M. T. Dove, D. A. Keen, S. A. Wells and J. S. O. Evans, Negative Thermal Expansion in ZrW₂O₈: Mechanisms, Rigid Unit Modes, and Neutron Total Scattering, Phys. Rev. Lett., 2005, 95, 255501,  DOI:10.1103/PhysRevLett.95.255501.
51. M. G. Tucker, D. A. Keen, J. S. O. Evans and M. T. Dove, Local structure in ZrW₂O₈ from neutron total scattering, J. Phys.: Condens. Matter, 2007, 19, 335215,  DOI:10.1088/0953-8984/19/33/335215.
52. A. Sanson, Toward an Understanding of the Local Origin of Negative Thermal Expansion in ZrW₂O₈: Limits and Inconsistencies of the Tent and Rigid Unit Mode Models, Chem. Mater., 2014, 26, 3716–3720,  DOI:10.1021/cm501107w.
53. F. Bridges, T. Keiber, P. Juhas, S. J. L. Billinge, L. Sutton, J. Wilde and G. R. Kowach, Local vibrations and negative thermal expansion in ZrW₂O₈, Phys. Rev. Lett., 2014, 112, 045505,  DOI:10.1103/PhysRevLett.112.045505.
54. G. Ernst, C. Broholm, G. R. Kowach and A. P. Ramirez, Phonon density of states and negative thermal expansion in ZrW₂O₈, Nature, 1998, 396, 147–149,  DOI:10.1038/24115.
55. R. Mittal, S. L. Chaplot, H. Schober and T. A. Mary, Origin of Negative Thermal Expansion in Cubic ZrW₂O₈ Revealed by High Pressure Inelastic Neutron Scattering, Phys. Rev. Lett., 2001, 86, 4692–4695,  DOI:10.1103/PhysRevLett.86.4692.
56. Y. Yamamura, N. Nakajima, T. Tsuji, M. Koyano, Y. Iwasa, S. Katayama, K. Saito and M. Sorai, Low-temperature heat capacities and Raman spectra of negative thermal expansion compounds ZrW₂O₈ and HfW₂O₈, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 014301,  DOI:10.1103/PhysRevB.66.014301.
57. J. N. Hancock, C. Turpen, Z. Schlesinger, G. R. Kowach and A. P. Ramirez, Unusual Low-Energy Phonon Dynamics in the Negative Thermal Expansion Compound ZrW₂O₈, Phys. Rev. Lett., 2004, 93, 225501,  DOI:10.1103/PhysRevLett.93.225501.
58. A. P. Ramirez and G. R. Kowach, Large Low Temperature Specific Heat in the Negative Thermal Expansion Compound ZrW₂O₈, Phys. Rev. Lett., 1998, 80, 4903–4906,  DOI:10.1103/PhysRevLett.80.4903.
59. R. Mittal and S. L. Chaplot, Lattice dynamical calculation of isotropic negative thermal expansion in ZrW₂O₈ over 0-1050 K, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 7234–7237,  DOI:10.1103/PhysRevB.60.7234.
60. K. Wang and R. R. Reeber, Mode Grüneisen parameters and negative thermal expansion of cubic ZrW₂O₈ and ZrMo₂O₈, Appl. Phys. Lett., 2000, 76, 2203–2204,  DOI:10.1063/1.126296.
61. A. Kojima, Y. Kuroiwa, S. Aoyagi, A. Sawada, Y. Yamamura, N. Nakajima and T. Tsuji, Charge density study of negative-thermal-expansion material ZrW₂O₈, J. Korean Phys. Soc., 2003, 42, 1257–1260.
62. V. Gava, A. L. Martinotto and C. A. Perottoni, First-Principles Mode Gruneisen Parameters and Negative Thermal Expansion in α-ZrW₂O₈, Phys. Rev. Lett., 2012, 109, 195503,  DOI:10.1103/PhysRevLett.109.195503.
63. M. K. Gupta, R. Mittal and S. L. Chaplot, Negative thermal expansion in cubic ZrW₂O₈: role of phonons in the entire Brillouin zone from ab initio calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 014303,  DOI:10.1103/PhysRevB.88.014303.
64. F. D. Vila, S. T. Hayashi and J. J. Rehr, Efficient Calculation of the Negative Thermal Expansion in ZrW₂O₈, Front. Chem., 2018, 6, 296,  DOI:10.3389/fchem.2018.00296.
65. L. H. N. Rimmer, M. T. Dove and K. Refson, Phonon mechanism for the negative thermal expansion of zirconium tungstate, ZrW₂O₈, Phys. Chem. Chem. Phys., 2023, 25, 564–569,  10.1039/D3CP01606D.
66. R. A. Ewings, A. I. Duff, K. Refson, T. G. Perring and J. Ollivier, Low-energy phonons in single crystal ZrW₂O₈, Phys. Rev. B: Condens. Matter Mater. Phys., 2025, 112, 014305,  DOI:10.1103/t84f-kg8t.
67. D. Taylor, Thermal-expansion data. 12. Complex oxides - AB₂O₆, AB₂O₇, A₂B₂O₇ plus complex aluminates, silicates and analogous compounds, Br. Ceram. Trans. J., 1988, 87, 39–45.
68. V. Korthuis, N. Khosrovani, A. W. Sleight, N. Roberts, R. Dupree and W. W. J. Warren, Negative thermal expansion and phase transitions in the ZrV_2−xP_xO₇7 series, Chem. Mater., 1995, 7, 412–417,  DOI:10.1021/cm00050a028.
69. N. Khosrovani, A. W. Sleight and R. Vogt, Structure of ZrV₂O₇ from −263 to 470 °C, J. Solid State Chem., 1997, 132, 355–360,  DOI:10.1006/jssc.1997.7474.
70. J. S. O. Evans, J. C. Hanson and A. W. Sleight, Room-temperature superstructure of ZrV₂O₇, Acta Crystallogr., Sect. B, 1998, 54, 705–713,  DOI:10.1107/S0108768198000962.
71. T. Chatterji, P. Henry, R. Mittal and S. Chaplot, Negative thermal expansion of ReO₃: neutron diffraction experiments and dynamical lattice calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134105,  DOI:10.1103/PhysRevB.78.134105.
72. T. Chatterji, T. C. Hansen, M. Brunelli and P. F. Henry, Negative thermal expansion of ReO₃ in the extended temperature range, Appl. Phys. Lett., 2009, 94, 241902,  DOI:10.1063/1.3155191.
73. M. Dapiaggi and A. N. Fitch, Negative (and very low) thermal expansion in ReO₃ from 5 to 300 K, J. Appl. Crystallogr., 2009, 42, 253–258,  DOI:10.1107/S002188980804332X.
74. B. K. Greve, K. L. Martin, P. L. Lee, P. J. Chupas, K. W. Chapman and A. P. Wilkinson, Pronounced negative thermal expansion from a simple structure: cubic ScF₃, J. Am. Chem. Soc., 2010, 132, 15496–15498,  DOI:10.1021/ja106711v.
75. T. A. Bird, M. G. L. Wilkinson, D. A. Keen, R. I. Smith, N. C. Bristowe, M. T. Dove, A. E. Phillips and M. S. Senn, Soft-mode anisotropy in the negative thermal expansion material ReO₃, Phys. Rev. B: Condens. Matter Mater. Phys., 2021, 104, 214102,  DOI:10.1103/physrevb.101.064306.
76. J. C. Hancock, K. W. Chapman, G. J. Halder, C. R. Morelock, B. S. Kaplan, L. C. Gallington, A. Bongiorno, C. Han, S. Zhou and A. P. Wilkinson, Large Negative Thermal Expansion and Anomalous Behavior on Compression in Cubic ReO₃-type A(II)B(IV)F₆: CaZrF₆ and CaHfF₆, Chem. Mater., 2015, 27, 3912–3918,  DOI:10.1021/acs.chemmater.5b00662.
77. L. Hu, J. Chen, J. Xu, N. Wang, F. Han, Y. Ren, Z. Pan, Y. Rong, R. Huang and J. Deng, Atomic Linkage Flexibility Tuned Isotropic Negative, Zero, and Positive Thermal Expansion in MZrF₆ (M = Ca, Mn, Fe, Co, Ni, and Zn), J. Am. Chem. Soc., 2016, 138, 14530–14533,  DOI:10.1021/jacs.6b08746.
78. F. Han, L. Hu, Z. Liu, Q. Li, T. Wang, Y. Ren, J. Deng, J. Chen and X. Xing, Local structure and controllable thermal expansion in the solid solution (Mn_1−xNi_x)ZrF₆, Inorg. Chem. Front., 2017, 4, 343–347,  10.1039/C6QI00483K.
79. L. Hu, J. Chen, A. Sanson, H. Wu, C. G. Rodriguez, L. Olivi, Y. Ren, L. Fan, J. Deng and X. Xing, New Insights into the Negative Thermal Expansion: Direct Experimental Evidence for the “Guitar-String” Effect in Cubic ScF₃, J. Am. Chem. Soc., 2016, 138, 8320–8323,  DOI:10.1021/jacs.6b02370.
80. D. Wendt, E. Bozin, J. Neuefeind, K. Page, W. Ku, L. Wang, B. Fultz, A. V. Tkachenko and I. A. Zaliznyak, Entropic elasticity and negative thermal expansion in a simple cubic crystal, Sci. Adv., 2019, 5, eaay2748,  DOI:10.1126/sciadv.aay2748.
81. M. T. Dove, J. Du, Z. Wei, D. A. Keen, M. G. Tucker and A. E. Phillips, Quantitative understanding of negative thermal expansion in scandium trifluoride from neutron total scattering measurements, Phys. Rev. B: Condens. Matter Mater. Phys., 2020, 094105,  DOI:10.1103/physrevb.102.094105.
82. S. U. Handunkanda, E. B. Curry, V. Voronov, A. H. Said, G. G. Guzmán-Verri, R. T. Brierley, P. B. Littlewood and J. N. Hancock, Large isotropic negative thermal expansion above a structural quantum phase transition, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 134101–134106,  DOI:10.1103/physrevb.92.134101.
83. S. U. Handunkanda, C. A. Occhialini, A. H. Said and J. N. Hancock, Two-dimensional nanoscale correlations in the strong negative thermal expansion material ScF₃, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 94, 214102,  DOI:10.1103/physrevb.94.214102.
84. C. A. Occhialini, G. G. Guzmán-Verri, S. U. Handunkanda and J. N. Hancock, Negative thermal expansion near the precipice of structural stability in open perovskites, Front. Chem., 2018, 6, 564,  DOI:10.3389/fchem.2018.00545.
85. Y. Liu, Z. Wang, M. Wu, Q. Sun, M. Chao and Y. Jia, Negative thermal expansion in isostructural cubic ReO₃ and ScF₃: a comparative study, Comput. Mater. Sci., 2015, 107, 157–162,  DOI:10.1016/j.commatsci.2015.05.019.
86. C. W. Li, X. Tang, J. A. Muñoz, J. B. Keith, S. J. Tracy, D. L. Abernathy and B. Fultz, Structural relationship between negative thermal expansion and quartic anharmonicity of cubic ScF₃, Phys. Rev. Lett., 2011, 107, 195504,  DOI:10.1103/PhysRevLett.107.195504.
87. Y. Oba, T. Tadano, R. Akashi and S. Tsuneyuki, First-principles study of phonon anharmonicity and negative thermal expansion in ScF₃, Phys. Rev. Mater., 2019, 3, 033601,  DOI:10.1103/PhysRevMaterials.3.033601.
88. M. T. Dove, Z. Wei, A. E. Phillips, D. A. Keen and K. Refson, Which phonons contribute most to negative thermal expansion in ScF₃?, APL Mater., 2023, 11, 041130,  DOI:10.1063/5.0147610.
89. N. Wang, J. Deng, J. Chen and X. Xing, Phonon spectrum attributes for the negative thermal expansion of MZrF₆ (M = Ca, Mn–Ni, Zn), Inorg. Chem. Front., 2019, 6, 1022–1028,  10.1039/C8QI01176A.
90. K. Zhao, Y. Jiao, Q. Sun, R. He, A. Sanson, Z. Zhong, E. Liang and Q. Gao, Octahedral tilt distortion in negative thermal expansion in the fluorides CaZrF₆ and ScF₃, Phys. Rev. B: Condens. Matter Mater. Phys., 2024, 110, 064322,  DOI:10.1103/PhysRevB.110.064322.
91. C. R. Morelock, J. C. Hancock and A. P. Wilkinson, Thermal expansion and phase transitions of α-AlF₃, J. Solid State Chem., 2014, 219, 143–147,  DOI:10.1016/j.jssc.2014.07.031.
92. V. A. Trounov, A. L. Malyshev, D. Y. Chernyshov, M. M. Korsukova, V. N. Gurin, L. A. Aslanov and V. V. Chernyshev, Temperature dependences of the parameters of atoms in the crystal structure of the intermediate-valence semiconductor SmB₆: investigation by high-resolution powder neutron diffraction, J. Phys.: Condens. Matter, 1993, 5, 2479,  DOI:10.1088/0953-8984/5/16/007.
93. L. Li, M. T. Dove, Z. Wei, A. E. Phillips and D. S. Keeble, Electronic origin of negative thermal expansion in samarium hexaboride revealed by X-ray diffraction and total scattering, Phys. Chem. Chem. Phys., 2024, 26, 7664–7673,  10.1039/D3CP05954E.
94. L. Li, K. Refson and M. T. Dove, The contribution of phonons to the thermal expansion of some simple cubic hexaboride structures: SmB₆, CaB₆, SrB₆ and BaB₆, Phys. Chem. Chem. Phys., 2023, 25, 10749–10758,  10.1039/D3CP01306E.
95. W. Tiano, M. Dapiaggi and G. Artioli, Thermal expansion in cuprite-type structures from 10 K to decomposition temperature: Cu₂O and Ag₂O, J. Appl. Crystallogr., 2003, 36, 1461–1463,  DOI:10.1107/S0021889803020818.
96. M. Dapiaggi, W. Tiano, G. Artioli, A. Sanson and P. Fornasini, The thermal behaviour of cuprite: an XRD-EXAFS combined approach, Nucl. Instrum. Methods Phys. Res., Sect. B, 2003, 200, 231–236,  DOI:10.1016/S0168-583X(02)01682-8.
97. G. Artioli, M. Dapiaggi, P. Fornasini, A. Sanson, F. Rocca and M. Merli, Negative thermal expansion in cuprite-type compounds: a combined synchrotron XRPD, EXAFS, and computational study of Cu₂O and Ag₂O, J. Phys. Chem. Solids, 2006, 67, 1918–1922,  DOI:10.1016/j.jpcs.2006.05.043.
98. M. Dapiaggi, H. Kim, E. S. Božin, S. J. L. Billinge and G. Artioli, Study of the negative thermal expansion of cuprite-type structures by means of temperature-dependent pair distribution function analysis: preliminary results, J. Phys. Chem. Solids, 2008, 69, 2182–2186,  DOI:10.1016/j.jpcs.2008.03.030.
99. K. W. Chapman and P. J. Chupas, Anomalous thermal expansion of cuprites: a combined high resolution pair distribution function and geometric analysis, Chem. Mater., 2009, 21, 425–431,  DOI:10.1021/cm802900t.
100. K.-P. Bohnen, R. Heid, L. Pintschovius, A. Soon and C. Stampfl, Ab initio lattice dynamics and thermal expansion of Cu₂O, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 134304,  DOI:10.1103/PhysRevB.80.134304.
101. M. K. Gupta, R. Mittal, S. L. Chaplot and S. Rols, Phonons, nature of bonding, and their relation to anomalous thermal expansion behavior of M₂O (M = Au, Ag, Cu), J. Appl. Phys., 2014, 115, 093507,  DOI:10.1063/1.4867437.
102. L. H. N. Rimmer, M. T. Dove, B. Winkler, D. J. Wilson, K. Refson and A. L. Goodwin, Framework flexibility and the negative thermal expansion mechanism of copper(I) oxide Cu₂O, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 214115,  DOI:10.1103/physrevb.89.214115.
103. A. L. Goodwin and C. Kepert, Negative thermal expansion and low-frequency modes in cyanide-bridged framework materials, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 140301,  DOI:10.1103/physrevb.71.140301.
104. D. J. Williams, D. E. Partin, F. J. Lincoln, J. Kouvetakis and M. O’Keeffe, The Disordered Crystal Structures of Zn(CN)₂ and Ga(CN)₃, J. Solid State Chem., 1997, 134, 164–169,  DOI:10.1006/jssc.1997.7571.
105. K. W. Chapman, P. J. Chupas and C. J. Kepert, Direct observation of a transverse vibrational mechanism for negative thermal expansion in Zn(CN)₂: an atomic pair distribution function analysis, J. Am. Chem. Soc., 2005, 127, 15630–15636,  DOI:10.1021/ja055197f.
106. J. W. Zwanziger, Phonon dispersion and Grüneisen parameters of zinc dicyanide and cadmium dicyanide from first principles: origin of negative thermal expansion, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 052102,  DOI:10.1103/PhysRevB.76.052102.
107. H. Fang, M. T. Dove, L. H. N. Rimmer and A. J. Misquitta, Simulation study of pressure and temperature dependence of the negative thermal expansion in Zn(CN)₂, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 104306,  DOI:10.1103/physrevb.88.104306.
108. A. Phillips, A. Goodwin, G. Halder, P. Southon and C. Kepert, Nanoporosity and Exceptional Negative Thermal Expansion in Single-Network Cadmium Cyanide, Angew. Chem., Int. Ed., 2008, 47, 1396–1399,  DOI:10.1002/anie.200704421.
109. Q. Gao, J. Wang, A. Sanson, Q. Sun, E. Liang, X. Xing and J. Chen, Discovering Large Isotropic Negative Thermal Expansion in Framework Compound AgB(CN)₄ via the Concept of Average Atomic Volume, J. Am. Chem. Soc., 2020, 142, 6935–6939,  DOI:10.1021/jacs.0c02188.
110. Q. Gao, Y. Jiao, Q. Sun, J. A. P. Sprenger, M. Finze, A. Sanson, E. Liang, X. Xing and J. Chen, Giant Negative Thermal Expansion in Ultralight [NaB(CN)₄], Angew. Chem., Int. Ed., 2024, 63, e202401302,  DOI:10.1002/anie.202401302.
111. L. Li, K. Refson and M. T. Dove, Negative thermal expansion of cubic silicon dicarbodiimide, Si(NCN)₂, studied by ab initio lattice dynamics, J. Phys.: Condens. Matter, 2020, 32, 465402,  DOI:10.1088/1361-648x/aba8cb.
112. P. Kroll, M. Andrade, X. Yan, E. Ionescu, G. Miehe and R. Riedel, Isotropic Negative Thermal Expansion in β-Si(NCN)₂ and Its Origin, J. Phys. Chem. C, 2012, 116, 526–531,  DOI:10.1021/jp2106583.
113. S. Adak, L. L. Daemen, M. Hartl, D. Williams, J. Summerhill and H. Nakotte, Thermal expansion in 3d-metal Prussian blue analogs – A survey study, J. Solid State Chem., 2011, 184, 2854–2861,  DOI:10.1016/j.jssc.2011.08.030.
114. Q. Gao, N. Shi, Q. Sun, A. Sanson, R. Milazzo, A. Carnera, H. Zhu, S. H. Lapidus, Y. Ren, Q. Huang, J. Chen and X. Xing, Low-Frequency Phonon Driven Negative Thermal Expansion in Cubic GaFe(CN)₆ Prussian Blue Analogues, Inorg. Chem., 2018, 57, 10918–10924,  DOI:10.1021/acs.inorgchem.8b01526.
115. N. Shi, Q. Gao, A. Sanson, Q. Li, L. Fan, Y. Ren, L. Olivi, J. Chen and X. Xing, Negative thermal expansion in cubic FeFe(CN)₆ Prussian blue analogues, Dalton Trans., 2019, 48, 3658–3663,  10.1039/C8DT05111A.
116. Q. Gao, E. Liang, X. Xing and J. Chen, Negative Thermal Expansion in Prussian Blue Analogues, Chem. J. Chin. Univ., 2020, 41, 388–400,  DOI:10.7503/cjcu20190631.
117. S. Park, R.-W. Große Kunstleve, H. Graetsch and H. Gies, in Progress in Zeolite and Microporous Materials, ed. H. Chon, S.-K. Ihm and Y. S. Uh, Elsevier, 1997, vol. 107 of Studies in Surface Science and Catalysis, 1989–1994 DOI:10.1016/S0167-2991(97)80664-7.
118. D. A. Woodcock, P. Lightfoot, P. A. Wright, L. A. Villaescusa and M. A. Camblor, Strong negative thermal expansion in the siliceous zeolites ITQ-1, ITQ-3 and SSZ-23, J. Mater. Chem., 1999, 9, 349–351,  10.1039/a808059c.
119. P. Lightfoot, D. A. Woodcock, M. J. Maple, L. A. Villaescusa and P. A. Wright, The widespread occurrence of negative thermal expansion in zeolites, J. Mater. Chem., 2001, 11, 212–216,  10.1039/B002950P.
120. C. Baerlocher and L. McCusker, Database of Zeolite Structures, https://www.iza-structure.org/databases/, 2025, Accessed: 19 February 2026.
121. H. Fang and M. T. Dove, Pressure-induced softening as a common feature of framework structures with negative thermal expansion, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 214109,  DOI:10.1103/physrevb.87.214109.
122. M. Ducamp and F.-X. Coudert, Systematic Study of the Thermal Properties of Zeolitic Frameworks, J. Phys. Chem. C, 2021, 125, 15647–15658,  DOI:10.1021/acs.jpcc.1c03975.
123. T. Carey, A. Corma, F. Rey, C. C. Tang, J. A. Hriljac and P. A. Anderson, The effect of extra framework species on the intrinsic negative thermal expansion property of zeolites with the LTA topology, Chem. Commun., 2012, 48, 5829–5831,  10.1039/C2CC30582H.
124. T. Carey, C. C. Tang, J. A. Hriljac and P. A. Anderson, Chemical Control of Thermal Expansion in Cation-Exchanged Zeolite A, Chem. Mater., 2014, 26, 1561–1566,  DOI:10.1021/cm403312q.
125. M. P. Attfield and A. W. Sleight, Strong negative thermal expansion in siliceous faujasite, Chem. Commun., 1998, 601–602,  10.1039/a707141h.
126. A. Dosen and B. Marinkovic, Negative thermal expansion and cationic migration in zeolite Y used in FCC catalysts, Bull. Mater. Sci., 2019, 42, 86,  DOI:10.1007/s12034-019-1778-4.
127. A. Shablinskii, S. Demina, Y. Biryukov, R. Bubnova, M. Krzhizhanovskaya and S. Filatov, A structural origin of both positive and negative thermal expansion in langbeinite-, arcanite- and metathenardite-type and related Rb₂SO₄ and Rb₂Ca₂(SO₄)₃ compounds, Ceram. Int., 2025, 51, 51342–51350,  DOI:10.1016/j.ceramint.2025.08.360.
128. M. Kahrizi and M. Steinitz, Phase transitions and thermal expansion in langbeinite type compounds, Solid State Commun., 1988, 66, 375–378,  DOI:10.1016/0038-1098(88)90860-5.
129. S. N. Marshenya, A. G. Scherbakov, A. D. Dembitskiy, A. A. Golubnichiy, I. A. Trussov, A. A. Savina, S. M. Kazakov, D. A. Aksyonov, E. V. Antipov and S. S. Fedotov, NaZr₂(PO₄)₃ – a cubic langbeinite-type sodium-ion solid conductor, Dalton Trans., 2024, 53, 15928–15936,  10.1039/D4DT02288B.
130. H. Wu, T. Yildirim, J. R. Simpson and A. R. H. Walker, Origin of the exceptional negative thermal expansion in metal-organic framework-5 Zn₄O(1,4-benzenedicarboxylate)₃, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 054114,  DOI:10.1103/PhysRevB.78.054114.
131. N. Lock, Y. Wu, M. Christensen, L. J. Cameron, V. K. Peterson, A. J. Bridgeman, C. J. Kepert and B. B. Iversen, Elucidating negative thermal expansion in MOF-5, J. Phys. Chem. C, 2010, 114, 16181–16186,  DOI:10.1021/jp103212z.
132. N. Lock, M. Christensen, Y. Wu, V. K. Peterson, M. K. Thomsen, R. O. Piltz, A. J. Ramirez-Cuesta, G. J. McIntyre, K. Norén, R. Kutteh, C. J. Kepert, G. J. Kearley and B. B. Iversen, Scrutinizing negative thermal expansion in MOF-5 by scattering techniques and ab initio calculations, Dalton Trans., 2013, 42, 1996–2007,  10.1039/C2DT31491F.
133. L. H. N. Rimmer, M. T. Dove, A. L. Goodwin and D. C. Palmer, Acoustic phonons and negative thermal expansion in MOF-5, Phys. Chem. Chem. Phys., 2014, 16, 21144–21152,  10.1039/C4CP01701C.
134. M. J. Cliffe, J. A. Hill, C. A. Murray, F.-X. Coudert and A. L. Goodwin, Defect-dependent colossal negative thermal expansion in UiO-66(Hf) metal–organic framework, Phys. Chem. Chem. Phys., 2015, 17, 11586–11592,  10.1039/C5CP01307K.
135. S. M. Vornholt, Z. Chen, J. Hofmann and K. W. Chapman, Node Distortions in UiO-66 Inform Negative Thermal Expansion Mechanisms: Kinetic Effects, Frustration, and Lattice Hysteresis, J. Am. Chem. Soc., 2024, 146, 16977–16981,  DOI:10.1021/jacs.4c05313.
136. C. Chen, H. E. Maynard-Casely, S. G. Duyker, R. Babarao, C. J. Kepert, J. D. Evans and L. K. Macreadie, Lowering the Energetic Landscape for Negative Thermal Expansion in 3D-Linker Metal–Organic Frameworks, Chem. Mater., 2023, 35, 9945–9951,  DOI:10.1021/acs.chemmater.3c01744.
137. Y. Wu, A. Kobayashi, G. Halder, V. Peterson, K. Chapman, N. Lock, P. Southon and C. Kepert, Negative Thermal Expansion in the Metal–Organic Framework Material Cu₃ (1,3,5-benzenetricarboxylate)₂, Angew. Chem., Int. Ed., 2008, 47, 8929–8932,  DOI:10.1002/anie.200803925.
138. C. Schneider, D. Bodesheim, M. G. Ehrenreich, V. Crocellà, J. Mink, R. A. Fischer, K. T. Butler and G. Kieslich, Tuning the Negative Thermal Expansion Behavior of the Metal–Organic Framework Cu₃BTC₂ by Retrofitting, J. Am. Chem. Soc., 2019, 141, 10504–10509,  DOI:10.1021/jacs.9b04755.
139. A. F. Sapnik, H. S. Geddes, E. M. Reynolds, H. H.-M. Yeung and A. L. Goodwin, Compositional inhomogeneity and tuneable thermal expansion in mixed-metal ZIF-8 analogues, Chem. Commun., 2018, 54, 9651–9654,  10.1039/C8CC04172E.
140. E. Grüneisen, Theorie des festen Zustandes einatomiger Elemente, Ann. Phys., 1912, 39, 257–306,  DOI:10.1002/andp.19123441202.
141. E. Grüneisen, in Zustand des festen Körpers, ed. C. Drucker, E. Grüneisen, P. Kohnstamm, F. Körber, K. Scheel, E. Schrödinger, F. Simon, J. D. van der Waals and F. Henning, Springer Berlin Heidelberg, Berlin, Heidelberg, 1926, pp. 1–59 DOI:10.1007/978-3-642-99531-6_1.
142. E. Grüneisen, The state of a solid body, NASA Tech. Rep., 1959, RE 2-18-59W, 1–76, https://ntrs.nasa.gov/api/citations/19980231989/downloads/19980231989.pdf.
143. T. H. K. Barron, Grüneisen parameters for the equation of state of solids, Ann. Phys., 1957, 1, 77–90,  DOI:10.1016/0003-4916(57)90006-4.
144. M. G. Tucker, M. T. Dove and D. A. Keen, Direct measurement of the thermal expansion of the Si–O bond by neutron total scattering, J. Phys.: Condens. Matter, 2000, 12, L425–L430,  DOI:10.1088/0953-8984/12/26/101.
145. Q. Hui, M. G. Tucker, M. T. Dove, S. A. Wells and D. A. Keen, Total scattering and reverse Monte Carlo study of the 105 K displacive phase transition in strontium titanate, J. Phys.: Condens. Matter, 2005, 17, S111–S124,  DOI:10.1088/0953-8984/17/5/012.
146. R. W. Robinett, Average value of position for the anharmonic oscillator: classical versus quantum results, Am. J. Phys., 1997, 65, 190–194,  DOI:10.1119/1.18747.
147. T. Claeys, B. Fahs, G. Lambert and C. Webb, How much can the eigenvalues of a random Hermitian matrix fluctuate?, Duke Math. J., 2021, 170, 2085–2235,  DOI:10.1215/00127094-2020-0070.
148. V. Heine, P. R. L. Welche and M. T. Dove, Geometrical Origin and Theory of Negative Thermal Expansion in Framework Structures, J. Am. Ceram. Soc., 1999, 82, 1793–1802,  DOI:10.1111/j.1151-2916.1999.tb02001.x.
149. T. H. K. Barron, Vibrational effects in the thermal expansion of noncubic solids, J. Appl. Phys., 1970, 41, 5044,  DOI:10.1063/1.1658595.
150. H. Fang, M. T. Dove and A. E. Phillips, Common origin of negative thermal expansion and other exotic properties in ceramic and hybrid materials, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 214103,  DOI:10.1103/physrevb.89.214103.
151. H. Fang, A. E. Phillips, M. T. Dove, M. G. Tucker and A. L. Goodwin, Temperature-dependent pressure-induced softening in Zn(CN)₂, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 144103,  DOI:10.1103/physrevb.88.144103.
152. Z. Wei, L. Tan, G. Cai, A. E. Phillips, I. da Silva, M. G. Kibble and M. T. Dove, Colossal Pressure-Induced Softening in Scandium Fluoride, Phys. Rev. Lett., 2020, 124, 255502,  DOI:10.1103/PhysRevLett.124.255502.
153. M. T. Dove, M. Gambhir, K. D. Hammonds, V. Heine and A. K. A. Pryde, Distortions of framework structures, Phase Trans., 1996, 58, 121–143,  DOI:10.1080/01411599608242398.
154. A. L. Goodwin, Rigid unit modes and intrinsic flexibility in linearly bridged framework structures, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 134302,  DOI:10.1103/physrevb.74.134302.
155. S. J. Hibble, G. B. Wood, E. J. Bilbé, A. H. Pohl, M. G. Tucker, A. C. Hannon and A. M. Chippindale, Structures and negative thermal expansion properties of the one-dimensional cyanides, CuCN, AgCN and AuCN, Z. Kristallogr., 2010, 225, 457–462,  DOI:10.1524/zkri.2010.1314.
156. T. Chatterji, M. Zbiri and T. C. Hansen, Negative thermal expansion in ZnF₂, Appl. Phys. Lett., 2011, 98, 181911,  DOI:10.1063/1.3588414.
157. G. Cai, Z. Wei, A. E. Phillips and M. T. Dove, Manuscript in preparation.
158. J. P. Stoppelman, A. P. Wilkinson and J. G. McDaniel, Equation of state predictions for ScF₃ and CaZrF₆ with neural network-driven molecular dynamics, J. Chem. Phys., 2023, 159, 084707,  DOI:10.1063/5.0157615.
159. M. K. Gupta, B. Singh, R. Mittal and S. L. Chaplot, Negative thermal expansion behavior in MZrF₆ (M = Ca, Mg, Sr): ab initio lattice dynamical studies, Phys. Rev. B: Condens. Matter Mater. Phys., 2018, 98, 014301,  DOI:10.1103/physrevb.98.014301.
160. J. C. Maxwell, On the calculation of the equilibrium and stiffness of frames, London, Edinburgh Dublin Philos. Mag. J. Sci., 1864, 27, 294–299,  DOI:10.1080/14786446408643668.
161. J. C. Phillips, Topology of covalent non-crystalline solids I: short-range order in chalcogenide alloys, J. Non-Cryst. Solids, 1979, 34, 153–181,  DOI:10.1016/0022-3093(79)90033-4.
162. J. C. Phillips, The physics of glass, Phys. Today, 1982, 35, 27–33,  DOI:10.1063/1.2914932.
163. M. F. Thorpe, Continuous deformations in random networks, J. Non-Cryst. Solids, 1983, 57, 355–370,  DOI:10.1016/0022-3093(83)90424-6.
164. H. He and M. F. Thorpe, Elastic properties of glasses, Phys. Rev. Lett., 1985, 54, 2107–2110,  DOI:10.1103/physrevLetters54.2107.
165. Y. Cai and M. F. Thorpe, Floppy modes in network glasses, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 10535–10542,  DOI:10.1103/physrevb.40.10535.
166. M. F. Thorpe, M. V. Chubynsky, D. J. Jacobs and J. C. Phillips, Non-randomness in network glasses and rigidity, Glass Phys. Chem., 2001, 27, 160–166,  DOI:10.1023/A:1011336511583.
167. K. D. Hammonds, A. Bosenick, M. T. Dove and V. Heine, Rigid unit modes in crystal structures with octahedrally coordinated atoms, Am. Mineral., 1998, 83, 476–479,  DOI:10.2138/am-1998-5-607.
168. L. H. N. Rimmer and M. T. Dove, Simulation study of negative thermal expansion in yttrium tungstate Y₂W₃O₁₂, J. Phys.: Condens. Matter, 2015, 185401,  DOI:10.1088/0953-8984/27/18/185401.
169. A. K. A. Pryde, K. D. Hammonds, M. T. Dove, V. Heine, J. D. Gale and M. C. Warren, Origin of the negative thermal expansion in ZrW₂O₈ and ZrV₂O₇, J. Phys.: Condens. Matter, 1996, 8, 10973–10982,  DOI:10.1088/0953-8984/8/50/023.
170. A. Sanson, F. Rocca, G. Dalba, P. Fornasini, R. Grisenti, M. Dapiaggi and G. Artioli, Negative thermal expansion and local dynamics in Cu₂O and Ag₂O, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 214305,  DOI:10.1103/PhysRevB.73.214305.
171. M. T. Dove, A. K. A. Pryde, V. Heine and K. D. Hammonds, Exotic distributions of rigid unit modes in the reciprocal spaces of framework aluminosilicates, J. Phys.: Condens. Matter, 2007, 19, 275209,  DOI:10.1088/0953-8984/19/27/275209.
172. K. Hammonds, H. Deng, V. Heine and M. T. Dove, How floppy modes give rise to adsorption sites in zeolites, Phys. Rev. Lett., 1997, 78, 3701–3704,  DOI:10.1103/PhysRevLett.78.3701.
173. K. D. Hammonds, V. Heine and M. T. Dove, Rigid-unit modes and the quantitative determination of the flexibility possessed by zeolite frameworks, J. Phys. Chem. B, 1998, 102, 1759–1767,  DOI:10.1021/jp980006z.
174. A. van Roekeghem, J. Carrete and N. Mingo, Anomalous thermal conductivity and suppression of negative thermal expansion in ScF₃, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 94, 020303–020305,  DOI:10.1103/PhysRevB.94.020303.
175. M. Liu, T. R. Finlayson and T. F. Smith, High-resolution dilatometry measurements of SrTiO₃ along cubic and tetragonal axes, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 3480–3484,  DOI:10.1103/PhysRevB.55.3480.
176. A. Okazaki and M. Kawaminami, Lattice constant of strontium titanate at low temperatures, Mater. Res. Bull., 1973, 8, 545–550,  DOI:10.1016/0025-5408(73)90130-x.
177. A. Okazaki and M. Kawaminami, Lattice constant measurements on the 106 K transition of SrTiO₃, Ferroelectrics, 1974, 7, 91–92,  DOI:10.1080/00150197408237957.
178. W. G. Stirling, Neutron inelastic scattering study of the lattice dynamics of strontium titanate: harmonic models, J. Phys. C: Solid State Phys., 1972, 5, 2711–2730,  DOI:10.1088/0022-3719/5/19/005.
179. C. R. Morelock, L. C. Gallington and A. P. Wilkinson, Solid solubility, phase transitions, thermal expansion, and compressibility in Sc_1−xAl_xF₃, J. Solid State Chem., 2015, 222, 96–102,  DOI:10.1016/j.jssc.2014.11.007.
180. C. R. Morelock, L. C. Gallington and A. P. Wilkinson, Evolution of Negative Thermal Expansion and Phase Transitions in Sc_1−xTi_xF₃, Chem. Mater., 2014, 26, 1936–1940,  DOI:10.1021/cm5002048.
181. Y. Qiao, S. Zhang, P. Zhang, J. Guo, A. Sanson, X. Zhen, K. Zhao, Q. Gao and J. Chen, Simple chemical synthesis and isotropic negative thermal expansion in MHfF₆ (M = Ca, Mn, Fe and Co), Nano Res., 2024, 17, 2195–2203,  DOI:10.1007/s12274-024-6445-4.
182. E. Bernhardt, G. Henkel and H. Willner, Die Tetracyanoborate M[B(CN)₄], M = [Bu₄ N]⁺, Ag⁺, K⁺, Z. Anorg. Allg. Chem., 2000, 626, 560–568,  DOI:10.1002/(SICI)1521-3749(200002)626:2<560::AID-ZAAC560>3.0.CO;2-E.
183. J. Guo and M. T. Dove, Unpublished calculations, more details are available from the corresponding author on request.
184. M. T. Dove, Unpublished calculations, more details are available from the corresponding author on request.
185. S. A. Wells, M. T. Dove and M. G. Tucker, Finding best-fit polyhedral rotations with geometric algebra, J. Phys.: Condens. Matter, 2002, 14, 4567–4584,  DOI:10.1088/0953-8984/14/17/327.
186. S. A. Wells, M. T. Dove and M. G. Tucker, Reverse Monte Carlo with geometric analysis – RMC + GA, J. Appl. Crystallogr., 2004, 37, 536–544,  DOI:10.1107/s0021889804008957.
187. S. A. Wells and A. Sartbaeva, GASP: software for geometric simulations of flexibility in polyhedral and molecular framework structures, Mol. Simul., 2015, 41, 1409–1421,  DOI:10.1080/08927022.2015.1032277.
188. J. Liu, J. Du, P. B. Wyatt, D. A. Keen, A. E. Phillips and M. T. Dove, Orientational order/disorder and network flexibility in deuterated methylammonium lead iodide perovskite by neutron total scattering, J. Mater. Chem. A, 2024, 12, 2771–2785,  10.1039/D3TA04586B.
189. J. Liu, A. E. Phillips, D. A. Keen and M. T. Dove, Thermal Disorder and Bond Anharmonicity in Cesium Lead Iodide Studied by Neutron Total Scattering and the Reverse Monte Carlo Method, J. Phys. Chem. C, 2019, 123, 14934–14940,  DOI:10.1021/acs.jpcc.9b02936.
190. M. T. Dove and G. Li, Review: pair distribution functions from neutron total scattering for the study of local structure in disordered materials, Nuclear Anal., 2022, 1, 100037,  DOI:10.1016/j.nucana.2022.100037.
191. D. A. Keen, M. G. Tucker and M. T. Dove, Reverse Monte Carlo modelling of crystalline disorder, J. Phys.: Condens. Matter, 2005, 17, S15,  DOI:10.1088/0953-8984/17/5/002.
192. H. Koiso, S. Yoshida, T. Nagai, T. Isobe, A. Nakajima and Y. Mochizuki, Thermal expansion and phase stability of BF₃ (B = Sc, Y, La, Al, Ga, In) from first principles, Phys. Rev. B, 2024, 110, 064104,  DOI:10.1103/PhysRevB.110.064104.
193. J. P. Rumsey, The secret behind crystals that shrink when heated, https://www.ornl.gov/news/secret-behind-crystals-shrink-when-heated, 2019, Accessed: 2026-02-17.
194. P. Li, P. Yao, G. Cai, T. Li, A. E. Phillips and M. T. Dove, Submitted for publication.
195. W. W. Schmahl, I. P. Swainson, M. T. Dove and A. Graeme-Barber, Landau free energy and order parameter behaviour of the α/β phase transition in cristobalite, Z. Kristallogr., 1992, 201, 125–145,  DOI:10.1524/zkri.1992.201.1-2.125.
196. I. Swainson and M. T. Dove, Low-frequency floppy modes in β-cristobalite, Phys. Rev. Lett., 1993, 71, 193–196,  DOI:10.1103/physrevLetters71.193.
197. G. L. Hua, T. R. Welberry, R. L. Withers and J. G. Thompson, An electron diffraction and lattice-dynamical study of the diffuse scattering in β-cristobalite, SiO₂, J. Appl. Crystallogr., 1988, 21, 458–465,  DOI:10.1107/s0021889888004637.
198. I. P. Swainson and M. T. Dove, Molecular dynamics simulation of α- and β-cristobalite, J. Phys.: Condens. Matter, 1995, 7, 1771,  DOI:10.1088/0953-8984/7/9/005.
199. M. Gambhir, V. Heine and M. T. Dove, A one-parameter model of a rigid-unit structure, Phase Trans., 1997, 61, 125–139,  DOI:10.1080/01411599708223733.
200. M. Gambhir, M. T. Dove and V. Heine, Rigid unit modes and dynamic disorder: SiO₂ cristobalite and quartz, Phys. Chem. Miner., 1999, 26, 484–495,  DOI:10.1007/s002690050211.
201. C. Ablitt, A. A. Mostofi, N. C. Bristowe and M. S. Senn, Control of Uniaxial Negative Thermal Expansion in Layered Perovskites by Tuning Layer Thickness, Front. Chem., 2018, 6, 455,  DOI:10.3389/fchem.2018.00455.
202. S. A. Mabud and A. M. Glazer, Lattice parameters and birefringence in PbTiO₃ single crystals, J. Appl. Crystallogr., 1979, 12, 49–53,  DOI:10.1107/S0021889879011754.
203. A. W. Sleight, Compounds That Contract on Heating, Inorg. Chem., 1998, 37, 2854–2860,  DOI:10.1021/ic980253h.
204. X. Xing, J. Deng, J. Chen and G. Liu, Novel thermal expansion of lead titanate, Rare Met., 2003, 22, 294–297,  DOI:10.1007/BF02672536.
205. J. Chen, K. Nittala, J. S. Forrester, J. L. Jones, J. Deng, R. Yu and X. Xing, The Role of Spontaneous Polarization in the Negative Thermal Expansion of Tetragonal PbTiO₃-Based Compounds, J. Am. Chem. Soc., 2011, 133, 11114–11117,  DOI:10.1021/ja2046292.
206. N. Shi, A. Sanson, Q. Gao, Q. Sun, Y. Ren, Q. Huang, D. O. de Souza, X. Xing and J. Chen, Strong Negative Thermal Expansion in a Low-Cost and Facile Oxide of Cu₂P₂O₇, J. Am. Chem. Soc., 2020, 142, 3088–3093,  DOI:10.1021/jacs.9b12442.
207. M. T. Dove, N. Shi, S. Wang, J. Song, J. Xu, W. Yin, H. Yang, J. Chen and G. Cai, Phase Transition and negative thermal expansion in copper pyrophosphate Cu₂P₂O₇ studied by neutron total scattering and the reverse Monte Carlo method, Phy. Rev. B: Condens. Matter Mater. Phys., 2026 , in press.
208. P. Lloveras, E. Stern-Taulats, M. Barrio, J. Tamarit, S. Crossley, W. Li, V. Pomjakushin, A. Planes, L. Mañosa, N. D. Mathur and X. Moya, Giant barocaloric effects at low pressure in ferrielectric ammonium sulphate, Nat. Commun., 2015, 6, 8801,  DOI:10.1038/ncomms9801.
209. G. Cai, J. Liu, H. Y. Playford, H. C. Walker, A. E. Phillips and M. T. Dove, Manuscript in preparation.
210. F. Qin, L. Hu, Y. Zhu, Y. Sakai, S. Kawaguchi, A. Machida, T. Watanuki, Y.-W. Fang, J. Sun, X. Ding and M. Azuma, Integrating abnormal thermal expansion and ultralow thermal conductivity into (Cd,Ni)₂Re₂O₇ via synergy of local structure distortion and soft acoustic phonons, Acta Mater., 2024, 264, 119544,  DOI:10.1016/j.actamat.2023.119544.
211. M. T. Dove, Comments on ‘Integrating abnormal thermal expansion and ultralow thermal conductivity into (Cd,Ni)₂Re₂O₇ via synergy of local structure distortion and soft acoustic phonons’, Materialia Scripta, 2026 , in press.
212. J. Guo, K. Refson and M. T. Dove, Submitted for publication.
213. Q. Gao, S. Zhao, S. Lorenzen, K. Zhao, Q. Sun, M. Finze, G. Cai, S. Kawaguchi, E. Liang and J. Chen, Biaxial zero thermal expansion in zinc tetracyanoborate, Phys. Rev. B: Condens. Matter Mater. Phys., 2026, 113, 064107,  DOI:10.1103/3sz4-rm67.
214. L. Wang, C. Wang, Y. Sun, S. Deng, K. Shi, H. Lu, P. Hu and X. Zhang, First-Principles study of Sc_1−xTi_xF₃ (x ≤ 0.375): negative thermal expansion, phase transition, and compressibility, J. Am. Ceram. Soc., 2015, 98, 2852–2857,  DOI:10.1111/jace.13676.
215. H. El Hyani, H. Habib, D. Tian and G. Li, Unveiling structural phase transitions and Negative Linear Compressibility of ScF₃ under high pressure, Phys. B: Condens. Matter, 2025, 701, 416976.
216. M. T. Dove, Structure and Dynamics: An Atomic View of Materials, Oxford University Press, 2003, p. 334.
217. M. T. Dove, Introduction to the theory of lattice dynamics, Neutrons et Simulations, 2011, 12, 123–159,  DOI:10.1051/sfn/201112007.
218. M. T. Dove, Theory of displacive phase transitions in minerals, Am. Mineral., 1997, 82, 213–244,  DOI:10.2138/am-1997-3-401.

---

Footnotes

† Note that in the mixed-cation fluorides, MM′F₆ discussed in Section 5.4.2, because of the doubling of the unit cell in each direction to give a face-centred Bravais lattice, the RUMs lie along the [100] directions, given that the point at the corner of the Brillouin zone in ScF₃ is at the centre of the Brillouin zone for the reciprocal space of the MM′F₆ lattice.

‡ The argument is more clearly articulated in the press release that accompanied publication of ref. 80.¹⁹³

---