Unraveling the ground-state structure of BaZrO3 by neutron scattering experiments and first-principle calculations

The all-inorganic perovskite barium zirconate, BaZrO3, is a widely used material in a range of different technological applications. However, fundamental questions surrounding the crystal structure of BaZrO3, especially in regard to its ground-state structure, remain. While diffraction techniques indicate a cubic structure all the way down to T = 0 K, several first-principles phonon calculation studies based on density functional theory indicate an imaginary (unstable) phonon mode due to the appearance of an antiferrodistortive transition associated with rigid rotations of ZrO6 octahedra. The first-principles calculations are highly sensitive to the choice of exchange-correlation functional and, using six well-established functional approximations, we show that a correct description about the ground-state structure of BaZrO3 requires the use of hybrid functionals. The ground-state structure of BaZrO3 is found to be cubic, which is corroborated by experimental results obtained from neutron powder diffraction, inelastic neutron scattering, and neutron Compton scattering experiments.


INTRODUCTION
Perovskite-type oxides, of the general chemical formula ABO 3 , where A and B denote different metal ions, constitute an extremely important class of materials, with properties such as electronic and/or ionic conductivity, 1 multiferroicity, 2 piezoelectricity, 3 magnetocalorimetry, 4 and luminescence. 5 These properties make them highly attractive for a range of technological applications. Various cation substitutions and/ or variation of the temperature are common means to alter the properties of the perovskites, but their underlying mechanisms are, in several cases, unclear. In particular, the closeness in energy between the various temperature-dependent phases represents a formidable challenge, both for experiments and for a theoretical description.
Barium zirconate (BaZrO 3 ) based ceramic materials have recently gained considerable attention for applications in protonic fuel cell applications 6−10 and hydrogen separation membranes. 11,12 Protons are introduced into the perovskite lattice by a process involving acceptor-doping to the Zr site, followed by hydration in a humid atmosphere at which the protons form covalent bonds with the oxygen ions of the perovskite lattice. 1 The basic mechanism for the proton mobility consists of rotational diffusion of the protonic defect (−OH) around the oxide ion and a proton transfer to a neighboring oxide ion. The proton transfer is generally believed to be the rate limiting process and is driven by a bending type motion whereby two oxygen ions come close to each other and thereby facilitate the transfer. 1,13−16 The flexibility and motion of the oxygen sublattice becomes important.
Fundamental questions surrounding the crystal structure of BaZrO 3 , especially in regard to its ground-state structure, remain. Also in this case, the motion of the oxygen sublattice becomes important. While diffraction techniques indicate a cubic structure all the way down to T = 0 K, several firstprinciples phonon calculation studies indicate an imaginary (unstable) phonon mode due to the appearance of an antiferrodistortive (AFD) transition associated with rigid rotations of the ZrO 6 oxygen octahedra corresponding to out-of-phase tilts of sequential octahedra. An accurate description of, in particular, the oxygen sublattice motion therefore becomes critical for both establishing the groundstate structure of BaZrO 3 and for modeling the proton mobility in acceptor-doped BaZrO 3 . Here, we will scrutinize the ground-state structure of BaZrO 3 , both experimentally and theoretically.
The high-symmetry structure of ABO 3 perovskites is simple cubic, with the O ions at the face centers, the A ion at the cube corners, and the B ion in the body center ( Figure 1a). The cubic structure is common at sufficiently high temperatures, while upon lowering of the temperature the ABO 3 perovskites generally display a number of structural phase transitions to phases of lower symmetry. The BaZrO 3 perovskite may be a very rare exception that stays cubic all the way down to T = 0 K. A very common transition is the out-of-phase AFD transition ( Figure 1b) which is manifested in theoretical spectra as the appearance of an imaginary (unstable) phonon mode at the R-point of the Brillouin zone of the cubic perovskite structure, and in the experimental spectra as, socalled, soft modes. The AFD instability is found to be well correlated with the Goldschmidt tolerance factor, t = (r A + r O )/√2(r B + r O ), where r A , r B , and r O are the ionic radii for the A, B, and O species, respectively. For t = 1, there is a perfect match to a cubic lattice in which rigid spheres with radii of the corresponding ions are perfectly stuck together. For t < 1, the A cation is small as compared with the cage of the surrounding oxygen ions and rigid rotations of the BO 6 octahedra commonly occur to optimize the A−O bond distances. For example, CaTiO 3 with t = 0.97 17 remains cubic down to about T = 1523 K where it makes an AFD transition to a tetragonal structure, 18 whereas SrTiO 3 , with a nearly optimal tolerance factor of t = 1.01, 17 remains cubic down to T = 110 K. 19 Barium zirconate, BaZrO 3 , also with a nearly optimal tolerance factor of t = 1.01 17 is a prominent candidate perovskite to stay cubic down to T = 0 K. Experimentally, both X-ray and neutron powder diffraction (XRPD and NPD, respectively) data, indeed suggest a cubic structure down to, at least, T = 2 K. 20 Theoretically, however, first-principles calculations based on density functional theory (DFT) show varying results, highly sensitive to the choice of exchangecorrelation (XC) functional in the calculations. Specifically, the use of the local density approximation (LDA) 20−23 and the PW91, 24 PBEsol 25 and Wu−Cohen 26 functionals predicts an instability at the R-point of the Brillouin zone, thus pointing toward an AFD transition, whereas the use of PBE, 27 RPBE, 27 and PBE0 28 functionals predicts a stable cubic structure at low temperature.
In cases where the instability, as observed in the DFT calculations, is weak, it has been suggested that quantum fluctuations may suppress the phase transition in BaZrO 3 and the material therefore stays cubic all the way down to T = 0 K. 20,21,24,25 Such quantum fluctuations have been shown to completely suppress the predicted ferroelectric (FE) transition in SrTiO 3 . 29 However, it has been argued that they are not sufficient in BaZrO 3 to suppress the transition. Instead, due to the closeness of the energies of the phases allowed from condensation of the R mode, a structural "glass state" may be formed upon cooling and the system would appear cubic on average. 23 A somewhat similar idea, an "inherent dynamical disorder", has been put forward to account for the apparent local deviation from a cubic structure identified by Raman spectroscopy of BaZrO 3 powder. 30 It is suggested that correlated tilts of ZrO 6 octahedra may occur on a local length-scale, shorter than what is probed by conventional diffraction methods where only the average cubic structure is observed. However, recent Raman studies on a single-crystal sample of BaZrO 3 show no direct evidence for such "nanodomains", and the spectra are instead explained by classical second-order Raman scattering. 31 In fact, as our knowledge about the real ground-state structure of BaZrO 3 has advanced, it has become clear that this seemingly simple material is a very challenging one. It is only through a more systematic, combined experimental and theoretical study that a clear mechanistic picture of the ground-state structure of BaZrO 3 will emerge.
Here, we present a systematic analysis of the temperature dependence of the structure and dynamics of BaZrO 3 , using a combination of high-resolution NPD, inelastic neutron scattering (INS), neutron Compton scattering (NCS), and first-principles DFT calculations. NPD gives information on the average structure but also the displacements of the nuclei from their equilibrium positions through the atomic displacement parameters (ADPs). The INS measurements are used to probe the vibrational properties of the powder sample for wave vectors close the R-point of the Brillouin zone, as a function of temperature, and NCS provides direct access to the nuclear momentum distributions. The quantal and thermal fluctuations of the ionic positions are therefore obtained as a function of temperature both in momentum space, through the NCS measurements, and in ordinary space, through the extraction of the ADPs from the NPD measurements. Six different wellestablished XC functional approximations are used in the DFT calculations to carefully investigate the accuracy of the theoretical predictions and, together with the experimental studies, establish the structure of BaZrO 3 .

