Soft mode anisotropy in negative thermal expansion material ReO$_3$

We use a symmetry-motivated approach to analyse neutron pair distribution function data to investigate the mechanism of negative thermal expansion (NTE) in ReO$_3$. This analysis shows that the local structure of ReO$_3$ is dominated by an in-phase octahedral tilting mode and that the octahedral units are far less flexible to scissoring type deformations than the octahedra in the related compound ScF$_3$. These results support the idea that structural flexibility is an important factor in NTE materials, allowing the phonon modes that drive a volume contraction of the lattice to occupy a greater volume in reciprocal space. The lack of flexibility in ReO$_3$ restricts the NTE-driving phonons to a smaller region of reciprocal space, limiting the magnitude and temperature range of NTE. In addition, we investigate the thermal expansion properties of the material at high temperature and do not find the reported second NTE region. Finally, we show that the local fluctuations, even at elevated temperatures, respect the symmetry and order parameter direction of the observed $P4/mbm$ high pressure phase of ReO$_3$. The result indicates that the motions associated with rigid unit modes are highly anisotropic in these systems.


I. INTRODUCTION
The phenomenon of negative thermal expansion (NTE) is an intriguing and unusual property for a material to exhibit. Broadly speaking, there are two families of NTE materials: those in which the anomalous thermal expansion behaviour arises solely from vibrational effects, and those in which it has an electronic origin [1][2][3][4] . In a typical material, we expect to observe positive thermal expansion (PTE) since the anharmonic shape of the interatomic potential leads to an increase in the equilibrium distance between two atoms as they gain more energy from an increase in temperature. The typical explanation for how vibrational effects can lead to a deviation from this behaviour is the "tension effect" 1,4 . A lot of structural NTE materials consist of a network of cationanion linkages, often with an anion linked to two cations in a straight line [5][6][7][8] . If the energy cost to expand these bonds is quite high, then the central ion can displace perpendicularly to the bonds, which has the effect of pulling the two outer ions towards each other. The bonds connecting the linked ions still show PTE, but the linkage as a whole shrinks. This hypothesis is supported by neutron and X-ray powder diffaction experiments showing significant transverse atomic displacement parameters 4 . The tension effect is often realised in NTE materials via rigid unit modes (RUMs) 4,9 -materials made from a network of rigid polyhedra often have a high energy cost to distort the polyhedral units but a low energy cost for cooperative rotations of the polyhedra. These distortions typically cause a contraction of the volume of the material and since they are low in energy will have a large contribution to the overall thermal expansion behaviour. Many of these cation-anion linked NTE materials give rise to a network of connected polyhedra, such as the archetypal NTE material ZrW 2 O 8 10 , SiO 2 11 and layered materials showing uniaxial NTE such as Ca 3 Mn 2 O 7 12 . Two materials that are often used to illustrate the RUM model are ReO 3 and ScF 3 , due to their relatively simple structure when compared to more complex materials like ZrW 2 O 8 . Both compounds consist of cornersharing octahedra and both exhibit NTE, up to around 220 K in ReO 3 13,14 and 1100 K in ScF 3 15 , although the exact range of NTE in the former is dependent on sample preparation 16 . There has also been an observation of a reappearance of NTE in ReO 3 between 600-700 K 14 . This overlaps with the temperature at which ReO 3 is known to decompose to Re 2 O 7 and ReO 2 17 . Both ReO 3 and ScF 3 undergo a phase transition via octahedral tilts with applied pressure. In ReO 3 , the octahedral tilts in successive planes along the tilt axis are in-phase tilts 18-21 , whereas for ScF 3 , the tilts are out of phase 15,22 .
Compared to typical ABX 3 perovskites, which also exhibit rigid rotations of the octahedral units, ReO 3 and ScF 3 both have a vacant A-site, allowing larger tilt angles 19 . Recent analysis of the local structure of ScF 3 23,24 demonstrated that flexibility of the octahedra themselves, i.e., distortions of the intraoctrahedral F-Sc-F bond angles away from 90 • , is a key contributor to the NTE. In perovskites, pure RUMs are restricted to the line M-R in reciprocal space, which in principle gives their mode Grüneisen parameters a vanishingly small contribution towards the mean Grüneisen parameter. On moving away from this line, the modes have an increasing component of octahedral distortion. These modes, termed quasi-RUMs, can still contribute to NTE if they have negative Grüneisen parameters. Simplistically, structures with greater flexibility will have a greater volume of reciprocal space occupied by quasi-RUMs with a negative Grüneisen parameter than those of lesser flexibility. In addition, these modes will have a greater contribution to the mean Grüneisen parameter, since their component of octahedral deformation will have a lower energy cost 25 . The oxygen anions in ReO 3 have an increased charge compared to the fluoride anions in ScF 3 and therefore an increased coulomb repulsion force between them, which one might expect to decrease the flexibility of the octahedra. Molecular dynamics simulations of an ReO 3 -like structural model have also shown a decreasing magnitude of NTE for increasing anion interaction strengths 26 . In particular, the sign of the coefficient of thermal expansion in ScF 3 has been shown to be highly sensitive to changes in the force constant governing flexing of the F−Sc−F right angle 24 . We investigate this hypothesis by using a symmetrymotivated approach to analyze neutron pair distribution function (PDF) data collected on ReO 3 across its entire temperature range of stability. Our study enables us to identify the characters of the dominant dynamic deviations away from the average structure as a function of temperature.

