First-principles study of cubic alkaline-earth metal zirconate perovskites

Zirconate perovskites are known for their ionic conductivity. First principles study are performed to analyze the structural, elastic, electronic and optical properties of zirconate perovskites. The effects of alkaline earth metal cation substitution on their overall properties are further investigated. The calculations are performed using the projector augmented wave (PAW) within GGA-PBE and HSE06 formalism. The obtained results not only predict the properties in accord with the experimental ones, but also compare the efficacy of these two functionals.


Introduction
Perovskite materials exhibit numerous functionalities such as piezoelectric, ferroelectric, ferromagnetic and pyroelectric. They are utilized in photovoltaic cells, LEDs, superconductivity, colossal magneto-resistance and topological insulators [1][2][3][4][5]. Generally, oxide perovskites display good dielectric properties, halide perovskites are good photonic materials and chalcogenide perovskites exhibit applications in energy harvesting, solid-state lighting and sensing [6][7][8]. In recent decades, there has been keen interest in zirconate perovskites as these materials are known for good proton conductivity and thus have potential uses for fuel cells or hydrogen sensors [9,10]. Moreover, due to their high thermal stability, zirconate perovskites are often regarded as potential candidates for thermal barrier coating materials [11].
Significant experimental work has been done to study different types of zirconate perovskites, establishing that injecting small dopant into these materials can result in ionic conduction behaviour [12,13]. From a theoretical perspective, commonly first-principles based DFT calculations, the study of zirconate perovskites has been implemented from the beginning of the twenty-first century. In 2005, Terki et al [11] have studied the structural, elastic and electronic properties of BaZrO 3 and SrZrO 3 using Perdew-Burke-Ernzerhof (PBE) functional [14,15]. Likewise, Hou (2008) [16] and Stoch et al (2011) [17] have investigated cubic and orthorhombic phases of CaZrO 3 , using PBE functional. Similarly, Shawahni et al (2018) [18] have studied SrRhO 3 and SrZrO 3 compounds, using PBE functional. Furthermore, one can notice from the literature that there are numerous DFT based studies performed on different zirconate perovskites. Using PBE functional for investigation, especially of optoelectronic properties, is always dubious as this functional underestimates the band gap by more than 40% [19,20]. This is because standard functionals such as PBE have intrinsic selfinteraction errors which treat electrons to be more itinerant. Therefore, it is often not suitable to use PBE functional for localized systems such as defects, d and f block elements. To overcome such shortcomings, usually Heyd-Scuseria-Ernzerhof hybrid functional (HSE) [21] approach is incorporated with standard calculations. The HSE functional rectifies the self-interaction error to some extent by mixing some portion of Hartree-Fock (HF) exact exchange with the exchange part of PBE. DFT calculations using HSE can thus give more reliable results other than standard functionals [22], though, the choice of functionals along with the accuracy of calculations depend on both atomistic systems and the availability of computing resources. HSE functional is an improvement over another hybrid functional PBE0 [23,24]. Due to problematic convergence of long-range mix of HF exact exchange in PBE0, HSE utilizes the screening parameter that allows the short-range mix between non-local HF exchange with local PBE exchange. Here, we have used standard HSE functional, called HSE06, which incorporates 25% Hartree-Fock (HF) exact exchange α, screening parameter w of 0.21 Å −1 and 100% correlation energy from PBE [25]. The magnitude of the screening parameter is crucial in predicting the accurate electronic structure of a material. It should be noted that HSE06 functional converges to PBE0 and PBE whenever ω → 0 and ω → ∞ respectively. Besides these functionals, we have also tested another functional PBEsol, which is quite popular for structure prediction [26,27].
The principal aim of this paper is to perform a comprehensive study of the cubic phase of AZrO 3 (A = Ca, Sr, Ba). Their structural, elastic, electronic and optical properties are predicted not only by using PBE functional but also with HSE06 functional. Throughout this paper, it is believed that one can easily trace the computational difference between these two functionals. Further, the effect of variation of cation A on the overall properties of these materials can also be examined. To the best of our knowledge, the simultaneous investigation of calcium zirconate, strontium zirconate and barium zirconate perovskites using both PBE and HSE06 functionals is the first of its kind in the literature.

