Anharmonicity and scissoring modes in the negative thermal expansion materials ScF$_3$ and CaZrF$_6$

We use a symmetry-motivated approach to analysing X-ray pair distribution functions to study the mechanism of negative thermal expansion in two ReO$_3$-like compounds; ScF$_3$ and CaZrF$_6$. Both average and local structure suggest that it is the flexibility of M-F-M linkages (M = Ca, Zr, Sc) due to dynamic rigid and semi-rigid"scissoring"modes that facilitates the observed NTE behaviour. The amplitudes of these dynamic distortions are greater for CaZrF$_6$ than for ScF$_3$, which corresponds well with the larger magnitude of the thermal expansion reported in the literature for the former. We show that this flexbility is enhanced in CaZrF$_6$ due to the rock-salt ordering mixing the characters of two of these scissoring modes. Additionally, we show that in ScF$_3$ anharmonic coupling between the modes responsible for the structural flexibility and the rigid unit modes contributes to the unusually high NTE behaviour in this material.


I. INTRODUCTION
Research into materials that contract upon heating, termed negative thermal expansion (NTE) materials, has been steadily increasing over the past 30 years. The significance of the phenomena was first underlined by Sleight et al. in 1996 1 by linking the large, isotropic NTE of ZrW 2 O 8 to the crystal structure of the material, opening up the field to synthesis of new compounds. Since then, this field has been expanded to a wider range of materials, including simple oxides (such as Cu 2 O 2 and ReO 3 3,4 ) and metal-organic frameworks 5,6 . The rigid unit mode (RUM) model is a common way to explain the origin of NTE 7 . Materials made from rigid polyhedra have a significant energy barrier to distortions of the polyhedra, but a low barrier to collective dynamics such as rotations. These modes are often low in energy so have a significant contribution to the coefficient of thermal expansion, and they can lead to NTE via the tension effect -if two linked bonds are straight, or nearly straight, and stretching the bonds would take a large amount of energy, a transverse displacement of the central atom would pull the two other atoms closer together, resulting in a local decrease in volume, the magnitude of which would increase when the temperature is raised 8 . ReO 3 , a material made from corner-sharing ReO 6 octahedra (hence can be thought of as an A-site deficient perovskite), is commonly used to illustrate this model due to the complexity of the motion in more typical NTE materials such as ZrW 2 O 8 . The octahedra in this material are expected to dynamically rotate in an out-of-phase manner with respect to their neighbouring units about their average positions, resulting in a contraction of the structure whilst the material remains, on average, cubic 3 . It is two compounds similar to ReO 3 that are studied herein; the isostructural ScF 3 , and the A-site deficient double perovskite CaZrF 6 . Metal trifluorides adopting the ReO 3 structure typically undergo a transition from the P m3m cubic structure to a rhombohedral phase (R3c) upon cooling, via long range ordering of the MF 6 octahedra (aaain Glazer notation). The dynamic motion of these tilts was expected to be the mechanism for NTE in ScF 3 9 supported by the fact that a phase transition to the rhombohedral tilt phase is observed under hydrostatic pressure of 0.7 GPa at ambient temperature 9,10 and in the related material CoZrF 6 , whose high temperature phase is isostructural to CaZrF 6 11 . NTE is observed at a range of temperatures above the phase transition, but below it, once the phonon mode associated with the RUM has been frozen in, strong PTE is observed. Previous studies of these materials have shown large displacements of the fluoride ions perpendicular to the M-F-M bonds (M = Sc, Ca, Zr) 12,13 , consistent with a polyhedral rocking mechanism for NTE. Other studies have challenged the RUM model, concluding that only certain bonds were rigid 14,15 , rather than entire polyhedra, and that bond bending could be a contributor to NTE 16,17 .
Several studies have been performed recently to try and ascertain the origin of NTE in these materials. X-ray pair distribution function (PDF) analysis of two materials in the cubic MZrF 6 (M = Ca, Ni) series has shown that differing degrees of flexibility in M-F linkages results in isostructural materials having very different thermal expansion properties 18 . Lattice dynamics calculations of ScF 3 performed by Li et al. 19 showed mostly soft lattice modes that distorted the ScF 6 octahedra, however a 3 x 3 x 3 grid of unit cells was chosen, which excludes the zone-boundary wavevectors which the RUMs are confined to. Molecular dynamics simulations on the general ReO 3 structure 20 , with variable interaction strengths, arXiv:1911.12695v2 [cond-mat.mtrl-sci] 3 Feb 2020 suggest a degree of flexibility in the octahedra enhances NTE. Another conclusion from these simulations was that a weaker anion-anion nearest neighbour interaction enhances NTE, which is supported experimentally by the greater magnitude of NTE in ScF 3 compared to ReO 3 . There is experimental evidence from Raman spectra and inelastic neutron scattering that the large NTE in these materials cannot be accurately predicted with the quasiharmonic approximation 15,19 , so subsequently lattice dynamics calculations have been done to elucidate the connection between NTE and phonon anharmoncity, since the relatively simple structure compared to other NTE materials allows for a more detailed analysis. These calculations show that cubic 21 and quartic 13,19 anharmonicity contribute significantly to the temperature dependence of the thermal expansion coefficient. Other simulations have shown that modes with quartic potential can have an enhanced NTE compared to a single-well potential 22 . ABO 3 perovskites exhibit a wide range of octahedral tilt phase transitions, as classified by Glazer 23 , yet do not generally display phonon driven NTE. However, we have recently demonstrated how, by using a symmetry motivated approach to analysing PDF data, we can gain extra information on disorder and dynamics 24 . Our study on BaTiO 3 showed that this method is very sensitive to soft phonon modes of RUM-like character. Here, we use this method to probe the character of the low-lying thermal excitations in the titled compounds, where the amplitudes of such vibrations are believed to be very large.