EXPERIMENTAL TECHNIQUES
2.1. Neutron Powder Diffraction. The NPD experiment was performed at the high-resolution two-axis diffractomer D2B at the Institut Laue Langevin (ILL), France. Data were collected on 6.3 g of polycrystalline BaZrO 3 powder, placed in a cylindrical sample holder of Al with a diameter of 8 mm inside a standard cryostat. The wavelength was set to λ = 1.051 Å using the (557) reflection of a vertically focusing germanium monochromator. Diffractograms were collected at T = 5, 100, and 300 K. The data were analyzed using the standard Rietveld refinement method as implemented in the FullProf software. 32 2.2. Inelastic Neutron Scattering. The INS experiment was performed on the thermal neutron three-axis spectrometer IN8 at the ILL, France. The sample, the same as used in the NPD experiment, was placed in an open cylindrical sample holder of Nba purely coherent scatterer, chosen to minimize the incoherent contribution of the elastic line. The measurements were performed at a constant neutron final energy, E f = 14.68 meV (corresponding to a final neutron wave vector of k f = 2.662 Å −1 ), using as focusing monochromator the (111) reflection of doubly bent Si single crystals Chemistry of Materials pubs.acs.org/cm Article and as analyzer the (200) reflection of Cu single crystals. This setup was associated with an energy resolution of 0.57 meV at the elastic line, and a momentum transfer (Q) resolution of 0.03 Å −1 . A pyrolitic graphite filter was used to suppress the third-order harmonics from the Si monochromator crystals. While the IN8 spectrometer may be used to measure phonon dispersions of single-crystal samples, it is used here to specifically measure the phonon modes at the Q-values corresponding to the Rpoints of the Brillouin zone, since all the directions are merged in a powder sample. We thus measured spectra at constant Q, referred to as ω-scans, at the Q-point (1.5 1.5 1.5), corresponding to |Q| ≈ 2.9 Å −1 . Measurements of the other accessible R-points were inconclusive, as the spectrum at (0.5 0.5 0.5) was contaminated by spurious signals and the spectra at (1.5 0.5 0.5) and (1.5 1.5 0.5) contain overlapping contributions from other high-symmetry points of the Brillouin zone. Spectra were collected between T = 5 and T = 500 K using a cryofurnace.
2.3. Neutron Compton Scattering. The NCS experiment was performed at the ISIS Pulsed Neutron and Muon Source in the U.K. on the VESUVIO spectrometer. 33 VESUVIO is an inverted-geometry spectrometer where the final energy of scattered neutrons is fixed to 4.9 eV by a nuclear resonance of a Au foil 34 and the incident energy is obtained by time-of-flight measurements. The sample, 14.15 g of BaZrO 3 powder, was loaded into a Sn sachet of the cross-sectional area 3 cm × 3 cm. NCS spectra were measured at T = 15 and T = 300 K using the VESUVIO backscattering detector banks, in the scattering-angle range between 130°and 170°and sample−detector distances between 0.4 m and 0.6 m. The temperature was controlled using a closed circuit refrigerator.
The formalism of NCS has been recently reviewed in ref 35. Importantly, as the energy transfer ℏω and the momentum transfer Q are in the range 1−800 eV and 30−200 Å −1 , respectively, the scattering process can be described within the impulse approximation. 36 By imposing the conservation of the total kinetic energy and momentum in the neutron-plus-atom system, it follows that the dynamical structure factor, S(Q,ℏω), simplifies to a linear superposition of contributions from each isotope of mass M in the sample, S M (Q,ℏω), centered along the recoil lines ℏω = ℏ 2 Q 2 /2M. Each function S M (Q, ℏω) is defined from the distribution of momenta p of the scattering atom of mass M in the direction of the momentum transfer QQ where n(p) is the momentum distribution, y M = (M/ℏ 2 Q)[ℏω − ℏ 2 Q 2 /2M], and J(y M ;Q̂) is the so-called neutron Compton profile (NCP). 35 For a powder sample, such as the one studied here, the momentum distribution n(p) can be measured as a function of the magnitude of p only, and the directional information on Q̂is lost. Moreover, in the approximation of an isotropic and harmonic potential, the NCP can be written as a normalized Gaussian of the form where σ M is the standard deviation of the momentum distribution of mass M. From the measurement of σ M , it is possible to define the mean kinetic energy of the atom of mass M as ⟨E kin ⟩ = 3ℏ 2 σ M 2 /2M. One should notice that the powder-averaged NCP can still provide information on the anisotropy of the local potential, as in the case of deuterium and oxygen NCPs in heavy water, 37 because the powderaverage of an anisotropic NCP corresponds to a Gauss−Hermite expansion rather than a simple Gaussian function. 38 In general, the standard deviations of the momentum and spatial distributions are related to each other via the Heisenberg uncertainty principle. Therefore, a direct measurement of σ M allows an estimate of the spatial delocalization of an atom of the order of 1/2σ M .