Computational details
The investigation of structural, elastic, electronic and optical properties of alkaline-earth metal zirconate perovskites was executed through Density Functional Theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP) [28,29]. The projector augmented wave (PAW) [30,31] method considers the valence states 3s 2 3p 6 4s 2 for Ca, 4s 2 4p 6 5s 2 for Sr, 5s 2 5p 6 6s 2 for Ba, 4s 2 4p 6 4d 2 5s 2 for Zr and 2s 2 2p 4 for O while the remaining core states are considered to be frozen with the ion's environment. The cut-off energy for plane-wave basis sets was set at 520eV (1.3 times the maximum cut-off energy) and the k-points grid for Brillouin-zone integration was set as 24 × 24 × 24 for all PBE calculations and 4 × 4 × 4 for HSE06 calculations. Only exception was made for the calculations of the dielectric functions where the k-points grid was set as 12 × 12 × 12. The criteria for the convergence tolerance of self-consistent calculations was set as 10 −6 eV for total energy and 0.015 eV/Å for force.

Structural properties
The cubic structure of alkali-earth metal zirconate perovskite AZrO 3 (A = Ca, Sr and Ba) with space group Pm3m is shown in supplementary material, figure S1. The initial structures of AZrO 3 are taken from the materials project databases [32]. The geometrical positions of the elements in the primitive unit cell are: A (1/ 2,1/2,1/2), Zr (0,0,0) and O (0,1/2,0), (0,0,1/2), (1/2,0,0). The obtained lattice constants and bond lengths computed using PBE and HSE06, along with their respective available literature values are shown in table 1. Both of these functionals overestimate the lattice constant. Lattice constants are also calculated using functionals PBEsol and PBE0. It can be seen that the PBEsol greatly reduces the overestimation of PBE, and can be used as an alternative to computationally expensive hybrid functionals. Moreover, the structure predicted by hybrid functional PBE0 indicates that the lattice parameters are insensitive to the screening factor. In the case of CaZrO 3 and SrZrO 3 , there is a significant difference in computed values with their respective available literature values. This may be due to the fact that these zirconate perovskites do not crystallize in cubic phase but rather in orthorhombic phase below room temperature [33]. Further, the larger size of cation A would have higher lattice constant as there is no effect of charge in the structure due to the same number of valence electrons. Since the DFT simulation is a zero temperature calculation, the cubic structural phase of AZrO 3 may not be dynamically stable at 0 K. This can be seen from the phonon dispersion curves and density of states curves, as shown in figures 1 and 2 respectively. These phonon spectra are calculated from PHONOPY [34], under harmonic approximation, by considering non-analytical term correction. The phonon calculations are done in a 2 × 2 × 2 super-cell, with a sampling mesh of 30 × 30 × 30. In the case of PBE functional, density functional perturbation method (DFPT) [35] is implemented whereas for HSE06 functional, finite displacement method (FDM) [36] is used for phonon calculations. It is ensured that the computed dynamical force constant matrix is symmetric. The presence of soft modes and negative (imaginary) frequency at the M and R points at the Brillouin zone in CaZrO 3 and SrZrO 3 indicates that the cubic structure is not their ground state structures. One can easily notice that the stability increases as the size of cation A increases. More accurately, the orthorhombic form transitions to cubic phase as the cation size A increases. Only BaZrO 3 shows the stable cubic structure at 0 K. Furthermore, it has been experimentally known that there is no phase transition in BaZrO 3 between 4 and 1200 K, indicating high stable cubic structure [37]. The negligible forbidden gap between acoustic and optical phonons illustrates that they are good for thermal barrier coating materials.