II. EXPERIMENTAL DETAILS AND DATA ANALYSIS
Scandium trifluoride was used as supplied by Strem Chemicals. Synchrotron radiation X-ray total scattering experiments were conducted at the synchrotron facility PETRA III (beamline P02.1 25 ) at DESY, Hamburg. A wavelength λ = 0.2070Å was used to collect data. Data were collected at temperatures of 125, 140, 147, 152 K and at intervals of 25 K from 175 to 450 K. The obtained 2D images were masked and radially integrated using the DAWN 26 software. G(r) and D(r) functions were computed using GudrunX 27 using a Q max = 21Å −1 . Gu-drunX was also used to perform background subtraction and sample absorption corrections.
The CaZrF 6 was that prepared via a standard solid state synthesis methods in Ref 18. The total scattering data was collected at 11-1D-C APS, Argonne National Laboratory, using a wavelength λ = 0.117 98Å between 25 and 400 K. The PDFs were computed using PDFGetX2 28 , which was also used for background subtraction and sample absorption corrections. A Q max = 28Å −1 was used for the analysis presented below.

A. Pair Distribution Function Analysis
Some form of modelling is usually required to extract information of interest, such as local distortions of atoms away from their high symmetry positions, from pair distribution functions. The method presented here involves expanding the possible degrees of freedom in terms of symmetry adapted displacements of the zone centre and zone boundary irreducible representations (irreps) of the P m3m A-site deficient perovskite structure. For this analysis we use a parent P m3m perovskite with the A-site at the origin. Symmetry-breaking displacements transforming as the same irrep can be further decomposed into symmetry adapted distortion modes by choosing a sensible basis that reflects the chemistry and crystallographic axes of the structure. The distortion modes have a 1:1 correspondence with phonon eigenvectors in the limit that only one set of atomic displacements transforms as the corresponding irrep. In cases where distortions from different Wyckoff sites transform as the same irrep, the character of the low lying excitations can still be ascertained through refining the relative amplitudes of the individual distortion modes. An overview of the displacements that enter into each irrep is tabulated in a recent paper by Popuri et al. 29 . For both compounds, ISODISTORT 30 was used to generate a model parameterised in terms of symmetry adapted displacements. A 2 x 2 x 2 P1 supercell was used for ScF 3 , since this allows phonon modes with propogation vectors k = [0 0 0], [1/2 0 0], [1/2 1/2 0] and [1/2 1/2 1/2] to be modelled. Whilst this is only a small fraction of possible wave vectors, these are both the ones that PDF data has the greatest sensitivity to and for which our symmetry motivated approach provides the greatest number of constraints. Furthermore, even if the exact wave vectors of the NTE driving phonons are of a longer wavelength, we still expect the character of those phonons to be reflected in our results which probe a shorter wavelength. To generate the parameterisation of CaZrF 6 , a 2 x 2 x 2 supercell of disordered Ca 0.5 Zr 0.5 F 3 was used. The cations were then set to be fully ordered to generate the rock-salt ordered structure. In all refinements, the breathing mode about the Ca/Zr (transforming as R − 2 ) was refined, making this description equivalent to the published F m3m structure 11 . The generated mode listings were output from ISODISTORT in .cif format and then converted to the .inp format of the TOPAS Academic software v6 31 . Modes transforming as the same irreducible representation (irrep) were tested simultaneously. An example of the best single-irrep refinement for each compound using this method is shown in Fig. 1. The results shown below ( Fig. 2) were performed with a fitting range of 1 (ScF 3 ) or 1.7 (CaZrF 6 ) to 10Å. The refinements were also done out to a higher radius, however the results were broadly similar for these larger fitting ranges. A comparison between the results for 10 and 30Å can be seen in the SI (Fig. S3).
The thermal parameters for each site were modelled with a simple quadratic, i.e b i = b i,low + ur + vr 2 where u and v are constant across all sites for each refinement, and b i,low is element dependent. Whilst this does not capture the true physical behaviour of the system, it was found to produce more robust fits to the data (more stable and fewer false minima) than other functional forms of b i , with the results still being consistent with our analysis performed using different functional forms of b i (see Fig. S2 in the SI). To get an unbiased view of how each irrep influences the local structure, the refinement for each irrep was initiated from randomised starting values of the relevant mode amplitudes. When a minimum was reached, the refined parameters were stored, re-randomised, and a new cycle was initiated. This process was repeated until 25000 iterations were reached (between 300 and 4000 refinements); this process was used to ensure a global minimum was reached for each mode. For refinements of atomic dis-placements transforming as the Γ 4 irrep, corresponding to ferroelectric type distortions, the amplitudes of modes affecting the metal cations were used to fix the origin; otherwise the mode amplitudes of this irrep would appear artificially high, due to the floating origin of the unit cell. Finally, we note that if the refined mode amplitudes are treated as the mean absolute value of displacement of an harmonic oscillator, then the amplitude of the harmonic motion will be a factor of √ 2 larger than the refined values.