II. EXPERIMENTAL DETAILS
Rhenium trioxide was purchased from Sigma-Aldrich and used as received. A 4.63965 g sample was loaded into a vanadium can of 6 mm diameter and mounted onto the Polaris instrument at ISIS Neutron and Muon Source (Rutherford Appleton Laboratory, U.K.) 27 . Data were collected in "short" 13 µA hr runs for Rietveld refinement and "long" sets of five 150 µA hr runs for pair distribution function analysis, equivalent to data collection periods of 5 minutes and 1 hour for each run, respectively. For low temperatures (between 4 K and 293 K) a helium flow cryostat was used, whereas for higher temperatures from 293 K to 750 K a furnace was used. The furnace data were collected in two separate experiments, due to an unscheduled beam shut-off at the facility during the first experiment. The furnace data for 600 K and above were taken from the second experiment. Both experiments used the same ReO 3 sample, stil loaded in the same sample can. During these experiments, data were also collected from the empty instrument, empty sample environment, and an empty vanadium can for Figure 1. (a) Temperature variation in the cubic lattice parameter of ReO3 from Rietveld refinement. The data points shown in red are from the longer collections which were used for total scattering measurements. (b) Plots of the PDF peak position ratios indicated in the legend. Peak positions for ScF3 were extracted from atomistic configurations generated using the reverse Monte Carlo method on neutron PDF data 24 and OriginLab was used to fit the peaks for ReO3. The error bars for ReO3 correspond to the uncertainties from the peak fitting algorithm. We note that values of (Re-Re)/2(Re-O) above 1 are not physically possible, but at worst they are less than 0.05 % too high. This is likely due to a systematic error and we believe that the trends with temperature are reliable. The uncertainty in these atom pair distances due to the r-space resolution (∆r ≈ 0.1 A) is an order of magnitude greater than this. background correction, as well as a solid vanadium rod for normalisation.
For average structure determination, data reduction was carried out using Mantid 28 . The lattice parameter for each temperature was determined by Rietveld refinement using the EXPGUI 29 The results of the SAPA on ReO3. The best individual fitting statistic is plotted for each irrep at each temperature. The Rw is shown relative to the Rw for the refinement with no symmetry adapted displacement modes active. The irreps are labelled as follows: colour denotes the k-point of the irrep, with blue referring to the M point, green to X, pink to R and yellow to Γ; marker shape denotes the irrep number, with a circle referring to 1, an upward-pointed triangle to 2, a star to 3, a square to 4 and a downward-pointed triangle to 5; linestyle denotes the parity of the irrep, with a solid line referring to a + irrep, and a dashed line referring to a − irrep. A representative fit of the average structure to the PDF at 293 K is shown in the SI. GSAS 30 against the data from Polaris detector banks 3-5, refining both the unit cell and atomic displacement parameters (noting that all atomic coordinates for the Re and O atoms are fixed by symmetry in the P m3m space group). A representative fit to the data at 150 K is included in the SI. The background was modelled using an 8-term shifted Chebyschev function. An absorption correction was refined to account for the neutron absorption of rhenium.
To ensure that the unit cell parameters were consistent between the cryostat and furnace environments, data were collected at 293 K in both environments. For the furnace data set, the unit cell parameter was fixed at the value determined from the cryostat data and the diffractometer constant DIFC for the backscattering detector bank was refined instead. This new value was held constant for subsequent furnace data sets to ensure selfconsistency across the entire temperature range (with an identical approach used to ensure consistency between the first and second furnace data sets).
For local structure analysis, Gudrun 31 was used to subtract the background, correct the data for self-shielding, absorption, and multiple scattering, normalise them to give the scattering function S(Q), and finally Fourier transform this to give the pair distribution function D(r). The scattering function S(Q) was determined over the Q range 0.6-50 A −1 in steps of 0.02 A −1 . This data normalisation is also dependent on the density of the powdered sample: the powder packing fraction was initially set to the value measured experimentally and then adjusted by hand for all data sets in order to set the limiting value of the total scattering structure factor, F (Q) 32 , as Q → 0 to its theoretical value of − i c ib 2 i 33 , and to set the coordination number for the first peak (Re-O) to its ideal value of six O atoms per Re atom.
Analysis of the pair distribution functions was carried out using the symmetry-adapted PDF analysis (SAPA) method described in ref. 23. For each sample, a 2 × 2 × 2 P 1 supercell of the P m3m aristotype ReO 3 with Re at (0.5, 0.5, 0.5) and O at (0.5, 0.5, 0) was generated and parameterised in terms of symmetry adapted displacements using the ISODISTORT software 34 . The modes modelled using this supercell expansion only represent a small fraction of possible phonon modes. However, even if the exact wave vectors of the soft phonon modes do not coincide with the Γ, X, M or R points of reciprocal space, since their characters should vary continuously between high symmetry points, the nature of the local symmetry breaking should still be manifested in our analysis. The generated mode listings were output in .cif format and then converted to the .inp format of the TOPAS Academic software v6 using the Jedit macros 35 . In total, there were 96 modes which transformed according to one of 19 irreducible representations. For each irreducible representation (irrep) at each temperature, refinements of the corresponding modes were started from random starting mode amplitudes. This was repeated 500 times for each irrep at each temperature to ensure the global minimum of the refinement was reached. For all temperatures, the refinements were carried out with a fitting range of 1.5 to 10 A. Refinements were carried out using the Topas Academic software v6 36 . The DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP) 37-40 , version 5.4.4. We employed the PBEsol exchange correlation potential 41 and projector augmented-wave (PAW) pseudopotentials 38,42 , as supplied within the VASP package. A plane wave basis set with a 900 eV energy cutoff and a 12 × 12 × 12 Monkhorst-Pack k-point mesh with respect to the parent cubic primitive cell (scaled accordingly for other supercells) were found suitable.