Elastic properties
The elastic properties of a crystal can be determined by computing its stiffness tensor [C ij ], which is actually the response of a material to external stress. The stiffness matrices C ij s [ ] of calcium zirconate, strontium zirconate and barium zirconate perovskites, as calculated using PBE and HSE06 functionals are shown in the Supporting information. The inclusion of HF exchange causes an increment in the values of the elements in stiffness tensor. As one can notice from the stiffness matrices that, in a cubic system, normal stress only produces normal strain and shear stress only produces shear strain. In contrast with the diagonal elements of the fourth quadrant, the diagonal elements of the first quadrant decrease and off-diagonal elements increase in the order of cation size A, as observed in these zirconate perovskites. This observation implies that the shear modulus/bulk modulus increases/decreases whenever the cation size increases. Using the Voigt-Reuss-Hill (VRH) [41,42] averaging scheme, the shear modulus G and the bulk modulus B are related to [C ij ] as [43], 3.4 ( ) where h, k, N A and ρ are Planck constant, Boltzmann constant, Avogadro number and density respectively; n and M are number of atoms and molecular mass per formula unit respectively; v m is the average wave velocity calculated from G and B. They are shown in table 2.
The mechanical stability of these zirconate perovskites can be known from the Born's stability criteria [52], i.e., C 11 > 0, C 44 > 0, C 11 > |C 12 | and (C 11 + 2C 12 ) > 0. Based on these criteria, the materials under study are mechanically stable. Similarly, the anisotropic factor A being less than unity for all zirconate perovskites conveys that the materials exhibit elastic anisotropy. The brittle/ductile nature of the material can be known from the Pugh criterion B/G and the Poisson's ratio [53]. If B/G exceeds the critical value of 1.75, the material shows ductile nature; otherwise it is brittle. Likewise, the Poisson's ratio v greater than 0.33 indicates the ductile nature of the material; otherwise, it is brittle. Moreover, according to Pettifor criterion [54], the materials develop covalent character and possess brittle nature for negative Cauchy pressure; otherwise, they exhibit metallic character with ductile nature. All these criteria points that strontium zirconate and barium zirconate show brittle nature. Except the Poisson's ratio criterion, the Pugh criterion and the Pettifor criterion indicate that calcium zirconate shows ductile character. The hardness test is important as it describes the material's ability to resist deformation. The computed Vickers hardness qualitatively conveys that the more brittleness the material, the more is its hardness. The Debye temperature correlates the highest normal mode of vibration. It can be seen from table 2 that the Debye temperature decreases as one goes down the group of alkali-earth metal. This is because the density of an atom increases down the group which in turn decreases the velocity of sound.

Electronic properties
The electronic band structures for calcium zirconate, strontium zirconate and barium zirconate perovskites, along with their density of states (DOS), are shown in figures 3 and 4. It can be seen that HSE06 functional predicts the band gap with better accuracy than PBE functional. The spin-orbit coupling (SOC) shows trivial influence on these band structures. The band gaps due to HSE06 are also calculated using the PBE-generated structures (PBE-g). The severe underestimation of band gap by PBE and PBEsol indicates that these functionals suffer intrinsic self-interaction errors tremendously. On the other hand, the hybrid functional PBE0 tends to overestimate the band gap. Though PBE0 and HSE06 utilize the same portion of the mixing parameter, it shows that the screening factor is responsible for affecting the band gap. One can easily notice that the band gap bears an inverse relation with the screening parameter. Furthermore, these zirconate perovskites contain 4d electrons and besides global exchange correction, the correlation effect does influence the determination of band gap. The correction of correlation errors along with the pre-knowledge of the amount of HF exchange and the screening parameter is therefore crucial for predicting the more accurate electronic structure of zirconate perovskites. All these zirconate perovskites show topological resemblance and possess indirect band gap with transition occurring at R − Γ symmetric points. However, the energy eigen values of the valence band are close to each other at M and R symmetric points. The computed band gaps along with their experimental values are shown in table 3. From the analysis of orbital contributions to DOS, it has been found that the valence band maxima (VBM) is dominated by O-2p states and the conduction band maxima (CBM) is dominated by Zr-4d states. These states also satisfy the selection rule for transition, that the change in angular momentum is unity, i.e., Δl = ±1. One can notice that the band gap is independent of the order of cation size A. Similarly, in order to know the nature and strength of bonds in these structures, the Born effective charges (BECs) [35,55] have been calculated and shown in table S1. The flow of a charge in an ion with respect to its static charge (SC) arises from the coupling between lattice displacement and electrostatic field and can be quantified as,

