Combined computational and experimental investigation of high temperature thermodynamics and structure of cubic ZrO2 and HfO2

Structure and thermodynamics of pure cubic ZrO2 and HfO2 were studied computationally and experimentally from their tetragonal to cubic transition temperatures (2311 and 2530 °C) to their melting points (2710 and 2800 °C). Computations were performed using automated ab initio molecular dynamics techniques. High temperature synchrotron X-ray diffraction on laser heated aerodynamically levitated samples provided experimental data on volume change during tetragonal-to-cubic phase transformation (0.55 ± 0.09% for ZrO2 and 0.87 ± 0.08% for HfO2), density and thermal expansion. Fusion enthalpies were measured using drop and catch calorimetry on laser heated levitated samples as 55 ± 7 kJ/mol for ZrO2 and 61 ± 10 kJ/mol for HfO2, compared with 54 ± 2 and 52 ± 2 kJ/mol from computation. Volumetric thermal expansion for cubic ZrO2 and HfO2 are similar and reach (4 ± 1)·10−5/K from experiment and (5 ± 1)·10−5/K from computation. An agreement with experiment renders confidence in values obtained exclusively from computation: namely heat capacity of cubic HfO2 and ZrO2, volume change on melting, and thermal expansion of the liquid to 3127 °C. Computed oxygen diffusion coefficients indicate that above 2400 °C pure ZrO2 is an excellent oxygen conductor, perhaps even better than YSZ.

Thus, we conclude that no direct experimental measurements of fusion enthalpy for ZrO 2 and HfO 2 have been performed to date.
In this work, we sought to fill this gap in the available data by measuring and computing the fusion enthalpies of ZrO 2 and HfO 2 . The combination of experimental and computational method offers a unique opportunity for corroboration that is essential given the challenges associated with each approaches. On the experimental side, the difficulties lie in thermal gradients unavoidable in conditions of uniaxial laser heating of aerodynamically levitated samples used for calorimetry 17,18 and X-ray diffraction [19][20][21] . On the computational side, the difficulties reside in reaching sufficiently large system sizes and sufficiently long simulation times while still using accurate electronic structure calculations as well as ensuring proper modeling of all forms of excited states (defect formation and diffusion, potential anharmonic phonons and electron excitations). Thermal expansions of cubic ZrO 2 and HfO 2 and volume change during the transition from the tetragonal phase were measured by high temperature X-ray diffraction experiments. The agreement between computed and measured values for fusion enthalpies and for thermal expansion supports the validity of the heat capacities, diffusion coefficients, and volume change upon melting obtained from the computation.

Results and Discussion
A summary of the results of ab initio computations is presented in Table 1 and Fig. 1. Results from high temperature X-ray diffraction are tabulated in Supplementary Information. Below, the thermodynamic data for cubic ZrO 2 and HfO 2 from computation and experiment are discussed together in the same order as in Tables 2 and 3 and are compared with literature values.