B. Constrained Order Parameter Directions
Some order parameters can have many degrees of freedom associated with them. The exact number is a function of the degeneracy of the propogation vectors, the dimensionality of the irrep and the number of distortions transforming as the irrep. All of these degrees of freedom are described by the collection of symmetry adapted displacements or "distortion modes" that can be labelled accordingly. For example, in the parent structure (P m3m) of ScF 3 there are 3 types of distortion that transform as X + 5 which is two dimensional and associated with the triply degenrate k-vector [1/2 0 0], which results in a total of 18 parameters, compared to just 3 for M + 2 (a triply degenerate single dimensional k-vector) and R − 5 (a nondegenerate k-vector with 3 dimensions). In our refinements, to facilitate a fairer comparison between irreps, the order parameter direction (OPD) associated with the 3 wave vectors for each distortion have been set to the same values, i.e the general OPD (a,b;c,d;e,f) has been set to (a,b;a,b;a,b). Different distortion modes of the same type associated with the a and b branches of the OPD are allowed to have different values. However to further reduce the degrees of freedom that ratio between a and b across all distortion types that transform as a single irrep are fixed to be constant across different temperature ranges. This reduces the number of parameters for X + 5 from 18 to 4. Physically, these approximations correspond to a harmonic approximation in which the order parameter directions with respect to the propogation vectors and irrep dimensionality are strictly degenerate in energy. An example of this implementation is given in the SI.