Optical properties
The electronic response of a material towards an incident photon can be described from its complex dielectric function ε(ω) = ε 1 (ω) + iε 2 (ω). Generally, the real part ε 1 (ω) represents the phase lag between the incident field and induced field due to polarization and the imaginary part ε 2 (ω) represents the measure of energy loss. The imaginary part ε 2 (ω) can be evaluated from the momentum matrix elements between the unoccupied and occupied electronic states as [60], where,  p ,  f k n ( ),  E k n are the momentum operator, Fermi distribution function and eigenvalue of state  k n respectively. The real part ε 1 (ω) can then be evaluated from the imaginary part ε 2 (ω) using the Kramers-Kronig relation as [61], The calculated ε 1 (ω) and ε 2 (ω) for calcium, strontium and barium zirconate perovskites using both PBE and HSE06 functionals are shown in figure 5. Both PBE and HSE06 computations show similar dielectric spectra and any variation can be attributed to the effect of short range electron-electron interaction. For the time dependent field, ε 1 (0) is the static dielectric constant. It is noteworthy that the variation of static dielectric constant is independent of the size of ion A and shows inverse relation to band gap as suggested by Penn model [62] as, where, S 0 is a dimensionless constant, ω p is a screened plasma frequency and E g is a Penn gap. Usually, the Penn gap is larger than a band gap, and refers to energy corresponding to the first reflectivity peak. The value of S 0 , as suggested by Penn himself, is 1 and the more accurate value, especially for homopolar semiconductors, is taken as 0.62 [63]. Recent study by Diego Julio Cirilo-Lombardo [64] has suggested the value to be 2/3. Our work on zirconate perovskites suggest the value to be 0.86. Table 4 shows the values of ε 1 (0), plasmon energy ÿω p , Penn gap E g and reflectivity peak E R . On the other hand, ε 2 (ω) is seldom negative and starts from zero and shows variation after the incident photon energy exceeds the threshold band width. The peaks of ε 2 (ω) are related to the interband transitions. From the knowledge of ε 1 (ω) and ε 2 (ω), the other optical properties such as reflectivity R, refractive index n, extinction coefficient k, energy loss spectra L and absorption coefficient α are calculated as [65][66][67], They are shown in figure 6. One can notice that not only ε 1 (0) but also R(0) and n(0) are independent of the size of cation A. The values of R(0) and n(0) indicate that BaZrO 3 shows higher metallic character for the timeindependent field. The reflectivity spectra R(ω) of all zirconate perovskites show minima at energies corresponding to the plasma edge. This decrease in reflectivity may be due to plasmonic excitation. Furthermore, our calculated static refractive index satisfies the relation n(0) = (ε(0)) 1/2 and its inverse dependence with band gap can be verified from Harve-Vandamme (H−V ) relation [68], Figure 5. Real and imaginary parts of dielectric function of calcium zirconate, strontium zirconate and barium zirconate perovskites, computed using functionals PBE and HSE06. Table 4. Calculated values of ε 1 (0), ÿω p , E g (Penn) and E R . where A = 13.6 eV and B = 3.4 eV are the constants. Utilizing H − V relation, the refractive indices are calculated with the aid of their experimental band gaps and then compared with the simulated values. They are shown in table 5. The peak of energy loss spectra L(ω) characterizes the plasma resonance. The absorption coefficient α is almost zero when the photon energy is below band gap, which is trivial for single particle model. At energy higher than plasmon energy, the absorption coefficient decreases to zero and refractive index approaches one, indicating the transparency of the material. All these materials show high absorption in the U.V region.

Conclusions
In summary, a comprehensive study of the cubic phase of CaZrO 3 , SrZrO 3 and BaZrO 3 was performed using both PBE and HSE06 functionals. It was found that the lattice constants of these zirconate perovskites are Figure 6. Reflectivity, refractive index, extinction coefficient, energy loss and absorption coefficient spectra. dependent on the size of their respective alkaline-earth metal cations. Moreover, except BaZrO 3 , the cubic phases of CaZrO 3 and SrZrO 3 are high temperature structures. The study of mechanical properties reveals that SrZrO 3 and BaZrO 3 show brittle nature whereas CaZrO 3 shows some ductile character. All these materials show high indirect band gap with A − X dominated by ionic bond and B − X by covalent bond. The optical absorption of these zirconate perovskites are higher in the U.V. regions. Finally, it is anticipated that this study will be helpful to tune the properties of these materials by altering the cation A.