Tetragonal -cubic transition and thermal expansion of cubic phases.
Temperatures for tetragonalcubic transition and melting points for ZrO 2 and HfO 2 were accepted from the WZA 6 assessment and were used in this work for the evaluation of the temperature of the diffracted volume of the laser heated samples. Cubic ZrO 2 and HfO 2 have a fluorite structure with space group Fm3m and 4 formula units per cell (Z = 4). Besides the mineral fluorite (CaF 2 ), which gives the name for the structure type, natural and synthetic uraninite (UO 2 ), thorianite (ThO 2 ), and cerianite (CeO 2 ) are found in this structure. Thermophysical properties of UO 2 and ThO 2 above 2000 °C were studied extensively for nuclear reactors safety assessments 22,23 , and a comparison of the high temperature structures for UO 2 with ZrO 2 and HfO 2 from this work is given at the end of this paper. In the tetragonal (P4 2 /mmc, Z = 2) and cubic phases, Zr and Hf are coordinated by eight oxygen atoms, but in the monoclinic structure (P2 1 /c), stable at room temperature, the cation coordination is 7. Unit cell parameters of the tetragonal and cubic ZrO 2 and HfO 2 at transition temperatures were refined from XRD patterns containing both phases (Fig. 2), giving volume change upon transition. There are a number of values in the ICSD database 24 for volumes of stable and metastable tetragonal ZrO 2 at temperatures below 1627 °C (see Supplementary Information). Our value for the volume of tetragonal ZrO 2 at the transition temperature is consistent with the trend of close to linear volume expansion of the tetragonal phase, yielding an average value for of volumetric thermal expansion (α v ) of 3.9·10 −5 K −1 in the 300-2311 °C range.
At the transition temperatures, refined unit cells (a × c) for tetragonal ZrO 2 and HfO 2 are 3.690 × 5.337 Å and 3.669 × 5.327 Å, respectively. The variation in cell parameters from Pawley refinements of individual patterns on the same sample is within ±0.001 Å (Tables S2 and S3). Based on refinements from different beads, accuracy is estimated to be within ±0.003 Å, and this uncertainty has been propagated to experimental density and thermal expansion values in Tables 2 and 3. Both ZrO 2 and HfO 2 show less than 1% volume increase during the tetragonal to cubic transformation. From X-ray diffraction, unit cell parameters for cubic ZrO 2 increases from 5.265 to 5.291 Å from the tetragonal-cubic transformation temperature to the melting temperature. The corresponding values for HfO 2 are 5.246 and 5.265 Å. The cell parameters from ab initio MD computations (Table 1) show good agreement with experiment with differences less than 0.5%. After propagation of uncertainties, the experimental values for volumetric thermal expansion are in good agreement with computations and within 4 (±1)·10 −5 K −1 for both cubic ZrO 2 and HfO 2 in their stability range (Tables 2 and 3). We did not locate any previous reports on experimental or computational values for the thermal expansion. There are few reported values for the cell parameters of high temperature cubic ZrO 2 and HfO 2 and all of them were measured in vacuum and thus on possibly somewhat reduced samples. In fact, even though cubic ZrO 2 was assumed in the early phase diagrams by Kelley 16 in 1936 for the assessment of the fusion enthalpy, the existence of pure cubic phases at high temperature was still questioned in 1962 25 , due to the lack of structural data in oxidizing conditions. Boganov et al. 26 studied high temperature transformations in ZrO 2 and HfO 2 in a vacuum of 5·10 −6 Torr with heating by the electron beam and reported the unit cell parameter for ZrO 2 as 5.256(3) Å at 2330 °C and for HfO 2 as 5.30(1) Å at ~2700-2750 °C. The latter value for cubic HfO 2 was cited in reviews by Glushkova 27 and Wang 12 . Considering experimental conditions, these values probably refer to oxygen deficient cubic ZrO 2−x and HfO 2−x , known to exist in Zr-O and Hf-O systems 28 , and thus the differences with the results of our work are expected. Passerini 29 derived room temperature cell parameters for cubic ZrO 2 and HfO 2 as 5.065 Å and 5.115 Å by extrapolation from their fluorite solid solutions with CeO 2 . Combining his values with cell parameters before melting from this work (5.291 and 5.265 Å) gives an average volumetric thermal expansion from room temperature to the melting points of ~5·10 −5 K −1 for ZrO 2 and ~3·10 −5 K −1 for HfO 2.  30 . The density change of cubic ZrO 2 and HfO 2 with temperature from high temperature XRD data is shown in Tables 2 and 3 and Fig. 3 and compared with the results from computations. The good agreement allows us to rely on ab initio MD results for volume change upon melting as well as density and thermal expansion of the liquid phases. ZrO 2 and HfO 2 show similar expansion upon melting, 11 ± 2% and 10 ± 2%, respectively. For the temperature range sampled by computation, volumetric thermal expansions of liquid ZrO 2 and HfO 2 fall within (8 ± 1)·10 −5 -twice that for the cubic phase. Despite known biases in the lattice parameters calculated via DFT methods 31 , calculated volume changes tend to be much more accurate, due to systematic error cancellations. To the best of our knowledge, the volume change on melting has not been previously quantified. There are some published values for the density of liquid phases, since it has to be measured or refined for the analysis of the liquid structure by the pair distribution function (PDF) method 10,11 . The density of liquid ZrO 2 was recently measured by Kohara et al. 11 from the dimensions of aerodynamically levitated liquid ZrO 2 spheroids. His values, 5.1-4.9 g/cm 3 for in 2710-3000 °C, are in good agreement with our results 4.86-4.74 g/cm 3 at 2827-3127 °C. The density of liquid HfO 2 was refined from PDF measurements by Gallington et al. 10  Enthalpy and entropy of fusion. Ab initio MD computations resulted in values for fusion enthalpies of (ΔH fus ) 54 ± 2 kJ/mol for ZrO 2 and 52 ± 2 kJ/mol for HfO 2 . They agree, within experimental uncertainties, with values from the drop and catch calorimetry from samples levitated in argon flow (Fig. 4). It must be noted, however, that calorimetry experiments performed in oxygen flow did not provide a well defined step for HfO 2 fusion and resulted in a larger value for ZrO 2 (see Supplementary Information). This cannot be related to the sample reduction during levitation in Ar flow, as the calorimeter is not enclosed in the chamber, there is enough air entering in the levitation stream through turbulence to prevent ZrO 2 and HfO 2 reduction, and the samples were white in color after the drop experiments in Ar. We attribute observed differences to possible oxygen dissolution in ZrO 2 and HfO 2 melts, an effect previously observed by Coutures 32 in a number of oxide melts. The possibity of oxygen dissolution in molten ZrO 2 has profound implications for Zr-O phase equilibria at high oxygen fugacities and deserves a separate in-depth study.
In most of the assessments of the thermodynamic functions of cubic ZrO 2 , reviewed in detail by Wang et al. 28 , fusion enthalpy was kept fixed to the value 87 kJ/mol from JANAF tables 14 ; it was optimized to 68 kJ/mol in Chen et al. 's 33 assessment for the ZrO 2 -YO 1.5 system, while Chevalier et al. 34 obtained 90 kJ/mol. The value 87 kJ/mol was also accepted in the WZA 6 assessment and used for the calculation of ZrO 2 fusion entropy (ΔS fus = ΔH fus /T m (K) = 29 J/mol/K). The fusion enthalpy for HfO 2 was then estimated 6 from the melting temperature based on the assumption that it has the same fusion entropy as ZrO 2 . It must be noted that the widely used value for ΔH fus ZrO 2 from the JANAF Thermochemical tables 14 can be traced back to the 1951 data compilation by Wagman et al. 15 The same value is reported in Glushko's 35 compendium of thermodynamic properties of individual substances with reference to Kelley 1936 16 who calculated heats of fusions from available freezing point data of the binary systems ZrO 2 with SiO 2 and MgO. This approach is limited by a lack of information concerning activities in the liquid and solid solutions and by its reliance on the accuracy of the phase diagram determination. It is impressive that Kelley's 16 early assessment held for 80 years without any challenge.  The values computed for the fusion enthalpies of ZrO 2 and HfO 2 in this work (54 ± 2 and 52 ± 2 kJ/mol) are substantially lower than Kelley's 87 kJ/mol value included in JANAF tables 14 and used in the most current thermodynamic assessments 6,36 . They agree, within experimental uncertainties, with our drop and catch calorimetry measurements. Using accepted melting temperatures, the entropy of fusion for ZrO 2 and HfO 2 calculated as 17 J/ mol/K, which is substantially lower than ΔS fus 29 J/mol/K value obtained from Kelley's estimate and used in the WZA 6 assessment. While our experiments were in progress, a ΔH fus for ZrO 2 was reported by Kim et al. 36 as 26-49 kJ/mol from classical MD simulations based on interatomic potentials. We did not locate any reports on the computation of the fusion enthalpy of HfO 2 .
Heat capacities. As previously discussed 17 , drop and catch calorimetry cannot yet provide reasonably accurate values for the heat capacity due to differences in heat loss by radiation from different temperatures. The heat capacities of cubic ZrO 2 and HfO 2 obtained from ab initio MD computations are 111 ± 7 J/mol/K and 126 ± 4 J/ mol/K, respectively. The values computed for liquid ZrO 2 and HfO 2 at the modeled temperatures are close to the values for the cubic phase when uncertainties are taken into account (see Tables 2 and 3). Five different thermodynamic models have been proposed in recent years to model cubic zirconia and the liquid phase in assessments of the Zr-O system 28 . In the absence of data on the thermodynamics of cubic ZrO 2, they relied mostly on the reproduction of ZrO 2−x -Zr(O) and ZrO 2−x -liquid phase boundaries 37 . Heat capacities of cubic and liquid ZrO 2 calculated from different assessments are reviewed by Wang et al. 28 and for most models, they are in the range of 75-90 J/mol/K for cubic ZrO 2 and 80-100 J/mol/K for the liquid, below 3727 °C. Our computed values are close to 15R (where R is the gas constant) and substantially higher than those used in the assessments and higher than the 9R high temperature limit of Dulong and Petit for the contribution of lattice vibration. With the exception of UO 2 , which melts at 2874 °C, there are no experimental data for the heat capacity of fluorite type oxides melting in a comparable temperature range. Ronchi et al. 38,39 reported measurements of heat capacity for UO 2 from 1600 °C to 5000 °C using a custom-designed laser flash instrumentation. Their results indicate that UO 2 heat capacity exceeds 20R before melting, decreases to 15R after melting and decreases further to the 9 R limit only above 4000 °C. The excess heat capacity in UO 2 at high temperature is attributed to both electronic transitions and to disorder on the oxygen sublattice. The latter is also known as the Bredig 40 transition, which is common among fluorite halides and oxides above 0.8·T m . Clearly, the high temperature heat capacity needs further study.   oxygen defects in the cubic phases cannot be ruled out. The quality of the diffraction data did not allow refinement of oxygen occupancies due to the strong correlation with atomic displacement parameters (ADP). Isotropic ADPs were refined from selected XRD patterns as mean square displacement amplitude U iso (Å 2 ) and estimated from snapshots of MD trajectories (see Supplementary Information). U iso for Zr and Hf and for oxygen in HfO 2