III. RESULTS AND DISCUSSION
Rietveld refinement of ScF 3 and CaZrF 6 powder patterns can be used to gain some insight into the NTE behaviour, but can also be misleading; the average structure of both compounds remains cubic over the temperature ranges used, however, this structure fits the pair distribution function quite poorly, with PDFGUI 33 refinements of both structures from 1 − 10Å having R wp s ≈ 18 and 20% for ScF 3 and CaZrF 6 respectively (see SI Fig S1). The average linear coefficient of thermal is reported to have a magnitude of NTE approximately two to three times that of ScF 3 9,11 for the temperature range 25 -400 K, whereas in these measurements they have quite similar values. The difference to literature reports are in part due to the differing temperature ranges over which CTEs are reported, but may also be due to different strain, morphology and thermal histories of samples 35,36 . The refined atomic displacement parameters (Fig 2, top) reveal that most thermal motion of the F ions is perpendicular to the M-F-M linkages (M = Sc, Ca, Zr), indicating that a tensioning of these linkages could be responsible for the observed NTE. Some information can be gained from the PDFs without any modelling. Firstly, the effect of the rock-salt ordering of Ca 2+ and Zr 4+ in CaZrF 6 can be seen in the presence of 2 peaks at r ≈ 2Å, compared to just one in ScF 3 ; the greater positive charge of Zr 4+ compared to Ca 2+ means the Fions do not sit at the midpoint of Ca-F-Zr bonds (Fig. 1). Secondly, the relative magnitudes of the shorter inter-atomic separations (Sc-F, Sc-Sc; Ca-F, Zr-F and Ca-Zr) means that the magnitude of the mean M-F-M angle (M = Sc, Ca, Zr) must deviate from 180 • . The magnitude of this deviation is larger for CaZrF 6 than for ScF 3 (see Fig S3). The first peak for ScF 3 and the first two for CaZrF 6 are noticably less broad than the other peaks, indicating that the M-F bonds are relatively stiff. In contrast, the broadness of the F-F peaks at circa 3Å indicate a propensity for bending of the bond angles within the MF 6 octahedra. Little further information can be gained from a simple inspection of the PDFs, hence analysis of the structures has been performed in terms of symmetry-adapted displacements, as described in section II.A.
The results for the symmetry-adapted analysis are shown in Fig 2 (middle and bottom). The distortions can be classed into 3 general types: rigid unit modes, consisting of coherent rotations of the octahedra; semirigid "scissoring" modes, where there is a scissoring of some of the M-F bonds within the octahedra, and bondstretching modes, where some M-F bond lengths change. Most irreps in this analysis only have one distortion associated with them, although there are a few with more. There is a good degree of consistency between the two compounds; both have two "bands" of modes, one that fits well and one that fits poorly. The band with a ), all of which have distortions with a bondstretching character. The rest of the irreps, in the band that fits the data well, have at least one distortion associated with them that has a rigid unit (M + 2 and R − 5 ) or scissoring mode character. There are 4 zone boundary irreps that consistently have the lowest weighted R-factor for both compounds for the majority of temperatures -X + 5 , X − 5 , M + 5 and M − 5 . All of these irreps have one distortion associated with them that is of scissoring mode character. A depiction of the effect of these modes on the structure of CaZrF 6 is shown in Fig 3. The Γ − 5 irrep also fits well, especially in the refinements that go out to 30Å. The displacements associated with this irrep are also of a scissoring mode character, however despite the low R wp , the mode amplitudes are consistently small, hence most of the analysis is focused on X + 5 , X − 5 , M + 5 and M − 5 . The weighted mean amplitude over all refinements at each temperature for these irreps have been calculated and are shown in Fig. 2 (bottom), the weighting being given by a Boltzmann distribution w = exp[(R global − R wp )/σ], where R global is the minimum weighted R-factor achieved across all refinements and all temperatures for the relevant compound and σ is the value of a meaningful difference in the weighted R-factor, taken to be 0.1%. R global is taken to be 9% for both compounds. The amplitude of these modes (X + 5 , X − 5 , M + 5 and M − 5 ) are consistently higher for CaZrF 6 than for ScF 3 -this coincides well with the more significant distortion away from the av- erage structure for CaZrF 6 , as seen in the mean M-F-M bond angles, and the greater magnitude of NTE reported in the literature. These modes also fit significantly better than the RUMs (M + 2 and R − 5 ). These best-fitting irreps (X + 5 , X − 5 , M + 5 and M − 5 ) are all two-dimensional and all have 3 k-vectors, therefore the OPDs have been constrained as described in section II.B to allow for a fairer comparison with the RUMs, which have fewer degrees of freedom associated with them. For the ScF 3 , the unconstrained R − 5 (which is associated with the out-of-phase octahdral tilts observed in other metal trifluorides) has a similar quality of fit to the constrained X + 5 , X − 5 and M + 5 at lower temperatures and consistently performs better than M − 5 (Fig 4). This suggests that a combination of both the rigid unit and scissoring modes are responsible for NTE, which agrees with previous molecular dynamics studies of these materials 20 . In this study the authors argue that correlated dynamics of flexible polyhedra result in a greater degree of NTE than purely rigid unit dynamics. However for CaZrF 6 , we find the constrained scissoring modes, with the exception of M − 5 , consistently perform better than the RUMs. The RUMs also start to perform increasingly poorly as the temperature is raised above 100 K, suggesting that the thermal expansion in CaZrF 6 at higher temperatures may well be dominated by contributions from these scissoring modes. The increasing R wp of the RUMs as temperature is increased, and the contrasting decrease in R wp seen for the scissoring modes, tallies well with the phonon dispersion curves of both compounds 13,19 . These show that the scissoring modes are slightly higher in energy than the RUMs, so the scissoring modes will become more active at higher temperatures.
As discussed earlier, the different charges on the two cations in CaZrF 6 result in a need to refine the octahedral breathing mode, transforming as the R − 2 irrep, alongside the other distortion modes in order to facilitate a more direct comparison to ScF 3 . In the average structure of CaZrF 6 , this breathing mode is frozen in, lowering the symmetry from P m3m to F m3m. This also has the ef-fect of mixing the characters of some of the irreps such that the associated atomic displacements now transform as the same irrep. For example, the X + 5 and M − 5 irreps of P m3m correspond to the X − 5 irrep of F m3m and X − 5 and M + 5 correspond to X + 5 . To determine whether this mixing of characters has any effect on the observed local structure of CaZrF 6 , the constrained OPD X + 5 and M − 5 modes were refined together (hereafter referred to as X + 5 ⊕M − 5 ). This gave a significant improvement to the quality of the fit (Figs 4 and 5). To determine whether this coupling is a significant effect, results are compared to a two-phase model, in which modes transforming as different irreps are refined in separate phases. Hereafter these two models will be referred to as the "coupled" model (denoted with ⊕) and the "2-phase" model (denoted with &). The coupled modes have a significantly better R-factor above 100 K, but fit worse than the 2-phase refinement below this temperature. The same comparisons are also done for ScF 3 , where any coupling between phonons of these characters should only arise from anharmonic interactions. In contrast to CaZrF 6 , which shows a clear preference for coupling between X + 5 and M − 5 , no evidence of such coupling and hence an anharmonic interaction is seen here for ScF 3 . This suggests that whilst these scissoring modes are important in determining the local structure of ScF 3 , any anharmonic coupling between them has little influence on the lattice dynamics that drive NTE. A similar comparison is made for X − 5 /M + 5 , however the 2-phase refinements consistently fit better than the coupled model for both compounds. This may be due to both distortions locally having the same character (Eu of point group m3m) with respect to the MF 6 octahedra, making coupling unfavourable.
Next we investigate if the similar quality of fits of the scissoring modes X + 5 , X − 5 and M + 5 and the rigid unit mode R − 5 could be indicative that the two types of distortion are cooperatively coupled to produce the observed NTE. To test this hypothesis, we explore 2 scenarios; whether this observation is simply due to the dynamic distortions occurring in different sample volumes or at different times to each other or a coupled model which implies that significant anharmonic coupling between these modes is occurring. For both materials, the X + 5 /R − 5 refinements show a similar sort of behaviour to the X + 5 ⊕M − 5 refinements in CaZrF 6 , in that the refinements of the coupled modes perform worse than the 2-phase refinements at lower temperatures, but soon crossover to show an improved fit, although the results for CaZrF 6 are not robust. Since by the symmetry lowering of the rock-salt ordering in CaZrF 6 , X + 5 and M − 5 are allowed to couple, and we have shown our analysis to be sensitive to this coupling, the results in Fig. 7 are indicative that there is coupling between the X + 5 and R − 5 modes. However, as by symmetry coupling in X + 5 ⊕R − 5 is not permitted on its own, we construct a coupled distortion that forms an invariant in the free energy expansion by inclusion of the M − 5 irrep. The X + 5 and M − 5 OPDs in this refinement are still restricted, resulting in 3 more parameters than the X + 5 ⊕M − 5 refinements, but much improved fits (Figs 4 and  7). This model results in a very good agreement with the data (Fig. 7).
A very recent analysis of ScF 3 neutron PDF data using the reverse Monte Carlo (RMC) method 37 similarly concludes that it is a combination between structural flexibility and RUMs that causes the NTE in the compound. Dove et al. argue that the flexibility of the structure allows RUMs and RUM-like modes to occupy a larger volume in reciprocal space, meaning they give a greater contribution to the overall thermal expansion behaviour, compared to entirely rigid structures. Our results here echo this conclusion and underline the dominant contribution of scissoring modes in describing the fluctuations from the average symmetry. Additionally in the work of Dove et al., geometric algebra is used to quantify the proportion of the motion of the atoms in ScF 3 originating from correlated whole-body octahedral motion, deformations of the F-Sc-F right angles and changes to the Sc-F bond length. This analysis resulted in a ratio of approximately 7:2:1 of bends:rotations:stretches. The X + 5 &R − 5 and X − 5 &R − 5 2-phase refinements desribed previously give a similar ratio of bends:rotations, approximately 8:2, although the contribution from stretches is negligible (< 1% of the total motion). The X + 5 &R − 5 refinements for CaZrF 6 give an approximately 7:3 ratio of bends:rotations, again with a negligible contribution from stretches. There is hence a high degree of consistency between results derived via big box RMC methods and those of our symmetry motivated approach here. A different analysis of neutron PDF data of ScF 3 , performed by Wendt et al. 38 , models the F atoms as being randomly positioned on a torus-shaped gaussian distribution around the F sites in the average structure, with no correlation between neighbouring F atoms, in a similar fashion to entropic elasticity in polymers. The model reproduces the observed NTE behaviour and F-F distribution up to ≈ 700 K. It shows how important flexibility of Sc-F-Sc linkages are in this material, a fact consistent with our findings here, however it fails to account for the full range of NTE in the material. The previously discussed RMC model shows that at least a small fraction of the motion of F atoms in the material can be accounted for with correlated rigid-unit type distortions, results which are compatible with our symmetry based analysis of the X-ray PDF data.
In summary, we have shown via a symmetry motivated real-space analysis of PDF data that the most significant distortions in these ReO 3 -like NTE materials are scissoring modes, which involve scissoring of the MF 6 octahedral bond angles. These modes have a greater amplitude in CaZrF 6 than ScF 3 , which corresponds well to the greater magnitude of NTE reported in the literature for the former. Coupling between these modes and the rigid unit modes has been shown to be active and the likely origin of unusually high NTE in these structures.