FIRST-PRINCIPLES CALCULATIONS
The first-principles DFT calculations were performed using the projector augmented wave (PAW) method 39,40 as implemented in the VASP 41,42 software. Six different approximations to the XC functional were used. The LDA and the two constraint-based semilocal generalized gradient approximations (GGAs), PBE 43 and PBEsol, 44 where the latter is specifically designed for solids, were also used. The consistent-exchange vdW-DF-cx functional, 45−47 henceforth abbreviated as CX, is a version of the van der Waals density functional method. 47−49 It captures truly nonlocal correlations and balances exchange and correlation by use of the Lindhard screening logic. 48 A fraction of Fock exchange is included in two hybrid functionals. The vdW-DF-cx0p, 50 henceforth abbreviated as CX0p, is a hybrid extension of the CX functional 51 and includes both nonlocal correlations and nonlocal Fock exchange. The fraction of Fock exchange is set at 20%, following a coupling constant scaling analysis 52 of binding in sparse matter. 50 HSE 53,54 (sometimes called HSE06) is a range-separated hybrid extension of the PBE functional using 25% Fock exchange for the description of the short-range Coulomb interaction. The range separation is described by a screening parameter μ = 0.2 Å −1 in an error-function weighting erf(μr)/r of the Coulomb interaction. 53 Convergence of the R mode turned out to be very sensitive to the oxygen PAW potential and the energy cutoff. VASP comes with several different PAW potentials for oxygen. Both the regular and the hard PAW potentials treat 2s 2 2p 4 as valence states but with different core radii. The standard PAW potential for oxygen has core radii r s = 0.635, r p,d = 0.804, and a nominal energy cutoff of 400 eV, while the hard has r s,p,d = 0.582 and 700 eV, respectively. For the LDA, the regular PAW potential for oxygen was used with an energy cutoff of 900 eV. The semilocal functionals PBE and PBEsol were treated using hard PAW potentials and an energy cutoff of 1200 eV, while for the functionals including nonlocal Fockexchange and/or nonlocal correlation, CX, and the hybrids HSE and CX0p, hard PAW potentials and energy cutoffs of 1600 eV were required for full convergence. Convergence turned out to be less sensitive to the k-point sampling and a 6 × 6 × 6 k-point mesh was deemed sufficient for the hybrid functionals, while 8 × 8 × 8 was used for all nonhybrids.
The phonon spectra were calculated using the frozen phonon method with the default displacement of 0.01 Å in 2 × 2 × 2 supercells containing 40 atoms and postprocessed in phonopy. 55 A 10 × 10 × 10 k-point mesh was used for sampling of the Brillouin zone for all phonon calculations. In order to compare to the experimental INS spectra, the theoretical neutron scattering cross sections for a powder sample of BaZrO 3 was derived, using the in-house PowderTAS code. 56 The theoretical spectra were obtained by powderaveraging of the coherent and incoherent, elastic, and onephonon emission cross sections, using phonon eigenvectors and eigenvalues from DFT based phonon calculations. The PowderTAS code performs, for a series of Q-values, the powder-average of the neutron scattering cross sections calculated for a large number of randomly oriented Q-vectors, hence forming a S(Q, ℏω) scattering map.  20 The corresponding parameters from the Rietveld refinement are reported in Table 1. We note, however, that the refined occupancy factors correspond to a deviation from the BaZrO 3 stoichiometry, an excess of BaO or equivalently, a loss of ZrO 2 ( Table 1). The lack of superstructure and secondary phase peaks indicates that deviation from the BaZrO 3 stoichiometry is accommodated by randomly distributed defects, which is consistent with the pronounced Lorentzian component of the peak shape (cf. Y in Table 1), which generally indicates a reduced size of the scattering domains.
The anisotropic ADPs are shown as mean-square displacement tensors with terms reported in Table 1 as ⟨u αβ 2 ⟩, where α and β correspond to Cartesian direction indexes. We observe that both Ba and Zr show isotropic ADPs, where the one for Ba is somewhat larger in magnitude. The O displacements are highly anisotropic, corresponding to an oblate ellipsoid, and the corresponding ADPs are considerably larger. Note that the ADPs also contain contributions, if any, of small static distortions, which may lead to its value being slightly overestimated with respect to the real mean-square displacement of the atoms.
4.2. Structure and Vibrations. Table 2 shows the equilibrium lattice constants for the six different XC functional approximations, as determined from the DFT calculations. The lattice constants range from the well-known underestimation in LDA to the likewise well-known overestimation of PBE, with about ±1%, respectively. The CX and PBEsol functionals, as well as the two hybrid functionals HSE and CX0p, predict values between the extremes and agree well with the experimental value.
The stability of the cubic phase can be investigated by determining the phonon dispersion curves. The result using phonopy 55 is shown in Figure 3, with the inset showing a closeup around the R-point. The three functionals LDA, PBEsol, and CX all show an AFD R-mode instability while PBE and the two hybrid functionals predict a stable cubic structure. This unstable phonon mode, associated with the octahedral rotation motion, is hereafter denoted R 25 . The corresponding frequencies are given in Table 2. The results using the hybrid functionals give a frequency 2−3 times larger compared with PBE.
The stability of the cubic phase can be further investigated by mapping out the potential energy surface of the AFD mode. A potential energy surface can be generated by displacing the oxygen atoms along the R 25 phonon eigenmode while keeping the lattice spacing constant. The surface generated in this fashion can be well described by the expression 21,25,57   The stiffness parameter κ and the anharmonicity coefficient α are then determined through a fit to the computed data. The result is shown in Figure 4 for the different functionals, and the corresponding data for κ and α are given in The existence of minima at finite values for the displacement parameter u indicates a possible phase transition into a tetragonal phase. The PBEsol functional shows a small minimum of −0.44 meV, and CX, a very small minimum of −0.033 meV. It is likely that quantum fluctuations at T = 0 K suppress the tetragonal phase transition for these two functionals. Using LDA, a minimum of −1.81 meV is obtained for u = 0.1381 Å (θ = 3.8°). By dropping the cubic constraint, the full relaxation leads to the tetragonal I4/mcm phase with c/ a = 1.0031 and the energy −2.13 meV per formula unit relative to the cubic phase. With regard to the role of quantum fluctuations, previous reports are conflicting in that one report argues that quantum fluctuations are not sufficient to suppress the transition using LDA, 23 whereas other reports state that it is indeed likely that quantum fluctuations will suppress the transition also for LDA. 20,21 This implies that, most likely, if quantum fluctuations are included, all functionals, maybe with the exception of LDA, will predict a cubic structure down to T = 0 K. This is consistent with the experimental finding.
The stiffness parameter κ, which is related to the instability, depends rather strongly on the lattice spacing as can be seen in Figure 5. Consider first the nonhybrid functionals. By comparing these functionals at the experimental lattice The experimental values for the lattice constant and the R 25 phonon mode frequency are also indicated (from sections 4.1 and 4.3, respectively).    and PBE now predicts the largest instability. This shows that the main effect of the inclusion of gradient corrections and nonlocal correlation on the R mode stability is secondary and is mostly a consequence of the adjusted lattice constant. For the hybrid functionals, the situation is different. By including the nonlocal Fock exchange the stiffness increases substantially and becomes positive at the experimental lattice constant. The two hydrid functionals behave very similar to each other, and once again, the main difference between the functionals is the different equilibrium lattice spacing. Both functionals increase almost linearly with increasing lattice spacing and with similar magnitude. By reducing the lattice spacing (by applying a pressure) the stiffness parameter κ will also become negative for the hybrid functionals. Indeed, a structural phase transition from the cubic to the tetragonal phase has been observed for BaZrO 3 at about 17 GPa. 58 For the two hybrid functionals, the ratio between the harmonic and the anharmonic terms in eq 3 is quite small for relevant values of u indicating that for the R 25 mode anharmonic effects should be relatively small. For the PBE functional the above ratio is considerably larger, and a more pronounced effect from the anharmonicity should be present. To study the anharmonicity and hence the temperature dependence of the R 25 mode frequency, we make use of the quasi-harmonic approximation. 55 By taking the lattice expansion due to the zero-point energy (ZPE) fluctuations into account the equilibrium lattice constant is up-shifted by about 0.007 Å. This shift introduces a slight hardening of the R 25 mode frequencies. In Table 2 ZPE corrected lattice constants a 0 ZPE and corresponding R 25 mode frequencies ℏω R ZPE are given. The temperature dependence is then determined using the quasi-harmonic approximation. For the two hybrid functionals we obtain only a minor increase, about 10% at 500 K (see Figure 8). For the PBE functional a larger temperature dependence is obtained, consistent with its larger anharmonicity, and at T = 500 K the R 25 mode frequency has increased with about 30%. Figure 6a shows the analytically calculated (Q, ℏω) scattering map of a powder sample of BaZrO 3 , with the Q-values corresponding to the R-point of the Brillouin zone indicated by the vertical dashed lines. We observe that only the (0.5 0.5 0.5) and (1.5 1.5 1.5) R-points are sufficiently separated from other high-symmetry points and Bragg peaks in order for the R 25 phonon mode to be clearly distinguishable as the lowestenergy inelastic contribution, appearing as an inverted bellcurve with a minimum at about 6 meV. Unfortunately, as the experimental spectra at the (0.5 0.5 0.5) R-point is found to be contaminated by a spurious signal in the energy region of interest, the subsequent analysis concerns the (1.5 1.5 1.5) Rpoint only (marked with the vertical full line in Figure 6a).