Atomic diffusion in cubic and liquid ZrO 2 and HfO 2 . Diffusion rates for Zr, Hf and O in cubic phases
and in the liquid obtained from simulations for cubic and liquid phases are shown in Fig. 1. The proximity of diffusion rates of oxygen in cubic and liquid phases explains high heat capacity in fluorite phase and suggests that the notion of "oxygen sublattice melting" is an accurate description of the Bredig transition. Note that tetragonal -cubic transformation in ZrO 2 and HfO 2 was suggested to be a second order transition 6 and occurs shortly after exceeding 80% of the melting temperature threshold for the Bredig transition in fluorite structure. Diffusion coefficients were calculated from the MD trajectories, according to equation 42 where D is diffusion coefficient, t is time, r is atomic position and N is number of atoms. Temperature dependent diffusion coefficients are summarized in the Supplementary Information. Our computations show negligible Zr and Hf diffusion rates in stability range of cubic phases: within 0.4-0.7·10 −6 cm 2 /s for Zr and 0.1-0.3·10 −6 cm 2 /s for Hf. Oxygen diffusion coefficients above tetragonal-cubic transition temperatures are an order of magnitude higher than those for cations, which suggests significant oxygen diffusion. Notably, modeling cubic HfO 2 200 °C below its stability field does not show noticeable difference in Hf diffusion coefficient but show 10 fold decrease in oxygen diffusion (Table S5). Kilo et al. 42 reported MD computations of oxygen diffusion in YSZ with 8 and 24 mol % Y 2 O 3 from 400 to 1600 °C. Figure 5 show oxygen diffusion coefficients in ZrO 2 as a function of temperature and as a function of Y content using values extrapolated from Kilo's study. It is remarkable that oxygen diffusion in pure ZrO 2 is higher than in YSZ at any temperature and linear dependence on Y content is observed.
Comparison with UO 2 , ThO 2 , and fluorite-related bixbyite and pyrochlore structures. Thoria and urania both retain a fluorite structure from ambient temperature to their respective melting points (2874 and 3367 °C). Both oxides are believed to exhibit Bredig transitons above 0.8·T m . ThO 2 is the only known Th oxide and expected to be more similar to cubic HfO 2 and ZrO 2 than UO 2 which is known to exhibit electronic transitions and substantial hypo-and hyperstoichiometry ranges with a fraction of U going into trivalent or pentavalent states. For UO 2 at above 0.8·T m the linear thermal expansion increases 23 to 30·10 −6 K −1 , Oxygen U iso to 0.12 Å 2 22 , and melting is accompanied by 10 ± 1% volume increase 43 . The data for ThO 2 at above 0.8·T m are scarce, hence the thermal expansion at the melting point was extrapolated 44 to be 14·10 −6 K −1 , and Oxygen U iso follows the trend for UO 2 22 , but was not measured above 0.8·T m . Yb 2 O 3 and Lu 2 O 3 melt at 2435 and 2490 °C, respectively, and are stable in bixbyite or C-type structure (Ia3, Z = 16), which is often described as a derivative of a defected fluorite structure having ordered vacancies. Their linear thermal expansion was studied 21 both in argon and oxygen and was reported to not exceed 8.5·10 −6 K −1 with U iso values for Yb and Lu below 0.05 Å 2 up to the melting temperature and U iso values for O less than 0.07 Å 2 . Lanthanum zirconate (La 2 Zr 2 O 7 or LZ) is an example of compound stable up to the melting temperature in the pyrochlore (Fd3m, Z = 8) structure, which is often described as a defected fluorite structure with ordering of both cations and oxygen vacancies. Neutron diffraction in Ar atmosphere indicates that it does not display an anomalous thermal expansion or oxygen mobility indicative of a Bredig transition. The linear thermal expansion of LZ was reported as ~7·10 −6 K −1 from above 1650 °C to the melting temperature of 2300 °C, with U iso values for O and La not exceeding 0.07 Å 2 and that of Zr remaining below 0.03 Å 2 20 .

