BaZrO3 stability under pressure: The role of nonlocal exchange and correlation

The ground-stale structure of BaZrO3 is experimentally known to be cubic down to absolute zero. However, there exist several measured properties and experimental characterizations that earlier computational works have failed to accurately describe and explain within this cubic symmetry. Among these properties and observations are the dielectric constant and the parallel mean-squared relative displacement value that tracks the fluctuations in distance for Ba-O atom pairs. Previous density-functional theory (DFT) studies have resolved the issue by assuming that BaZrO(3 )undergoes a phase transition from cubic to tetragonal I4/mcm symmetry, possibly while forming a glasslike state that reflects cubic symmetry on average. In this paper, we show that the set of experimental results can indeed be satisfactorily explained by DFT entirely within the cubic symmetry. We find that past theory limitations arose from the choice of exchange-correlation-functional approximations and that the inclusion of Fock exchange in hybrids significantly improves the DFT performance. We also find that the inclusion of nonlocal correlation effects is beneficial. We conclude by making a prediction for the phase-transition pressure for the transition from cubic to tetragonal symmetry at zero kelvin.


I. INTRODUCTION
Perovskites comprise a large family of materials with the general formula ABO 3 , few limitations on A and B, and hence a large variety in crystal structure. In its idealized cubic Pm3m symmetry, the perovskite structure is built up of two interpenetrating simple cubic lattices of A and B cations with oxygen anions on the faces of the cube forming an inscribed octahedron. When the temperature is lowered, most perovskites are prone to some kind of distortion away from the high-symmetry cubic structure. Often, this distortion can be described either as a displacement of the B cations causing a ferroelectric transition or a (rigid) rotation of the oxygen octahedron. The former is associated with a softening of a phonon mode at the point, whereas the latter is associated with a softening of one of the zone boundary points such as the R or M point [1]. From the family of perovskites, BaZrO 3 stands out as a material which experimentally does not show any indication to a phase transition away from the cubic Pm3m, at least as far down as 2 K [2][3][4][5][6].
The symmetry-reducing distortions are often attributed to the ionic-size mismatch, described by the Goldschmidt tolerance factor [7]. This factor is sometimes used as a rule of thumb when determining the likelihood of a phase transition. A tolerance factor larger than unity indicates a cation * hyldgaar@chalmers.se displacement while a tolerance factor smaller than unity suggests a rotation of the inscribed oxygen octahedron. With a tolerance factor of 1.004, it is likely that BaZrO 3 remains cubic. However, a tolerance factor near (or slightly above) unity is only a necessary condition for cubic-phase stability. It is not sufficient to exclude the possibility of a phase transition.
Indeed, experimental data has been interpreted as reflecting inconsistencies with cubic symmetry. Several theoretical studies, based on the local density approximation (LDA) to density-functional theory (DFT), report antiferrodistortive (AFD) instabilities in BaZrO 3 [6,[8][9][10]. They find other symmetries with lower ground-state energies than the cubic in LDA [9,10]. A lower symmetry has also been assumed to explain a so-called anomalous Debye-Waller factor. This term refers to extended x-ray absorption fine structure (EX-AFS) measurements of a large mean-square relative deviation (MSRD) for the nearest-neighbor distance off the Ba atoms [9]. A lower symmetry structure has furthermore been proposed as the solution to the discrepancy between previously calculated and experimental values for the dielectric constant [10]. The discrepancies between previous LDA-based firstprinciples-theory calculations and experiments has not been addressed in detail.
In our previous work [11], we concluded that BaZrO 3 is cubic down to zero kelvin by combining neutron measurements and DFT. The calculations ranged from LDA to hybrids based on both the generalized gradient approximation (GGA) and the truly nonlocal correlation formulations of the van der Waals density functional (vdW-DF) method [12][13][14]. We documented experimentally that the R 25 mode, sometimes discussed as a soft mode, exhibits no softening with temperature. This is also what is predicted by the two hybrid functionals that we used, namely, HSE [15], and vdW-DF-cx0p [16]. The latter is henceforth abbreviated CX0p. It is a parameter-free hybrid formulation of the consistent-exchange vdW-DF-cx version [17,18], here abbreviated CX.
In that previous study [11], we also found that the PBE version of the GGA predicts the cubic phase down to T = 0, although the R 25 mode is significantly softer in PBE. By taking the effect of vibrational zero-point energy (ZPE) into account, the raw DFT lattice constant a 0 was seen to increase by a ZPE = a ZPE − a 0 = 0.007 Å in PBE and CX0p. For HSE, the shift is a ZPE = 0.008 Å [19]. Applying a similar shift, a ZPE = 0.008 Å, stabilizes the cubic phase in a direct determination of the R-mode frequencies described even in CX. Hence PBE, CX, HSE, and CX0p can, in principle, all be seen as predicting cubic-phase stability. However, this finding also raises a number of questions: (1) Why was there no stability found in previous LDA computational studies? (2) Is there one of these semilocal or nonlocal functionals that we can trust to accurately describe general BaZrO 3 properties? (3) Are there some BaZrO 3 measurements that, nevertheless, suggest that BaZrO 3 must be modeled as a noncubic materials?
In this paper, we show that there is no need to go beyond the assumption of cubic symmetry for an accurate description of a wide range of measurable BaZrO 3 properties. Thepoint frequencies, the dielectric constant, and the EXAFS measurements of the parallel MSRD values (tracking fluctuations in the Ba-O distances), can be accurately described computationally within the cubic symmetry. By the inclusion of Fock exchange (in HSE and CX0p), the calculated values are in very good agreement with the experimental values.
Also, we explore the role that an increasingly nonlocal exchange and an increasingly nonlocal correlation description can play in providing accurate materials characterizations. Multiple types of BaZrO 3 measurements motivate theoretical characterizations of these properties for the set of PBE/CX/HSE/CX0p functionals. We note that LDA does not give an accurate prediction of BaZrO 3 stability. Changing from LDA to first GGA and then to HSE or CX0p hybrids means an increasingly nonlocal exchange description that, generally, improve the description of charge transfers [20]. Meanwhile, changing from GGA to CX or from HSE to CX0p provides an increasingly nonlocal-correlation description. This change can help in setting the restoring forces [14,[21][22][23]. Our cross-functional characterization allows us to track how these factors affect the theory accuracy.
Finally, we include a theoretical prediction of the pressureinduced phase transition in BaZrO 3 . An independent experimental confirmation of the stability of the cubic ground state of BaZrO 3 is given by the observed pressure-induced phase transition at 17.2 GPa at 300 K [24]. The phase transition was hypothesized to be from cubic to the tetragonal I4/mcm symmetry, similar to the antiferrodistortive phase transition exhibited by SrTiO 3 at 105 K and ambient pressure. The SrTiO 3 phase transition has been studied theoretically as a function of both temperature and pressure [25,26]. Here we deliver a zero-temperature theory description for BaZrO 3 , confirming that the transition happens to the tetragonal I4/mcm symmetry by tracking both relaxations and the FIG. 1. BaZrO 3 structure and illustration of the Glazer angle rotation φ, taking the cubic Pm3m structure into the tetragonal I4/mcm structure. The Ba-O distance u is also marked. Note that the barium and oxygen ions are in different planes.
R-mode vibrational frequency under compression. We note that accuracy of phonon descriptions [22,23,27] is a strong discriminator of modeling quality since the lattice dynamics directly reflects the electron density variation, an essential quality measure of DFT performance [28]. We hope to inspire low-temperature experiments on the pressure-induced phase transition in BaZrO 3 for independent checks on our modeling.

II. UNDERSTANDING BaZrO 3
BaZrO 3 exhibits imaginary phonon modes at the R point in some functional approximations [11]. The unstable mode has the irreducible representation R 25 . It consists of essentially rigid rotations (often called tilts) of the inscribed oxygen octahedra in a cog-wheel fashion, with rotation in successive layers in the opposite direction. The oxygen octahedral rotations have been systematized by Glazer [29] and we term the rotation Glazer rotation. The tetragonal I4/mcm phase is formed after a Glazer rotation around one axis as illustrated in Fig. 1. By combining rotations around the three crystallographically equivalent Cartesian axes, different crystal symmetries can be attained. An equal rotation around two axes leads to the orthorhombic Imma symmetry and rotation around all three axes leads to the rhombohedral R3c.
An instability on the M point would lead to rotations in the same direction in successive layers. Since no functional exhibits instabilities on the M point, we have not considered these [30]. There are also three ferroelectric modes, associated with a -point instability, which we have considered but have found not relevant.

A. Structure, stability, and thermal expansion
Imaginary phonon modes are found as solutions in the diagonalization of the dynamical matrix when the potential energy landscape is not positive definite. When the configuration represents a saddle point in the potential energy landscape, a lower energy configuration is attainable by displacing the atoms along the unstable phonon mode. This stability criterion of positive definite, however, is only a sufficient condition. It is not a necessary condition since it does not take quantum fluctuations into account. Quantum fluctuations can stabilize a phase with only a weak instability [31,32].
The simplest way of taking some ZPE effects into account is through the quasiharmonic approximation [33] (QHA). We obtain the ZPE-corrected lattice constant a ZPE and the ZPEcorrected vibrational frequencies ω ZPE , as well as temperature effects, within the QHA. As we have reported before [11], and is shown in Fig. 2, the vibrational frequency for the R 25 mode is very sensitive to the lattice constant. Figure 3 shows the full thermal expansion as computed in the QHA. The same calculations were also implicitly used in Ref. [11] to obtain the ZPE-corrected lattice constant. The thermal expansion shows an excellent agreement with the measurements in Ref. [4] for the CX0p hybrid functional. The PBE functional only exhibits nonimaginary phonon modes in a relatively narrow span of lattice constants from a = 4.235 to a = 4.335. At larger lattice constants, the PBE exhibits a phase transition into a tetragonal Amm2 phase. Although the lattice constant stays well within this interval for temperatures below 1000 K, the uncertainty in the prediction of the thermal lattice expansion for PBE must be considered large, especially at higher temperatures.
The CX functional requires a separate discussion. At its equilibrium lattice constant, the CX functional exhibits imag- inary phonon modes. To perform thermal modeling in the QHA on the CX, we only use lattice constants exhibiting real frequencies, i.e., larger than 4.2073 Å. This corresponds to a temperature of about 50 K in the QHA. In Fig. 3, we only show values for CX at temperatures above 125 K, when a > 4.208 Å. Zero-point motion helps stabilize the cubic phase and we treat BaZrO 3 as cubic also in CX.

B. Dielectric constant
The dielectric constant has been measured accurately, both as a function of temperature and as a function of frequency. Measurements, performed at both low frequency and low temperature, set the dielectric constant values as 47 [6] and 43 [34]. Previous theoretical works [10] have failed to obtain this value by means of DFT. In cubic symmetry, the dielectric constant of BaZrO 3 has been calculated to a value of 65 using LDA [10], which is significantly higher than the experimental values. By proposing a lower symmetry structure, a dielectric constant of 50 was obtained [10] for LDA, a value much more in line with experiments. This explanation is, however, unsatisfactory as no phase transition away from the cubic high-symmetry phase has been observed. We have confirmed the cubic nature of BaZrO 3 in Ref. [11] using hybrid functionals in combination with inelastic neutron scattering experiments. The remaining question is if the same functionals that predict the cubic nature of BaZrO 3 can also resolve the question of the value of the dielectric constant.

C. Mean-square relative displacement
The fluctuations in the relative distance between two atoms in a crystal can be measured using EXAFS [35]. We compute these parallel MSRD values u 2 for the set of atom pairs. The distance u for the first-neighbor shell relative to barium is marked in Fig. 1. The MSRD differs from the ordinary mean-square displacement (MSD). While measured MSD values only reflect the fluctuations of an atom around its equilibrium position on a macroscopic level, the MSRD values give detailed information on a local level.
As a rule of thumb, u 2 increases with increasing interatomic distance. However, this does not seem to be the case for increasing distances relative to the barium atom in BaZrO 3 [9]. Here the first-neighbor shell exhibits significantly higher fluctuation than both the second and third interaction shells. This is to some extent expected if the Zr-O octahedron is assumed rigid. The R 25 mode has a very low frequency and is easily populated with increasing temperature. However, when the octahedron rattles, the Ba-Zr distance, as well as the Zr-O distance, will fluctuate significantly less than the Ba-O distance.
Lebedev and Sluchinskaya suggest, based on an LDAbased analysis of the EXAFS findings of large BaO u 2 values, the presence of actual deformations off of the cubic structure. The deformations are assumed to take the form of rigid rotations of the octahedrons at angles up to 4 • [9]. For this account of the so-called anomalous Debye-Waller factor, it is also assumed that the closeness in phase enthalpies makes BaZrO 3 form a glass state [9]. The assumed rotations are larger than the 2 • exhibited by SrTiO 3 at 50 K, that is, the value that holds well below the SrTiO 3 phase-transition temperature.
Since we have shown [11] that BaZrO 3 indeed remains cubic, an outstanding question is if the EXAFS data can instead be accurately described within the cubic phase by using one of the hybrid functionals.

D. Phase transitions
Contrary to SrTiO 3 , there is no observed temperatureinduced phase transition for BaZrO 3 [11]. There is, on the other hand, an observed pressure-induced phase transition for BaZrO 3 at 17.2 GPa [24] and 300 K into what is hypothesized to be a tetragonal I4/mcm phase. In this paper, a major focus is to predict the transition pressure at zero temperature. Figure 4 shows the energy-volume curves of the cubic phase and the three phases attainable from an R 25 -mode instability. The curves for the different phases follow almost on top of each other with very small energy differences. The inset shows the difference between the different phases and the cubic phase. The shaded area in Fig. 4 indicates the average magnitude of the observed energy differences between the cubic phase rendered in different simulation cells. Energy differences on this scale are usually considered smaller than the accuracy of DFT. We have therefore taken extra precautions by a symmetry adaptation of the phases, as described in Sec. IV below.
The important observation in Fig. 4 is that we find no discrete volume changes between the respective minima of the different phases. This suggests that the phase transition is of second order. However, Fig. 4 is computed for PBE which shows cubic-phase stability only because PBE also overestimates the cubic-lattice constant, Fig. 2. A more complete characterization of the pressure-induced phase transition is motivated.
FIG. 4. Energy as a function of volume (left) and enthalpy as a function of pressure (right) for the three different phases attainable from an R-mode instability computed using the PBE functional. Here the energies and volumes are given relative to the reference values E 0 and V 0 defined as the energy minimum and corresponding volume of the cubic phase.

A. Structure, stability, and thermal expansion
A sufficient condition for stability of a material phase is that the potential energy landscape be positive definite. For BaZrO 3 , the most crucial direction is that of the R 25 mode. The potential energy surface along the R 25 mode can be described by an energy expansion, where E is the energy per formula unit (five atoms) of BaZrO 3 , ω R is the vibrational frequency, and Q R is the generalized phonon coordinate for the R 25 mode. Q R can be related to the Glazer tilt angle φ for the rotation of the oxygen octahedra around the [001] axis (cf. Fig. 1 where m O is the mass of oxygen and a is the lattice constant.
The stability criterion for the cubic phase effectively means that ω 2 R > 0. That is, we check if the R-mode frequencies are real and positive when computed in DFT, using the standard assumptions of classical nuclear dynamics.
By varying the volume of the unit cell, assuming the harmonic approximation is valid at each volume, the thermal expansion can be obtained through the QHA in which the volume dependence of the vibrational free energy is accounted for through F (T, V ) ≈ U el (V ) + F vib (T, V ). The Gibbs free energy is obtained by minimizing the free energy with respect to the volume, As a result, the volume as a function of temperature V (T ) and the temperature-dependent frequencies ω(V (T )) can be extracted.

B. Dielectric constant
The dielectric constant is an interesting quantity as it is sensitive not only to the vibrational frequencies ω s but also to the eigenvectors e s,iα , and thus to the entire dynamical matrix 224105-4 at q = 0. These quantities can be readily obtained from a phonopy calculation [33]. In the eigenvector e s,iα , in Eq. 2 and below, i denotes the atom, α (or β) the Cartesian direction, and s is the band index. The ionic contribution to the dielectric constant can be computed through [36] The Born effective charges Z * i,αα can be obtained from a linear response calculation or a Berry's phase calculation. The calculation of the dielectric constant is discussed in further detail in Appendix A.

C. Mean-square relative displacement
The MSRD values are sensitive to the dynamical matrix not only at the point but in the entire Brillouin zone. The MSRD can be computed [35] from a phonon calculation through where Here R ab is the vector connecting atoms a and b andR ab the corresponding vector of unit length, m a and m b are the respective masses, μ ab the corresponding reduced mass.

D. Phase transitions
Many phase transitions can be expressed in terms of a continuous order parameter, often denoted q, which is a quantitative measure of the extent to which the phase transition has changed the structure [37]. A first, and necessary, requirement for this formalism to be applicable is that the two phases have a group-subgroup relation, i.e., all the symmetry elements in the lower symmetry structure are contained in the higher symmetry structure [1]. The phase transition from cubic Pm3m to tetragonal I4/mcm is an example of a continuous phase transition.
For a continuous phase transition, the free energy can be expressed as an expansion around q = 0: When the temperature is increased from below to above T C , the double well potential shape is replaced by a parabola and the higher symmetry phase with q = 0 as the equilibrium state is stabilized. This formalism often goes under the name of Landau expansion or excess free energy [37,38] and describes the increase (decrease) of the free energy. The Landau-theory formalism has been applied successfully to the cubic-tetragonal phase transition in SrTiO 3 [39], where it has also been used to determine the phase diagram under both pressure and temperature [25,26].
The term Landau expansion is often reserved for the temperature dependence of the excess free energy. However, since the temperature and the pressure occur on equal footing in the free energy, the excess free energy can also be written in terms of the pressure [25,26,40]: In fact, the transition pressure is proportional to the transition temperature [40]: Here K is the bulk modulus and λ * 2 is the (renormalized) coupling constant between the order parameters and the elastic strain. By defining c 0 = 1 2 aK/λ * 2 , the zero-temperature critical pressure can be written as P 0 = −c 0 T C . Equation (7) can thus be rewritten: At equilibrium ∂G ∂q = 0, the squared order parameter is given by Positive solutions indicate finite values of the order parameter and hence an instability of the cubic phase. Negative solutions result in imaginary order parameter values and indicate stability for the cubic phase. Note that we here view the problem from the tetragonal side. In the limit C → 0, the phase transition is termed second order and the equilibrium order parameter is given by We use three discriminators for the phase transition: (1) The Glazer or tilt angle φ R . This is perhaps the most natural order parameter. The pressure at which φ R reaches zero marks the critical pressure P C .
(2) The tetragonal strain: Here e i = a i /a 0 is the strain in the Cartesian direction i. The tetragonal strain e t measures the deviation from cubic symmetry and will reach zero at P C .
(3) Minus the square of the R 25 mode frequency −ω 2 R 25 , see Fig. 2. This is equivalent to measuring the curvature of the potential energy landscape. A value −ω 2 R = 0 identifies P C . For both (1) and (2), the critical pressure is determined by allowing strongly distorted configurations to relax fully, subject to the constraint of an increased pressure. The equilibrium tilt angle and tetragonal strain is then obtained from the relaxed configuration. We therefore treat these two discriminators together, in Sec. V D below.
The use of the R 25 -mode frequency as a discriminator is rationalized as follows. The Glazer angle is very closely related to the normal coordinates Q sq , which diagonalizes the Hamiltonian. In the harmonic approximation, the full potential is approximated up to second order by By diagonalizing the dynamical matrix, the harmonic potential can then be written as Expressing the Landau potential to second order, with the order parameter chosen as the normal coordinate of the R 25 mode (q ≡ Q sq | sq=R 25 ≡ Q R 25 ), we obtain the relation That is, the square of the harmonic vibrational frequency is directly proportional to the applied pressure: Pressures larger than P C lead to negative ω 2 R 25 and instability of the cubic phase. The ω 2 R 25 discriminator has the advantage that it can be obtained directly from a phonon calculation in a cubic cell.

IV. COMPUTATIONAL STRATEGY AND DETAILS
DFT calculations were performed using the projector augmented wave (PAW) method [41,42] as implemented in the VASP [43,44] software. Out of the six different approximations to the exchange-correlation (XC) functional used in our previous work [11], we here retain only the four which predict a stable cubic ground-state structure for BaZrO 3 . Figure 5 illustrates how the remaining four span a twodimensional space with varying degrees of nonlocality in the description of exchange and correlation effects.
The semilocal GGA in the constraint-based PBE version [45] can be viewed as a starting point. Truly nonlocal correlations are represented by the consistent-exchange vdW-DF-cx functional [12,17,18], here abbreviated CX. This functional is a version of the vdW-DF method [12][13][14]. It balances exchange and correlation by use of the Lindhard screening logic [13,46].
Nonlocal Fock exchange is included in two hybrid functionals. The vdW-DF-cx0p [16], here abbreviated as CX0p, is a hybrid extension of the CX functional [47] and includes both nonlocal correlation and nonlocal Fock exchange. The fraction of Fock exchange is fixed, set at 20% following a coupling constant scaling analysis [48] of binding in sparse matter [16].
HSE [15], (sometimes called HSE06) is a range-separated hybrid extension of the PBE functional. It uses 25% Fock exchange for the description of the short-range Coulomb interaction. The range separation is described by a screening parameter μ = 0.2 Å −1 and an error-function weighting erf (μr)/r of the Coulomb interaction [15].
The effect from the nonlocal exchange and nonlocal correlation can be analyzed separately. This is done by comparing how well this set of functionals (Fig. 5) performs in characterizations of known BaZrO 3 properties. Furthermore, by here making predictions (of the T = 0 transition pressure), we invite measurements that could provide a strong validation of the functionals.
The inclusion of nonlocal exchange can sometimes significantly improve the XC functional approximations. Similarly, the inclusion of a truly nonlocal correlation is important in sparse matter. BaZrO 3 is not sparse and actual van der Waals forces may be believed less important. However, the vdW-DF method captures more general nonlocal screening effects and the consistent versions (CX and CX0p) are designed to balance exchange and correlation [13,16,17,46,49,50]. Figure 6 shows that there are differences as we turn on the description of nonlocal correlation effects. Figure 6 illustrates the changes that arise in the dynamical matrix with the inclusion of first nonlocal exchange and then nonlocal correlation. The change in the dynamical matrix is significantly larger when nonlocal correlation effects are included.
The calculations were converged with respect to the R 25mode frequency. This 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 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. The hard PAW potential has been used throughout this paper. The calculations with the semilocal PBE were conducted at an energy cutoff of 1200 eV, while for the functionals including nonlocal effects, CX and the hybrids CX0p and HSE, energy cutoffs of 1600 eV were required. 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.
There is one important factor to consider with the k-point sampling of different crystal phases. The energy differences between the cubic phase and the distorted phases are too small to be compared directly as we track continuous phase transitions. While the cubic structure can be rendered in a five-atom cell, the smallest possible cell in which the other phases can be represented contains ten atoms. This causes slightly different k-point meshes in the different setups. To compare the very small energy differences (of the lowersymmetry phases and the cubic phase), the cell has to be symmetry adapted. That is, when comparing the cubic and a potentially competing lower-symmetry phase, we treat the cubic system in the unit cell that is natural for the lowersymmetry phase. In the lower-symmetry cell, we then adjust the parameters minimally to describe the cubic symmetry in that larger cell. This means that the k-point meshes become nearly identical.
The phonon spectra were all calculated using the frozenphonon method with a displacement of 0.01 Å in 2 × 2 × 2 supercells containing 40 atoms and postprocessed in phonopy [33]. A 10 × 10 × 10 k-point mesh was used for sampling of the Brillouin zone for all phonon calculations.

V. RESULTS AND DISCUSSION
In this paper, we assert the fundamental performance and the prediction accuracy of DFT by looking at a range of measurements available for BaZrO 3 . Below we present and discuss our comparison. Table I shows the equilibrium lattice constants obtained from Birch-Murnaghan fits. We have confirmed the implied precision for lattice constants as 0.001 Å for HSE (better for PBE/CX/CX0p). This was done both by varying the fitting scope (what range of lattice constants are included) and by testing the Birch-Murnaghan results against those obtained using a fourth-order polynomail fit, ensuring a focus on fitting the minimum of the potential-energy landscape [56].

A. Structure and stability
The computational values agree well with experiments for all functionals except the PBE, which exhibits the well-known overestimation by about 1%. We obtain the ZPE corrected lattice constants a ZPE by taking anharmonic effects into account through the QHA.
By applying the ZPE correction obtained from the QHA, the lattice constant increases by about 0.007 Å. The ZPE correction for CX could not be directly obtained since CX exhibits imaginary modes at the equilibrium lattice constant, as discussed in connection with Fig. 3. Table I also shows the vibrational frequencies at the point, and generally shows good agreement with the experimental values. Once again, PBE is less accurate, and generally predicts lower vibrational frequencies. The PBE problems coincide with the overestimation of the lattice constant, see Fig. 2.
The difference between CX and PBE can to a large extent be attributed to the different equilibrium lattice constants. Calculated at the experimental lattice constant [11], a = 4.188, the differences in the vibrational frequencies are negligible, with the exception of the R 25 mode, which shows a slightly larger dependence on the functional.
We note as an interesting aside that the 25 mode is very insensitive to the lattice constant. The silent 25 mode has atom motion perpendicular to the bond direction. It crosses the TO 2 mode at about a lattice constant of 4.200 Å such that their mutual order is changed at smaller lattice constants. Overall, the -point frequencies show a smaller variation with functional than the R point. It seems that the softer the frequency, the less accurate PBE will be.
The immediate effect of the inclusion of Fock exchange is to reduce the lattice constant. This increased stiffness in both HSE and CX0p naturally causes a general hardening of the frequencies. There is also an additional hardening of the modes containing motion mainly perpendicular to the bond direction, i.e., for the modes 25 (silent), R 25 (AFD), R 15 (scissor), but also the 15 modes TO 1 and TO 2 . This has been observed before [57,58] and is attributed to the better electron density obtained by the inclusion of Fock exchange [28].

B. Dielectric constant
Table I furthermore summarizes our results for the dielectric constant. Compared with the experimental data from Refs. [6,34], performed at low frequency and low temperature, the PBE functional performs poorly. The same is also true to some extent for CX. The reason behind this is not only the accuracy in the computed vibrational frequencies. It is also a problem that reflects the quality of the eigenvectors obtained 224105-7 from diagonalizing the dynamical matrix. This is discussed further in Appendix A. The Born effective charges and the electronic contribution to the dielectric constant differ very little between the different functionals. We find that the difference in the overall dielectric constant between the different functionals can almost entirely be attributed to the difference in the ionic contributions.
It should be noted that both hybrid functionals predict dielectric constants in a cubic symmetry in good agreement with the measured values. There is no need to assume a lower symmetry structure for a correct value for the dielectric constant.
Finally, it is interesting to note that the quite large differences in the dynamical matrix, illustrated in Fig. 6, do not result in any similarly large differences in the dielectric constant. In contrast, the inclusion of nonlocal exchange has only a small impact on the dynamical matrix but a significant effect on the dielectric constant that is derived from this dynamical matrix. Table II shows the computed parallel-MSRD results for the first-, second-, and third-interaction shells relative to barium and zirconium. It is compared with experimental values from Refs. [9,59,60]. At the barium edge, the MSRD values for oxygen u 2 BaO are about two to three times those of the second and third interaction shells. This is contradictory to the rule of thumb that the MSRD values should increase with increasing interatomic distance.

C. Mean-square relative displacement
Table II also shows that the unusually large MSRD value for BaO is a consequence solely of the harmonic approximation. It holds without the addition of any anharmonic effects, save for the ZPE adjustment of the lattice constant.
At elevated temperature, the calculated MSRD values can be compared with experimental values in Ref. [9]. As in the case of the dielectric constant, PBE overestimates slightly compared with experimental values. This is again due to the underestimated vibrational frequencies.
The calculated MSRD values are in good agreement with the experiments for the two hybrids. We have, in effect, demonstrated that the observation of cubic symmetry and of large MSRD BaO values are mutually consistent.
The temperature dependence for u 2 BaO is very strong compared to that of u 2 ZrO . This can be viewed as a confirmation of the impact of the R 25 mode. This mode is very easily populated as the temperature increases, due to its very low frequency [61]. The oxygen octahedra itself remains rigid. This causes the Zr-O distance to remain essentially unchanged while the Ba-O fluctuations increases with temperature.
For completeness, a comparison of measured [59,60] and computed temperature variations in the zirconium edge is also included. Here, in contrast to the MSRD values for the barium edge, all functionals perform well. We also see that the u 2 BaBa and u 2 ZrZr values are very similar but that u 2 BaBa has a much stronger temperature dependence than u 2 ZrZr . We attribute this to the lower frequency of the Last mode [51], with barium motion against the lattice, becoming populated at lower temperatures, while the zirconium motion is mainly found in the Slater mode [52], which has a higher frequency. Figure 7 compares the enthalpy variation with pressure for the set of phases that can result from an R-mode instability. The inset shows the enthalpy differences of these phases related to the symmetry-adapted cubic phase. The shaded area in the inset of Fig. 7 indicates the average magnitude of the observed energy differences between the cubic phase in the fiveand ten-atom simulation cells, respectively. The enthalpies are all calculated relative to their symmetry-adapted reference level.

D. Pressure-induced phase transition
A phase transition at T = 0 between two phases is often determined by tracking the intersection of the enthalpies H = U − PV . As can be seen in Fig. 7, there is no clear intersection for BaZrO 3 . At high pressure, all three lower symmetry phases are more stable than the cubic phase, with the tetragonal I4/mcm phase being the most stable for all functionals. As the pressure is reduced, the energy difference between the tetragonal and the cubic phase decreases until it finally disappears completely. The energy of the lower symmetry phase has become identical to the energy of the cubic phase. This can be understood by a symmetry analysis: the lower symmetry TABLE III. Phase transition pressures in GPa, computed at T = 0 K, using the three different discriminators: the Glazer angle squared φ 2 R , the tetragonal strain e t , and the R 25 -mode frequency squared ω 2 R , using the different functionals. The value marked by an asterisk is determined using a 2-4-6 potential, Eqs. (6) and (9) phases simply relax into the cubic phase, except at pressures larger than the critical pressure P C .
To further confirm that the cubic phase is the most stable phase at zero pressure, we performed full relaxations of structures in various symmetries starting from various pressures and various Glazer angles. The cubic structure was obtained as the zero-pressure structure in all cases, albeit after very many ionic-relaxation steps.
Since there is no clear intersection, a simple determination of the transition pressure becomes problematic [62]. However, the phase transition between the cubic and tetragonal phases can be treated as a continuous phase transition [1,25,40,63,64]. Since the tetragonal I4/mcm phase is the most stable of the lower symmetry phases, we will discuss only phase transition from the tetragonal phase to the cubic. Table III summarizes our predictions of the phasetransition pressure for the set of investigated functionals. We generally provide and cross-check three discriminators, see Sec. III D. However, for CX0p we were only able to converge an analysis in terms of discriminator 3, i.e., the R 25 -mode frequency squared. There is full agreement between the three discriminators for PBE, CX, and HSE (to be discussed below). We expect our prediction of the CX0p phase-transition pressure to be accurate even if we only use one discriminator. Figure 8 shows the determination of the transition pressure for discriminators 1 and 2. The calculations were done by setting up supercells in the tetragonal symmetry and letting these relax at constant volume. Since the tilt angle and the tetragonal strain are obtained from the same fully relaxed structure (and thus do not constitute independent discriminators), we henceforth focus on the tilt angle. For a purely second-order phase transition, the square of the order parameter, the Glazer angle φ, will be a linear function of the pressure which intersects the abscissa at the critical pressure P C . The tetragonal strain, being directly proportional to the square of the Glazer angle [40,63], is then also a linear function of the pressure.
For PBE (upper panel), the tilt angle does not exhibit a perfectly linear relation to the pressure. The data have instead been fitted with a 2-4-6 potential from Eq. (6), with the magnitude of the C parameter about 1000 times smaller than A. The tetragonal strain shows a more linear trend with increasing pressure for all functionals. For the hybrid functionals, full convergence was an issue and the uncertainty in the data is larger. Therefore, a straight-line fit has been chosen. The third discriminator, the square of the soft mode frequency, should exhibit a linear dependence of the pressure in a second-order phase transition, see Eq. (15). identical slope values. For PBE and CX, the slope is about −7.5 meV 2 /GPa. HSE and CX0p have a slightly less negative slope of −6.9 meV 2 /GPa.
The transition pressures determined from the three discriminators are found in Table III. The predicted phase transition pressure is closely related to the magnitude of the R 25 -mode frequency at zero pressure. The larger the R 25 -mode frequency at zero pressure, the larger the transition pressure. The R 25 -mode frequency is only barely positive for PBE. The phase transition pressure is correspondingly low, only about 0.5 GPa.
The hybrid functionals predict a stable cubic phase at absolute zero within a finite pressure range. This may to some extent also be true of a full PBE characterization of phase stability, i.e., when also including corrections from zero-point vibrational-energy effects as has previously been done for SrTiO 3 [25,26]. Still, our here-predicted transition pressure is significantly higher for the hybrids than for PBE. We have shown that both HSE and CX0p can predict a range of measured properties in BaZrO 3 , including the R 25 -mode frequency [11]. The HSE and CX0p predictions for the transition pressure will be in better agreement than the PBE predictions. In fact, we expect that the HSE and CX0p predictions will be in fair agreement with future measurements of the phasetransition pressure at low temperatures.
The existing experimental value for P C for BaZrO 3 is given at room temperature in Ref. [24]. Using Landau theory, the critical pressure P C and temperature T C are proportional, Eq. (8). We assume this relation to hold over a finite temperature and pressure range. With the critical pressure at zero kelvin P 0 taken from Table III, the parameter T C can be obtained by fitting the P C (T = 300 K) value to the experimental value of 17.2 GPa [24]. When positive, T C could be interpreted as the transition or Curie temperature at zero pressure. For HSE, we obtain T C = −131 K and for CX0p T C = −102 K. The set of T C , P C (T = 0 K), and P C (T = 300 K) values gives us a simple picture of the delineation between the cubic and tetragonal phases.
Our computationally based values for T C appear to be in good agreement with the values obtained in Ref. [6]. Akbarzadeh et al. investigated the dielectric function at a range of low frequencies and characterized the temperature variation [6]. Using the Barret relation, they obtained a value for the T C parameter of −114 K. However, we do not obtain any signature of a phase transition in our calculations of the dielectric constant in the perfect lattice. It is therefore unclear if the agreement can be taken as a corroboration of our T = 0 theory prediction. This is discussed further in Appendix C.

VI. SUMMARY AND OUTLOOK
We have computed the vibrational spectrum of cubic BaZrO 3 using four different approximations to the XC functional. The inclusion of Fock exchange as implemented in HSE and CX0p stiffens the atomic bonds and increases the force constants, in particular those relating to motion perpendicular to the bond. This leads to vibrational frequencies that are in good agreement with the experimentally reported frequencies. We have also computed a series of observable properties of BaZrO 3 in the assumed cubic symmetry and compared these to available experiments. The calculated values for the dielectric constant, the MSRD, the thermal expansion, as well as the R 25 -mode frequency are all well predicted by the hybrid functionals HSE and CX0p. Figure 10 shows a summary of the results in terms of the relative error in percent compared to the experimental values for each functional. Where two or more experimental results are available, the value first mentioned in our tables above has been used. Starting from the left, Fig. 10 shows the relative error for (1) the lattice constant a (multiplied with 10), (2) the R 25 -mode frequency ω R , (3) the dielectric constant , (4) the parallel MSD for the Ba-O bond u 2 BaO , and (5) the root mean square for the thermal expansion. Finally, the total rootmean-square averages 1 for the five previous bars are shown in brighter colors as an overall assessment for each functional. This gives an estimate of the average error for the computed properties.
Both nonhybrids perform rather poorly with overall relative errors of above 20%. The poor performance of PBE (at 48.1% relative error) is mainly due to its inaccuracy in the description of thermal expansion. If the thermal expansion is removed from the total score, PBE and CX end up with similar scores. The average performance of the HSE hybrid is also strongly affected on the relative error that it makes on the thermal expansion coefficient α, see Fig. 3. The average HSE relative error drops from 28.5% to 9.4% when α is removed from the overall assessment of HSE.
The CX0p functional performs well with an overall relative error of σ = 4.7%, dropping to 2.8% if the α assessment is omitted from the average. For several properties, the CX0p error cannot be discerned in the comparison, Fig. 10.
Overall, the values calculated within a cubic-symmetry assumption using either of the two hybrids are in good agreement with the experimental values, except for the HSE description of the thermal expansion coefficient. Hence, we conclude that measurements for the set of here-investigated BaZrO 3 properties are fully commensurate with the assumption of cubic symmetry. Past discrepancies between measurements and DFT calculations can be resolved within the assumption of cubic symmetry by the inclusion of Fock exchange in the XC functional.  Table IV shows that the value of the predicted dielectric constant differs greatly between the different functionals. This is partly because the vibrational frequencies are computed to different values using the different functionals. However, this is not the full explanation. Indeed, the calculation of the dielectric constant proves to be a good test for the accuracy of the different functionals since the computed value of the dielectric constant depends on the exact details of the full -point dynamical matrix. The static dielectric tensor can be computed as the sum of the electronic contribution ε (el) αβ and the ionic contribution ε (ph) αβ [36]. The electronic contribution as well as the Born effective charges Z * i,αα are obtained from a linear response calculation or a Berry's phase calculation. The calculation of the ionic contribution also requires the eigenvectors e s,iα and eigenfrequencies ω s of the dynamical matrix at q = 0. These quantities are obtained from a phonopy calculation [33].
We define the mode effective charge as [36] Here an additional sum over degenerate modes is implicit, i.e., if a mode is degenerate (two-or three-dimensional irreducible representation), they should be summed together. From the mode effective charge, we also define the mode dielectric tensor (ionic contribution), where 0 is the unit cell volume. The average mode dielectric constant (ionic contribution) is Here again, the sum over α also implies a sum over degenerate modes s. Finally, we compute the total dielectric constant as By symmetry, the dielectric response is a constant times the identity matrix for a cubic crystal such as BaZrO 3 . For PBE, the generally lower vibrational frequencies is one major reason for a too-high value of ε. However, it is not the only contribution. Table IV shows that while the electronic contribution ε ∞ and the Born effective charges Z * i are very similar for all functionals, the mode effective charges Z * sα are not. These are computed from Eq. (A1) without taking the vibrational frequencies into the picture. Since the Born effective charges Z * i are not the culprits (and since the masses remain identical between the functionals), the difference must stem from the eigenvectors. This conjecture is verified by the data in Table IV. The largest discrepancy between the functionals is found for the Last mode, which leads to a greatly overestimated Z * 1α in PBE relative to the other functionals. We find that the eigenvector for the Last mode overestimates the oxygen motion. This has an impact because of the small mass of oxygen, causing an overestimation of the mode effective charge, cf. Eq. (A1). In combination with the low frequency of the Last mode, the Z * 1α error affects the total dielectric constant in PBE. Thus it is important to determine both the mode effective charge and the frequency accurately. The hybrid functionals do that.

APPENDIX B: TEMPERATURE DEPENDENCE OF MODES AND DIELECTRIC CONSTANT
The frequencies in the QHA generally decrease relative to the raw DFT lattice constant value already at zero kelvin. This FIG. 11. The temperature dependence of (a) the dielectric constant, (b) the Last-mode frequency (TO1) computed using the quasiharmonic approximation. Values for CX are only given at temperatures above 125 K (see Sec. II A). Experimental values are taken from Refs. [6,34].
is because the ZPE effects cause larger lattice constants. Using the QHA, we find about a 0.007 Å increase at zero kelvin compared to the raw DFT value.
At higher temperatures, the frequencies generally decrease further as the lattice expands. An exception is the R 25 mode, which shows the opposite trend. Figure 5 of Ref. [11] shows the temperature dependence for the R 25 -mode frequency together with experimental values. There is close agreement with CX0p. The trend is a small increase in the frequency with increasing temperature. This is expected since a compression of the lattice decreases the R 25 -mode frequency, see Fig. 2. Figure 11(b) shows the temperature dependence for the lowest frequency at the point, the 15 (TO1) mode. Here the trend is decreasing frequency with increasing temperature. The experimental data on the temperature dependence of the frequencies are taken from Fig. 4 in Ref. [34]. In contrast to the calculated values, the experimental values for the 15 (TO1) increase rather than decrease. The difference in trends may mean that higher order effects are important. It is possible that one might need to rely on the so-called pseudoharmonic approximation [37], the self-consistent phonon approximation [65,66], or effective harmonic models [67] rather than the 224105-12 presently used QHA. It is also possible that the measurements are effected by the presence of defects. Further investigations are needed to resolve this.
The same discrepancy between theory and experiments are found for the temperature dependence of the dielectric constant in Fig. 11(a). This is natural since the dielectric constant is intimately connected to the -point vibrational frequencies. Since the calculated vibrational frequencies have a different trend than the experiments, so does the dielectric constant. Figure 12 shows the dielectric constant calculated with PBE at a range of lattice constants, including both the cubic and tetragonal phases. The transition from the cubic to the orthorhombic phase with Amm2 symmetry at a = 4.335 is marked by a significant increase in the dielectric constant. This is indicative of a ferroelectric phase transition at negative pressure.

APPENDIX C: PHASE-TRANSITION SIGNATURES IN THE DIELECTRIC CONSTANT
The AFD phase transition at a = 4.235 does not lead to a signature in the dielectric constant in the absence of defects. For the cubic crystal, this can be easily understood because only -point modes contribute to the dielectric constant, cf. Eqs. (A1) to (A4). A soft mode at a zone boundary does not affect the dielectric constant.
On the tetragonal side of the phase transition, it can be understood as follows. The R 25 mode is split into an E g and a A 1g mode. The frequencies of these modes go to zero at the phase transition and could potentially lead to a divergence of the dielectric constant. However, both these modes are also even and do not respond to a static electric field. Hence, no divergence of the dielectric constant is seen at the phase transition between the tetragonal and the cubic phase in ideal crystals.
We note that the observed dependence of the dielectric constant has been interpreted as an incipient ferroelectric phase transition [6,34]. Akbarzadeh et al. [6] and Helal et al. [34] obtained values for the parameter T C from the temperature dependence of the dielectric constant. The negative values of T C were interpreted as a signature of an incipient ferroelectric phase transition. We have not found stability for a ferroelectrically distorted phase at positive pressures in any of the four functionals.
The numerical agreement of measured and predicted T C values discussed in Sec. V D cannot be taken as a corroboration of our predictions for P C (T = 0 K). This is because the AFD phase transition studied in this paper does not exhibit a similar response in our modeling of the dielectric constant, shown in Fig. 12. It is possible that the differences arise from the presence of the defects discussed in Refs. [6,34], but such a discussion is beyond the present scope. Instead, we hope that our T = 0 characterization can inspire new experiments that might validate or falsify our prediction.