Inelastic Neutron Scattering Measurements.
A slice of the scattering map at the (1.5 1.5 1.5) R-point is shown in Figure 6b, detailing the contributions of each crosssection to the calculated scattering intensity, as well as the sum of these contributions convoluted to the experimental resolution. We observe that the elastic intensity is solely due to the incoherent elastic cross-section, and the inelastic spectra to the coherent one-phonon emission cross-section. By comparison between Figure 6a and b, we distinguish four regions of the inelastic spectra. The lowest-energy band (A) is, intensity-wise, mostly due to the R 25 phonon mode. The scattering intensity progressively increases with energy transfer in region (B) as the relatively weakly intense acoustic phonon branches, from the nearby Bragg peaks (211) and (220), disperse up to a Q-value matching the (1.5 1.5 1.5) R-point. The most intense inelastic contributions (C), constitute the first optic band, appearing as a broad asymmetric band centered at about 11 meV due to resolution effects. A dip in Chemistry of Materials pubs.acs.org/cm Article intensity is observed at higher energy (D), in a region where phonon branches are mostly dispersive, in between the 11 and 20 meV nondispersive optic bands.
In Figure 6c, we show the experimental INS spectrum measured by ω-scans at the (1.5 1.5 1.5) R-point, as a function of temperature. It is featured by the incoherent elastic line at ℏω = 0 meV, a broad asymmetric band in the range 10−13 meV, and a series of lower-energy bands appearing as shoulders to the optic band, including a relatively well-defined band at about 6 meV. By comparing with the theoretical spectrum in Figure 6b, we identify the 6 meV band to the R 25 phonon mode (A), and the 10−13 meV band to the first optic band (C).
In order to extract the temperature dependence of the R 25 phonon mode, we performed a peak fit analysis of the spectra over the energy region 4.7−8.7 meV, using as a fitting model two Gaussian functions corresponding to the bands A and B. The optic band has been omitted to reduce the number of free parameters, as its asymmetry would require up to three Gaussian functions to model, and also to use a single fitting model, as the optic band is not consistently measured in all data sets. Furthermore, the B band is used to model the tail coming from the optic band, with its center fixed at ℏω = 8.2 meV and its fwhm fixed to 2.1 meV. Fitted parameters are reported in Figure 7. The integrated intensity of the phonon mode (Figure 7a) decreases exponentially while increasing the temperature, which is consistent with the expected decrease of the intensity due to the Debye−Waller factor. The phonon energy (Figure 7b) and its fwhm (Figure 7c) are found to be (within error) temperature independent, with ℏω = 5.88 (15) meV and fwhm = 1.7(3) meV. The large phonon fwhm, compared to the estimated energy resolution of about 0.7 meV at ℏω = 6 meV, is rationalized by broadening due to the powder-state of the sample, as it is not a unique phonon line but a bell-shape from the powder-average that is measured (see A in Figure 6a). Furthermore, the phonon line may be broadened due to the presence of defects associated with the sample nonstoichiometry that limits the size of the scattering domains, which would also contribute to the large fwhm. The overall effect of the presence of such defects appears to be limited to finite-size effects of the scattering domains (Y parameter in Table 1, possible broadening of the phonon line), and as such should have no significant impact on the stability of the R 25 phonon. In overall, provided that the 5.9 meV band is indeed the R 25 phonon, we find it stable down to T = 5 K. This would indicate that the presence of a ground-state of lower symmetry than cubic, based on the collapse of the R 25 phonon mode, is unlikely.
In Figure 8 we compare the calculated and measured temperature dependence of the R 25 phonon mode. The calculated temperature dependence is determined using the quasi-harmonic approximation 55 and the results for all three functionals that have positive values for the stiffness parameter κ; PBE, HSE, and CX0p, are presented. The two hybrid functionals give the best result for frequency at T = 0 K, but they also reproduce the weak temperature dependence better compared with PBE. In particular, CX0p gives a very accurate description of the R 25 phonon mode, both its magnitude and temperature dependence.