Conclusion
The performed computations and experiments fill gaps in the available thermodynamic data for pure ZrO 2 and HfO 2 at temperatures where the cubic fluorite phase is stable and into the liquid range, thereby facilitating future assessments. The experimental confirmation of thermal expansion and fusion enthalpies validate the accuracy of computational approaches and open the way for further computational studies of the high temperature thermodynamics of more complex systems. Our combined approaches are easy to generalize from HfO 2 and ZrO 2 to a broader range of systems. Indeed, we have applied the same combined experimental and computational methods to a wide range of oxides (e.g., Y 2 O 3 18 , La 2 Zr 2 O 7 41 , and several rare earth oxides). In addition, the computational method has been employed to study dozens of systems 45 , including oxides 46 , carbides, such as the Hf-Ta-C-N system 47 , and metals 48 .
The electronic temperature was accounted for by imposing a Fermi distribution of electrons on the energy level density of states. The electronic temperature was set consistently with the ionic temperature. We used automated k-meshes generation with a k-point density of 15 3 /Å −3 in the Brillouin zone. First-principles molecular dynamics (MD) techniques were utilized to simulate atomic movements and trajectories. The MD simulations were carried out under a constant number of atoms, pressure and temperature condition (NPT, isothermal-isobaric ensemble) with a time step around 2fs. The thermostat was conducted under the Nosé-Hoover chain formalism 53,54 . The barostat was realized by adjusting the volume every 80 steps according to average pressure. Although this did not formally generate an isobaric ensemble, this approach has been shown 55 to provide an effective way to change volume smoothly and to avoid the unphysically large oscillation caused by commonly used barostats. MD simulations were carried out with 90 Zr (or Hf) and 180 O atoms in a periodic cell. Employing periodic boundary conditions is a completely standard way to model extended condensed phases in these types of calculations. The cell size is as large as 16 Å to reduce the finite-size effect. The liquid phase was prepared by heating the solid up to 6000 K (about twice the melting temperature) for 0.5 picoseconds. The liquid is then cooled to the simulation temperature. MD simulations were performed for a sufficiently long time to achieve convergence. The length of MD trajectory varies from 14 to 62 picoseconds, depending on convergence, but generally, 30-50 picoseconds were sufficient. On average, computations took about 25,000 CPU hours per data point, which required around two weeks on 64 cores of a computer cluster. Theoretical X-ray diffraction calculations were carried out using the AFLOW package 56 . MD trajectory was sampled every 80 ionic steps, which formed a set of snapshots that were used to generate X-ray diffraction patterns averaged for the final analysis.
Experiments. X-ray diffraction (XRD) and calorimetry experiments were performed on polycrystalline ZrO 2 and HfO 2 beads, 2-3 mm in diameter, prepared by melting of powders purchased from Alfa Aesar (99.98% or higher metals purity) with a 400 W CO 2 laser. Samples were first melted into oblate spheroids in a copper hearth, followed by melting in an aerodynamic levitator, as described in detail elsewhere 17 . High temperature XRD experiments were performed with the aerodynamic levitator 57 at beamline 6-ID-D at the Advanced Photon Source (APS) at Argonne National Laboratory. Diffraction images were collected with a Perkin Elmer XRD1621 amorphous silicon detector in transmission through the upper part of laser heated beads freely rotating in oxygen flow through a levitator nozzle. The X-ray beam (λ = 0.12359(7) Å) was collimated to 0.5 mm wide, 0.2 mm tall rectangular shape. All images were recorded as a sum of 120 0.1 s exposures. The diffraction images at room temperature with the laser off were recorded first; then the sample was heated by a 400 W CO 2 laser in 50-100 °C increments as monitored with a Chino IR-CAS8CS pyrometer with 1 mm spot size set to 0.92 emissivity and 0.85 window transmission corrections. Image calibration, integration, and sequential Pawley and Rietveld refinements of XRD patterns were performed with the GSAS-II software 58 , backgrounds were fitted manually for each pattern and were not refined. NIST CeO 2 SRM674b powder standard was used to calibrate detector tilt and rotation angles, beam center position and sample to detector distance (1036.2 mm). Unit cell parameters for laser melted monoclinic ZrO 2, and HfO 2 were refined using conventional powder XRD with internal NIST Si640C standard and Bruker D8 instrument. In a sequential refinement of high temperature patterns, sample displacement was refined for each sample bead at room temperature from calibrated cell parameter and fixed for all refinements of high temperature patterns. In Rietveld refinements, sample absorption and oxygen occupancy were not refined to avoid correlation with atomic displacement parameters. Fusion enthalpies for ZrO 2 and HfO 2 were measured using drop and catch calorimetry. The technique and apparatus were described in detail elsewhere 17,18 . Schematic diagrams and photographs are provided in Supplementary Information together with data from all experiments.