III. RESULTS AND DISCUSSION
Rietveld refinement of the a unit cell parameter from the Bragg reflection positions observed during the total scattering experiment shows we observe the same low temperature negative thermal expansion range as literature reports for ReO 3 13,14 ( Fig. 1 (a)). The high tem-perature measurements for both long and short runs do not show the second region of NTE observed by Chatterji et al 14 . Our conjecture is that the original observation of this phenomenon was likely due to sample decomposition, since above 673 K, ReO 3 starts to decompose via disproportionation 17 . We can gain some insights into the local structure from the PDFs of ReO 3 without performing any modelling. In  Fig 1 (b) To gain a fuller understanding of the character of the soft modes in ReO 3 , we turn to the results from our SAPA method. In Fig 2, we can see there is a group of irreducible representations whose modes consistently give the best improvement in R w compared to the average structure refinements at each temperature. These are M + 2 , M + 5 , M − 5 , X + 5 , X − 5 and Γ − 5 . All but M + 2 have symmetry-adapted displacements associated with them which are of a scissoring character: the Re-O octahedral bond lengths remain unchanged, but some of the bond angles are distorted via transverse displacements of the O anions. The M − 5 and X + 5 irreps both also support bond-stretching. However, this type of distortion is relatively high in energy so contributes negligibly both to the overall coefficient of thermal expansion and to the refined distortion. The amplitudes of these scissoring modes (Fig  3 (a)) are smaller than for ScF 3 23 which further supports the hypothesis that a lower flexibility is responsible for the reduced magnitude of NTE in ReO 3 compared to ScF 3 .
The other local symmetry breaking which is consistently amongst those which show the most improvement in R w from the average structure, and at some temperatures shows the most improvement overall, transform as irrep M + 2 . The distortions belonging to this irrep correspond to an in-phase tilt of the octahedra, with no other type of distortion associated with it, so it is of pure RUM character. This is a significant distortion for ReO 3 which undergoes a phase transition with applied pressure (ca. 5 kbar at 300 K) during which the M + 2 mode softens and the tilts are "frozen in" to the structure. To compare the relevant prevalence of this pure RUM mode with the scissoring modes, we refined a two-phase model, with displacements from the X − 5 irrep in one phase and from the M + 2 in another. The X − 5 irrep was chosen here since all the displacements associated with it are of a scissoring mode character and it provides a better fitting statistic than other irreps for which this is also true (M + 5 and Γ − 5 , Fig 2). At temperatures below 225 K, in the negative thermal expansion region, the RUM dominates as evidenced by the refined phase fraction in Fig 4 (a). At higher temperatures, there is a balance between the two phase fractions, indicating that both scissoring and RUM type motions account for a substantial proportion of the dynamic distortions.
At first glance, this increase in the proportion of motion arising from scissoring modes, coincident with the onset of PTE, is perhaps contradictory to the quasi-RUM mechanism for NTE which has been previously discussed for ScF 3 . In ScF 3 , the octahedral flexibility means a significant proportion of quasi-RUMs, modes which are of mixed scissoring and rigid unit mode character, will have a negative Grüneisen parameter. This expands the volume occupied in reciprocal space by NTE phonons, increasing the overall contribution of modes with negative mode Grüneisen parameters to the mean Grüneisen parameter. Previous analysis of both X-ray and neutron PDF data of ScF 3 finds that the majority of motion of fluorine atoms comes from scissoring modes, with approximate scissoring:RUM ratios of 4:1 (0.8 as a phase fraction) from SAPA analysis 23 and 3.5:1 (0.78) from reverse Monte Carlo analysis 24 at all temperatures. Our analysis in Fig 4 shows that in ReO 3 , the scissoring:rotation ratio is close to 1:1 at 675 K, but drops to 1:10 by 4 K. These ratios suggest that there is significantly more resistance to octahedral deformations in ReO 3 than ScF 3 . However, within renormalized phonon theory, we would expect the phase fractions 43 in Fig 4 to remain constant with temperature, as we observe in our RMC analysis of ScF 3 . The fact that the proportion of the scissoring to RUM phase fraction increases at high temperatures is more likely to reflect a hardening of the RUM than a sudden increase in the flexibility of the structure. This supposition is supported by the observation that the pressure at which the first phase transition occurs, increases with temperature 20 . The increased energy of the phonon modes with RUM character on warming lowers their contribution to the mean Grüneisen parameter at a given temperature. This hardening also has a knock-on effect on the quasi-RUMs, since their RUM component will also have a higher energy. It is of course tempting to suggest that this hardening of the modes with RUM character above 200 K is responsible for the concurrent change in sign of the bulk thermal expansion coefficient. However, it is conversely evident from the negative sign of the Grüneisen parameters of these modes that an in-  crease in volume would result in their hardening. Hence, we can not establish a causal connection with the current set of observations. Next, we consider the sequence of phase transitions in ReO 3 under pressure, since these are likely to be indicative of the character of the soft modes observed in the present study. The nature of the initial phase transition with pressure has come under question. The consensus in the literature is that at 5 kbar and 300 K, ReO 3 undergoes a transition to a tetragonal P 4/mbm structure involving an in-phase octahedral tilt along one pseudocubic axis, with only one arm of the propogation vector active (M + 2 (a; 0; 0)). At 5.3 kbar it undergoes a further transition to a cubic Im3 structure, with all 3 arms of the propogation vector active (M + 2 (a; a; a)) and a tilt of equal amplitude along all three pseudo-cubic axes. In Glazer notation, this corresponds to a transition from an a + b 0 b 0 tilt system to a + a + a + . Some experiments, however, report that there is no transition to the P 4/mbm structure 20 , only observing the transition to the Im3 phase.
Since we expect the distortion which is responsible for the phase transition to be a soft mode at ambient pressure, we can interrogate our PDFs to see if a precursor signature of this phase transition is already present. To do this, we parameterised the three dimensional M + 2 order parameter direction (OPD) (a; b; c) in terms of spherical polar coordinates with a = r cos φ sin θ, b = r sin φ sin θ and c = r cos θ. Refinements were performed at fixed values on a grid covering the range of values of θ and φ, whilst allowing the amplitude of the mode to vary. In this parameterisation, the OPD corresponding to P 4/mbm symmetry occurs if both θ and φ are integer multiples of π/2, and for all values of φ when θ = 0, π. The Im3 OPD would occur when φ = π/4, 3π/4, 5π/4 or 7π/4 and θ = arctan(± √ 2). In the Landau theory of phase transitions, the free energy expansion is written as a linear combination of sets of polynomials in the components of the order parameter. These polynomials must be invariant under all of the symmetry operations of the parent space group, P m3m in this case. Since the invariant polynomial truncated at the second (harmonic) order is of the form a 2 + b 2 + c 2 , which in the spherical coordinate parameterisation is equivalent to r 2 , the anisotropic R w distribution we observe over the spherical surface (Fig 5), is indicative of significant anharmonicity.
At low temperatures, the lowest R w refinements are clustered around points corresponding to the P 4/mbm OPD (Fig 5 (a)), while the worst fitting refinements are clustered around points corresponding to Im3 symmetry. Halfway in between points described by P 4/mbm and Im3 symmetry, corresponding to the OPD (a; a; 0) with I4/mmm symmetry, a small improvement in the quality of fit is observed compared to the OPD with Im3 symmetry. These areas correspond to distortions involving tilts about two orthogonal axes. This likely reflects a saddle point in the energy between the best (P 4/mbm) and worst (Im3) fitting OPDs. The anisotropy in the fitting statistics, that becomes more evident at higher temperatures and at large mode amplitudes, points towards significant anharmonicity, since, as discussed above, at the harmonic (quadratic) level, all OPDs must be equivalent with respect to the free energy expansion. Comparing OPDs against our PDF data for Im3 and P 4/mbm phases also shows a clear preference for P 4/mbm at all temperatures, with an R w of 7.28 % for P 4/mbm and 7.83 % for Im3 at 293 K, respectively ( Fig S3 in the SI). This supports the consensus that this distortion is the one first reached on the application of moderate pressure. However, it is surprising that even at high temperatures, such a pronounced anisotropic signature of this OPD is evident. We investigate the origin of this anisotropy below.
For small mode amplitudes (< 0.5 A) of the M + 2 (a; 0; 0) and (a; a; a) OPDs, such as we observe here, the distortions give rise to an almost identical volume strain in our DFT calculations (Fig 6 (a)). The two distortions are also found to be equally soft in energy over the range of amplitudes we expect to sample dynamically (Fig 6 (b&c)), hence these calculations do not explain the observed anisotropy with respect to the dynamic displacements transforming as M + 2 . Possible strain coupling is accounted for in the DFT calculations since the lattice parameters of each distorted structure were relaxed during the energy calculations. The only factor not accounted for by the DFT is that the M + 2 (a; a; a) OPD can couple to displacements that transform as the M + 1 irrep, but this would only serve to decrease the energy of the Im3 phase further, and so cannot explain why we observe the P 4/mbm distortion dynamically at ambient pressure. Since the ground state DFT calculations fail to account for the anisotropy in the OPDs, we must consider these effects to be either due to entropy or anharmonic couplings between phonons. We discount the latter proposal as the anisotropy persists even down to OPDs with P 4/mbm and Im3 symmetries for increasing 2 × 2 × 2 supercell-normalised mode amplitudes. These values are obtained using DFT calculations. Lattice parameters were relaxed, respecting the space group symmetry. The vertical line shows the mode amplitude for M + 2 at 293 K (Fig 3) The energy difference between the two is plotted (c). The lines shown in (b) are 4 th order polynomial fits to the data points shown. The line in (c) is the difference between the fits in (b). The line in (a) is shown as a guide to the eye. the lowest temperature, where we have shown that any dynamic deviations are small and are accounted for almost exclusively by the M + 2 RUM (Fig 4). Hence, in the following paragraph, we explore the prior suggestion that entropy might dictate the anisotropy of the dynamic fluctuations.
By virtue of our experimental observation of the anisotropy of these modes, we have already shown them to be anharmonic in nature. For the OPD M + 2 (a; 0; 0) with P 4/mbm symmetry, the structural fluctuation, taken in the static limit, has additional degrees of freedom that may be realised as dynamic tilts in the directions perpendicular to the spontaneously condensed tilt, corresponding to a line in the phonon dispersion curve between the R and M points in aristotypical perovskite symmetry. The Im3 structure has tilts along all three pseudo-cubic axes, and consequently will have no additional degrees of freedom of this manner. This leads to the P 4/mbm phase, or indeed an anharmonic distortion of this character, having the greater vibrational entropy of the two. This means that it is favoured over other OPDs for the smaller distortion amplitudes that are realised at lower pressures or during dynamic, anharmonic fluctuations of the system.
The above scenario implies that we should be sensitive in our PDF analysis to low amplitude, harmonic distortions of RUM character that are orthogonal to the dominant anharmonic one. Indeed, we can see a signature of this from the spherical polar coordinate plots of the M + 2 OPD at 600 K, in the form of "rings" around the (a; 0; 0) OPD (Fig 5 (b)), corresponding to smaller amplitude tilts about the axes orthogonal to the propagation vector of the anharmonic RUM. This is further supported by a fit to the 600 K PDF data with a two-phase model. Each phase contained a large amplitude (0.6 A) M + 2 (a; 0; 0) distortion, the "anharmonic part", and one smaller amplitude distortion (< 0.3 A) of M + 2 (0; b; c) or R − 5 (0; b; c) to mimic deviation in the data due to the line of mainly dispersionless harmonic RUMs running from k = [1/2 1/2 0] to [1/2 1/2 1/2]. All non-zero ratios of in-phase to out-of-phase tilts resulted in a slight improvement to the fit compared to an all in-phase model ( Fig S6 in the  SI).
The importance of entropy in determining the sequence of soft mode phase transitions in Ruddlesden-Popper perovskites has been highlighted by us recently 44 . In this instance, we discussed how the layering in this structure effectively affords phonon modes with an octahedral tilt character (RUMs about axes perpendicular to the layering axis) a greater vibrational entropy than those with an octahedral rotation character (about the layering axis). The sequence of temperature induced soft mode phase transitions we observed in these compounds were consistent with the idea that the entropic cost of ordering a tilt is higher than that of ordering a rotation. Given that the change in entropy associated with the ordering of a single mode is in principle vanishingly small, it is maybe surprising that the phase transition pathway is dictated in this manner. However, it is the associated renormalisation of phonon modes with RUM and quasi-RUM character, occupying a significant volume in reciprocal space, that provides the non-vanishing contribution to the Gibbs free energy, directing the soft mode transition pathway. Our present results tentatively go beyond these ideas in two respects. Firstly, ReO 3 is not a layered perovskite, so there is no distinction between octahedral tilts and rotations. However, it is clear that the same arguments about entropy and phonon renormalisation apply when considering the ordering of tilts about one axis compared to tilts about three. Secondly, our results imply that the anisotropic character of the phonon mode is essentially determined by the entropic cost of the anharmonic fluctuation itself, even without this fluctuation reaching the static limit required to initiate a soft mode phase transition.
In conclusion, a symmetry-motivated analysis of the pair distribution functions of ReO 3 has shown that the presence of a rigid unit mode allows this material to exhibit NTE, but a lack of flexibility of the structure limits the magnitude and extent of the NTE behaviour. The rigid unit mode has been shown to be anisotropic, displaying a clear preference for an (a; 0; 0) order parameter direction, even at elevated temperatures, which is consistent with the P 4/mbm space group the structure achieves after its phase transition with pressure. We tentatively suggest that the anisotropy we observe in the tilt direction of the M + 2 RUM is effectively determined by the entropic cost of the fluctuation itself.