Neutron Compton Scattering Measurements.
The experimental data of the NCS experiment is summarized in Figure 9. Figure 9a shows an example of the mass-resolved Compton profiles in the time-of-flight spectra, where the oxygen peak is well separated from the peaks of the heavier elements Ba and Zr. Here, each peak position, centered along a recoil line, depends on the geometry of the instrument, thus on the source−sample and sample−detector distances, and the scattering angle. For this purpose we show the sum of the spectra from all backscattering detectors after the subtraction of empty-can and closed-circuit-refrigerator backgrounds. Figure 9b shows the sum over all backscattering spectra after the conversion from time-of-flight to the y-space of the oxygen atom, y O , both at T = 15 and T = 300 K. In this representation,  Chemistry of Materials pubs.acs.org/cm Article all oxygen peaks from different detectors share the same position around y O ≃ 0 Å −1 , i.e., along the oxygen recoil line, with no dependence upon the instrument geometry. The peak at y O ≈ − 40 Å −1 is the superposition of NCPs from Ba and Zr, whose masses are too large to be resolved in two separate contributions. This latter peak is treated here as a background, that would not appear in the y O -space if the impulse approximation was completely fulfilled, i.e., for infinite values of Q. The amplitudes of these two peaks are defined by the stoichiometry of the system and the values of the neutron scattering cross sections of each atom. In the following analysis, the ratio of these amplitudes (1.29) was therefore fixed. As a first qualitative result, we note that each NCP is broader at higher temperatures as a result of the Boltzmann population factor. Recently, it has been shown that the analysis of the difference of NCPs can be directly related to the widths σ M of an NCP at two different temperatures automatically getting rid off any multiple scattering contribution or other background signal that is not dependent on the sample. 38 In this procedure, the change in temperature of a Gaussian NCP defined by standard deviations σ and σ + δ can be expressed as a Taylor expansion in powers of the ratio δ/σ: This model has been fitted to the oxygen NCP and to the effective peak from heavier masses in the difference of spectra in Figure 9c. By doing so, the width of the momentum distribution at T = 15 K was found to be σ O = 9.0 ± 0.5 Å −1 . This value can be related to the width of the square of the wave function in the real space of 0.056 ± 0.003 Å. At T = 300 K the width of the momentum distribution has increased to σ O = 10.7 ± 0.7 Å −1 .
The momentum distribution is related to the mean kinetic energy for the oxygen atom according to ⟨E kin ⟩ = 3ℏ 2 σ O 2 /2M. ⟨E kin ⟩ increases from 32 ± 3 meV at T =15 K to 45 ± 6 meV at T =300 K. At T = 15 K the kinetic energy is mainly due to the ZPE experienced by the oxygen atom. Using the classical expression, ⟨E kin ⟩ = (3/2)k B T, we can associate the isotropic experimental NCS value with an effective temperature T* affecting the oxygen dynamics. At T = 15 K, we obtain T* = 2⟨E kin ⟩/3k B ≃ 250 K. The ZPE effect is substantial also for oxygen, and it will affect the dynamics at low temperatures. Within the harmonic approximation the mean square displacement projected onto the individual atom i and Cartesian direction α as a function of temperature T can be computed according to 55 (5) where N is the number of unit cells, M i is the mass of atom i, ℏω sq is the eigenfrequency for phonon mode s with wavevector q, and e sq i,α is the corresponding eigenvector component. Our results at T = 0 K are given in Table 3 and compared with the experimental data at T = 5 K. We find that the oxygen displacement is very anisotropic. u O 2 11 ⟨ ⟩ corresponds to the displacement of the oxygen atom along the Zr−O−Zr bond direction and is quite small corresponding to a high frequency motion and a stiff bond. u O 2 22 ⟨ ⟩ is double degenerate and describes the motion in the plane perpendicular to the Zr−O− Zr direction. It corresponds to the low-frequency R mode, the tilting type of low-frequency motion of the oxygen octahedra. The mean square displacements for the heavier atoms Ba and  ⟨ ⟩, ⟨u Ba 2 ⟩, and ⟨u Zr 2 ⟩ are somewhat larger, but by taking the uncertainty in the Rietveld refinement into account, we conclude that the agreement between theory and experiments is very good. The temperature dependence of ⟨u 2 ⟩ is also computed and compared with the experimental data in Figure  10. All three functionals, PBE, HSE, and CX0p, reproduce the experimental results well, both at low and high temperatures.
In analogy with ⟨u 2 ⟩, the mean square momentum can be expressed as . In Table 3 we present our computed data for the various atoms. In this case the relation between larger and smaller values for the oxygen atom becomes reversed with y O 2 11 ⟨ ⟩ being larger than y O 2 22 ⟨ ⟩, consistent with the Heisenberg uncertainty relation. A stiffer bond gives rise to a higher kinetic energy and smaller atomic displacement in real space.
The more robust observable in the NCS measurements is the mean square momentum averaged over the three directions α 11 22 σ ≡ ⟨ ⟩ + ⟨ ⟩ The computed data at T = 0 K and the measured data at T = 15 K for this average are given in Table 3, and the temperature dependence is given in Figure 11. We find very good agreement for the average mean square momentum between measurements and calculations ( Figure 11). The PBE gives a slightly better agreement, but generally all functionals are (within error) in agreement with the experimental results. In sum, the agreement between experiments and theory, both for the mean square displacement and mean square momentum, is very satisfactory and gives a consistent picture of the oxygen anisotropic vibrational motion, both at low and high temperatures.

SUMMARY AND CONCLUSIONS
We have performed a study of the temperature dependence of the structure and dynamics of BaZrO 3 , using a multitechnique approach combining neutron powder diffraction, inelastic neutron scattering, neutron Compton scattering, and firstprinciples density functional theory calculations. The density functional theory calculations are performed using six different well-established exchange-correlation functional approximations. The small energy differences between different bulk phases make the theoretical modeling a challenge. While diffraction techniques indicate a cubic structure all the way down to T = 0 K, several first-principles phonon calculation studies based on density functional theory indicate an imaginary (unstable) phonon mode due to the appearance of an antiferrodistortive (AFD) transition associated with rigid rotations of ZrO 6 octahedra. We show that only by using hybrid functionals, where some Fock exchange is included, the theoretical modeling leads to predictions fully consistent with experiments.
Specifically, by using the experimental value of the lattice constant, we show that the four functionals LDA, PBE, PBEsol, and CX show a very similar behavior for the AFD mode, i.e. an unstable AFD mode with similar magnitude. The inclusion of truly nonlocal correlations in CX has only a minor effect on the AFD mode. The well-known under-and overestimate of the lattice constant for LDA and PBE, respectively, implies that by using the theoretical lattice constants PBE predicts a stable AFD mode while LDA predicts a highly unstable AFD mode. However, by using the two hybrid functionals, HSE and CX0p, the theory predicts that the cubic phase is stabilized down to T = 0 K and the AFD instability is eliminated, both using the experimental and theoretical lattice constants.
The combined analyses of neutron diffraction and Compton scattering data provide complementary information on the atomic displacements in real and momentum space. The diffraction data show that the oxygen displacements are anisotropic, and the data are in very good agreement with the theoretical results, both at low and high temperature. The three different functionals PBE, HSE, and CX0p, which all predict a stable cubic structure, give similar results for the mean squared atomic displacements in real and momentum space and agree within the experimental uncertainties with the neutron Compton data, both at low and high temperatures. This illustrates that neutron Compton scattering is a powerful technique, not only for the momentum distribution for the very light atoms as hydrogen, but also for heavier atoms as oxygen in complex materials.
Further, the analysis of variable temperature inelastic neutron scattering data shows that the low frequency R 25 mode is quite temperature independent and stable down to at least T = 5 K. This is accurately predicted by the two hybrid functionals, HSE, which is based on the constraint-based, semilocal PBE and CX0p based on a current-conserving implementation of the vdW-DF nonlocal-correlation method. Both functionals predict a stable AFD mode but CX0p performs better for the value of the R 25 -mode frequency and the ZPE-corrected lattice spacing. The result from CX0p compares excellently with the corresponding experimental data. Furthermore, the R 25 -mode shows a weak temperature dependence, which is predicted by both hybrid functionals. As will be detailed elsewhere, the CX0p functional also performs better than HSE for predictions of thermal expansion and of extended X-ray absorption fine structure data for BaZrO 3 , while the two functionals have matching performances for the dielectric constant of BaZrO 3 . The CX0p has the exact same nonlocal correlation as CX and vdW-DF1. 47 The nonlocal exchange helps set charge transfer effects while the CX nonlocal correlations helps set restoring forces at deformations. 59 We are pleased to find that CX0p provides an accurate description of the structure and phonons in the BaZrO 3 , a nontrivial theoretical modeling task.
To conclude, BaZrO 3 is one of the very few perovskites that stays cubic down to T = 0 K. To accurately describe the structure and vibrations of BaZrO 3 , hybrid functionals should be used. Quantum fluctuations are present but they are per se not responsible for the absence of an AFD transition at low temperatures.
assistance during the NCS experiment. The calculations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC).