Evaluation of thermodynamics, formation energetics and electronic properties of vacancy defects in CaZrO3

Using first-principles total energy calculations we have evaluated the thermodynamics and the electronic properties of intrinsic vacancy defects in orthorhombic CaZrO3. Charge density calculations and the atoms-in-molecules concept are used to elucidate the changes in electronic properties of CaZrO3 upon the introduction of vacancy defects. We explore the chemical stability and defect formation energies of charge-neutral as well as of charged intrinsic vacancies under various synthesis conditions and also present full and partial Schottky reaction energies. The calculated electronic properties indicate that hole-doped state can be achieved in charge neutral Ca vacancy containing CaZrO3 under oxidation condition, while reduction condition allows to control the electrical conductivity of CaZrO3 depending on the charge state and concentration of oxygen vacancies. The clustering of neutral oxygen vacancies in CaZrO3 is examined as well. This provides useful information for tailoring the electronic properties of this material. We show that intentional incorporation of various forms of intrinsic vacancy defects in CaZrO3 allows to considerably modify its electronic properties, making this material suitable for a wide range of applications.

conductivity at elevated temperatures. In fact, In 2 O 3 -doped CZO has been practically used as an electrolyte of galvanic cell-type hydrogen sensor for molten metals [20][21][22] . Similarly, several reports suggest that disturbing the stoichiometric of CZO can bring about novel mixed p-type and ionic conduction behavior that can be controlled by varying the synthesis conditions [23][24][25][26] .
Although a lot of experimental work has been undertaken for studying a wide range of properties of CZO, fewer theoretical studies of calcium zirconate are available in literature. First-principles calculations for investigating the structural, mechanical, and electronic properties of cubic and orthorhombic phases of CZO have been carried out by Hou 27 and Stoch et al. 28 , while Brik et al. 29 studied the (001) surfaces of CZO for exploring their electronic properties and energetic stability. On the other hand, the potential for using Yb 3+ -, Nd 3+ -, In 3+ -, Ga 3+and Sc 3+ -doped CZO for protonic conduction has lead many researchers to investigate the mechanisms of doping and trapping of H + ion in CZO using quantum mechanical techniques [30][31][32] . Since the intentional incorporation of intrinsic vacancy defects in doped transition metal perovskite oxides promises profound enhancement in their ionic conductivity 30,[33][34][35] , it is imperative to have a thorough theoretical insight in to the chemistry of intrinsic vacancy defects containing un-doped CZO under various growth conditions 3,36,37 . The motivation for present study is further strengthened by the fact that tuning the electronic transport properties of CZO by manipulating vacancy defects has not been carried out so far. To this end, we employ the full-potential linearized augmented plane-wave (FP-LAPW) method within the framework of DFT for investigating the influence of vacancy defect on the electronic structure of CZO. It is expected that the detailed study of vacancy defects and clustering of charge neutral O vacancies reported in this work may stimulate future experimental studies to identify and use non-stoichiometric CZO for advanced device applications.

Method of Calculation
To evaluate the relative stability of V Ca q , V Zr q and V O q vacancies (where q represents the charge state) in CZO and their corresponding electronic properties, we employ the all-electron FP-LAPW method as implemented in the WIEN2k code 38 . Throughout the whole of this study the exchange-correlation functional is modeled by the Perdew, Burke and Ernzerhof (PBE) 39 generalized gradient approximation (GGA) parametrization scheme. The FP-LAPW method requires partitioning the crystal by non-overlapping muffin-tin spheres. These spheres are centered at the calcium, zirconium and oxygen sites and are given radii (R MT ) of 2.11, 1.99 and 1.80 (in a.u.), respectively. FP-LAPW basis functions are constructed from spherical harmonics inside the muffin tin spheres, connected to plane waves in the interstitial region in between. The size of the basis set is controlled by the plane-wave cut-off (K max ). In order to pinpoint a basis set size that ensures sufficient precision, we exploit the fact that a determination of atomic chemical potentials for formation energy calculations depends entirely on the precise computation of enthalpies of formation. Since enthalpies of formation are calculated as the difference of total energies of a compound and its constituent atoms in their standard reference states, we have tested the precision of our results by comparing the enthalpies of formation computed using different choices for the K max and for the k-mesh used for numerical integration. Our results indicate that using K max = 4.444 and 6 × 4 × 6/4 × 4 × 4/12 × 12 × 12 k-meshes for CZO/ZrO 2 /CaO the precision in the calculated enthalpies of formation is ±1 meV/ atom. Therefore, the total energy calculations of different sizes of supercells used in the present work are performed using R o × K max = 8 (where R o is the R MT of O atom), while the maximum values of angular momentum of partial waves inside the muffin-tin spheres, l max , and the magnitude of vector for Fourier expansion of charge density, G max , are set at 10 and 18, respectively. As shown in ref. 40, all modern DFT codes make essentially identical predictions for properties that depend on total energies. The conclusions obtained in this work are therefore not affected by our choice for the FP-LAPW method.
CZO crystallizes experimentally in an orthorhombically distorted perovskite structure (space group # 62, Pnma) with 20 atoms and lattice parameters a exp. = 5.594 Å, b exp. = 8.021 Å and c exp. = 5.761 Å 41 , as shown in Fig. 1. This experimental unit cell has been used as the starting point for a complete DFT-GGA optimization of unit cell volume, c/a and b/a ratios and internal geometry using a 6 × 4 × 6 k-mesh. The GGA optimized cell is subsequently used for a 2 × 1 × 2 supercell with composition Ca 16 Zr 16 O 48 (80-atoms), which is sufficiently large for formation energy calculations and electronic properties of defected CZO 42 . A 3 × 4 × 3 k-mesh was used for this supercell and the self-consistent procedure to find the ground state density and the self consistent cycle was terminated when the difference in calculated energies between subsequent iterations was less than 10 −5 Ry. For all the defective supercells it was ensured that the atomic positions are fully relaxed, indicated by forces on each atom being below 1 mRy/a.u. The k-meshes for vacancy containing supercells of CZO adopting different symmetry structures and the basis set sizes for elemental solids are chosen with reference to the values of k-mesh and K max used for performing calculations for bulk unit cell of CZO 38 . The effect of spin-orbit coupling was neglected, after having tested it gave negligible effects on the energetic and electronic properties. In

Results and Discussion
Stability Diagram. For studying the chemical stability of pristine CZO, we have computed the limits of atomic chemical potentials (Δ μX ≤ 0) by assuming that the chemical potential (μ X ) of an isolated atom of species X is less than or equal to the chemical potential of that atomic species in its stable solid/gaseous state (µ X solid gas / ). The determination of valid limits of Δ μX is necessary for computing stability ranges shown in Fig. 2 and enables one to compute formation energies of vacancy defects under different growth conditions 45 . For stable production of CZO without the presence of secondary phases like CaO and ZrO 2 , we vary the values of Δμ Ca , Δμ Zr and Δμ O by satisfying the following equations 46 . The Gibbs free enthalpy, G, of crystalline solid can be distinguished from its total energy, E, through the term −TS + pV. Since the DFT calculations are performed at T = 0 and p = 0, the formation energy of a solid (e.g ∆E f CaZrO 3 ) coincides with its enthalpy of formation (e.g ∆H f CaZrO 3 ). In above equations, E's are the calculated minimum total energies of GGA optimized crystal structures of orthorhombic calcium zirconate, cubic calcium oxide, monoclinic zirconia, face centered cubic (fcc) calcium and hexagonal close packed (hcp) zirconium. The total energy of an O 2 molecule, on the other hand, has been computed by adding the GGA calculated cohesive energy of oxygen 47 to twice the total energy of a single oxygen atom inside a sufficiently large unit cell using only Γ-point sampling. For the sake of analyzing the error bars in the present DFT calculations, the cohesive energies, E c , of Ca and Zr are also computed. The calculated bare PBE equilibrium volume per atom (V PBE-bare ), zero point corrected equilibrium volume per atom (V corr. ), lattice parameters corresponding to corrected equilibrium volume (a corr. , b corr. and c corr. ), E c and ΔH f are presented in Table 1. Except for the case of fcc calcium, one can easily notice that all V PBE-bare values are overestimated, while the calculated PBE E c and ΔH f are underestimated. For the stable crystal structures of the compounds and elemental solids under study, V corr. have been computed by considering a systematic GGA deviation of 3.6% in the calculated volume and using zero-point correction 48 , ΔV, of 0.010 Å 3 / atom, 0.016 Å 3 /atom, 0.006 Å 3 /atom, 0.032 Å 3 /atom and 0.007 Å 3 /atom for CZO, CaO, ZrO 2 , fcc Ca and hcp Zr, respectively. Compared to experimental data, all the V corr. and the cohesive energies reside well within the intrinsic error bars proposed by K. Lejaeghere et al. 49 . It can also be seen that satisfies the fundamental criteria for the stable production of CZO 50 .
The chemical stability diagram (CSD) of CZO ( Fig. 2) has been obtained from the calculated values of ΔH f presented in Table 1 and by satisfying Equations 1-3. Stable production of CZO can be achieved for all points lying inside the quadrangle ABCD and on the lines joining points A (−6.203, −11.145, 0.000), B (−6.521,  Pristine and Neutral Vacancy Defects Containing CaZrO 3 . In the ideal cubic perovskite structure the oxygen atoms form an octahedral and cuboctocahedral coordination with Zr and Ca atoms, respectively, with Zr-O and Ca-O having constant values of bond lengths throughout the cubic unit cell. Under ambient conditions the orthorhombically distorted perovskite-like structure of CZO is found to be stable, which persists up to 1900 °C 28 . As evident from the calculated equilibrium unit cell of CZO shown in Fig. 1(a), the ZrO 6 octahedra are tilted in the ac-plane, which results in a deviation of both bond length and Zr-O-Zr bond angle from their ideal values. Table 2 provides a comparison between the atomic positions, bond lengths and bond angles computed in this study, and results from previous experimental and theoretical work. In the GGA optimized unit cell of orthorhombic CZO we found the Zr-O-Zr bond angle to be 146.290° and 144.632° with the equatorial (O 2 atoms residing in the ZrO 2 layer) and the apical (O 1 atoms residing in the CaO layer) oxygen atoms, respectively. Due to this tilt of the ZrO 6 octahedra the ideal 12 fold coordination of Ca with O in cubic perovskite structure is reduced to an 8 fold coordination in the orthorhombic structure, as shown in Fig. 1(b). The Ca atom is bonded with 8 oxygen atoms having bond lengths that range from 2.345 Å to 2.882 Å. On the other hand, the bond lengths of the Zr atom with its octahedrally coordinated oxygen atoms range from 2.108 Å to 2.118 Å.
Upon introducing non-stoichiometry in CZO, in form of vacancy defects, the bulk coordinations of O 1 and O 2 atoms with Ca and Zr atoms are partially eliminated. In order to obtain a stable atomic configuration, it is necessary to minimize the non-zero forces acting on atoms neighboring the vacancy site. Figure 3 displays the positions of atoms around the vacant Ca, Zr and O sites before (coloured spheres) and after (black circles) the relaxation of the internal geometry of defective supercells of CZO. The off-centering of a black ring from a coloured sphere demonstrates the movement of atoms from their ideal atomic positions of pristine CZO when the structure is fully relaxed. Figure 3(a) clearly shows that in the case of a calcium or zirconium vacancy the neighboring Zr and Ca atoms, respectively, are attracted towards the vacancy site due to the eliminated repulsive electrostatic forces between the two cation sites. Contrary to that, the oxygen  atoms exhibit outward relaxation when a Ca or Zr site is vacated. It is also evident that the vacant oxygen site in Fig. 3(c) results in lowest movement of the neighboring atoms.
In order to understand the above-mentioned changes in atomic positions resulting from vacant cation and oxygen site, it is useful to compare the bonding properties of Ca/Zr atoms with the O atoms in pristine and intrinsic vacancy defect containing CZO. For this reason we have computed the 3-dimensional (3D) electron-density distribution in pristine orthorhombic unit cell of CZO which is displayed in Fig. 4(a). In addition, the effective Bader charges (e) 51 are also calculated for a quantitative analysis of boning properties. Table 3 presents the effective Bader charges for the Ca, Zr and O atoms of pristine CZO and the atoms neighboring a vacancy site in defective supercell. It is clear from Fig. 4(a) that the bond between Zr and O atoms has a strong covalent nature, while ionic character of bonding prevails in the Ca-O bond. The effective Bader charges listed in Table 3 support these From the calculated electronic band structure of pristine CZO (Fig. 4(b)), we find that both the CBM and valence band maximum (VBM) are located at the Γ symmetry point; such that the occupied O-2p orbitals at the upper edge of the valence band and unoccupied Zr-4d orbitals at the lowest edge of the conduction band are separated by a direct fundamental energy band gap E g = 4.003 eV. The calculated band gap is close to a previous GGA study 27 , however it is underestimated when compared with experimental value (1.497 eV less than the band gap of CZO reported in ref. 53). In spite of the underestimation of the band gap by the GGA functional, these results do not restrict us from a qualitatively exploration of the electronic properties of non-stoichiometric CZO 54 . Analysis of the partial density of states (PDOS) plots shown in Fig. 4(c) reveal that the upper valence band is predominantly made up of occupied anion 2p orbitals with only minor differences between the contribution of O 1 -2p and O 2 -2p orbitals. The valence band of CZO also shows hybridization of the O-2p and Zr-4d orbitals which increases on going to lower (more negative) energies. On the other hand, the conduction band of CZO is dominated by the unoccupied Zr-4d orbitals. In case of defective supercells, the incorporation of V Ca 0 or V Zr 0 in CZO results in acceptor-like defect levels just above the bulk VBM, having an unoccupied O-2p character (Fig. 5(a and b). It is clear from Fig. 5(a and b) that a maximum localization of charge is achieved at occupied O site near the V Zr 0 . This suggests that presence of electrons in acceptor-like defect levels (e.g. a − V Zr 4 vacancy) would raise these levels closer to the CBM of bulk CZO as compared to charged Ca vacancies. This finding is in accord with the data listed in Table 3 where e of oxygen atom near V Zr 0 shows large reduction. The band structures of Ca and Zr deficient supercell presented in Fig. 5 The minimum total energies of bulk CZO and V X q containing supercells are represented by E CaZrO 3 and E nX [ ], q respectively, n is the number of X atoms removed from a supercell and μ x is the chemical potential of the species X. For comparing the stability of different types of vacancy defects in CZO, we use the chemical potential coordinates defined in Fig. 2. The term E F + E VBM in Equation 7 references the Fermi level with respect to energy of the VBM in defective supercell 46 . Since finite size of the supercell under periodic boundary conditions makes the E VBM in defective supercell is different from E VBM of the pristine CZO, we have computed the average electrostatic potentials, V avg. , at the occupied sites of vacancy type atom in pristine (pr.) and vacancy containing (de.) supercell (far away from the vacancy site) and the difference of these electrostatic potentials (ΔV avg. = V pr.
avg. ) has been added into E VBM . Moreover, the formation energies of oxygen vacancies have been corrected for the underestimation of band gap resulting from the use of GGA functional. This has been accounted for by using the band gap correction 54 where difference of experimental and GGA band gap (ΔE g = 1.497 eV) times m (number of electrons) has been added into the calculated Ω nO [ ] q values computed using Equation 7. This simple bad gap correction in the formation energies calculated using semilocal functionals has been found to provide reasonable agreement between theoretically computed defect levels for isolated oxygen vacancies and experimental observations 46,54 . It is worth pointing out here that hybrid DFT functionals can improve the evaluations of electronic properties in wide band gap complex oxides. However, the extremely high computational costs of hybrid  58 , ε is the experimental dielectric constant of CZO 59 and L is edge of the structure used for simulating charged vacancy defects in CZO. The calculated formation energies at points A, B, C, D and X for neutral and fully-charged vacancy defects are displayed in Fig. 6. The formation energies of intermediate charge is less than the formation energies of all other vacancy defects for point A through point D, showing that the fully charged oxygen vacancy can be easily incorporated in CZO. The formation energies of fully ionized Ca and Zr vacancies, on the other hand, are found to be larger than their charge neutral counterparts. These larger formation energies of negatively charged cation defects can be attributed to their electron scavenger nature as these vacancies tend to lower the energy of the supercell. Contrary to that, the creatn of ionized oxygen vacancies adds additional electrons to the system, compensating the occupied donor defect levels shown in Fig. 5  and we have obtained the full and partial Schottky reaction energies (ζ) 61 . The average Ca-partial and Zr-partial Schottky ζ are found to be 2.632 eV and 3.643 eV, respectively, where larger values of Zr-partial Schottky ζ supports easy incorporation of Ca vacancy in CZO 26 . On the other hand, the full Schottky ζ is found to be 3.839 eV. Figure 7   transport properties of transition metal perovskite oxides can not be attributed to isolated intrinsic vacancy defects 64 . In these situations, other forms of native point defect (such as self-interstitial or anti-site defects) are the possible suspects, however, large distortion in the crystal geometry and the resulting decrease in chemical stability 63 are generally not responsible for enhancing/degrading electronic and transport properties of bulk CZO. To this end, clustering of oxygen vacancies in perovskite oxides has been found responsible for greatly influencing the electronic properties. In fact, transmission electron microscopy and complementary independence spectroscopy studies have confirmed the role of oxygen vacancy clustering in the unusual properties of transition metal perovskite oxides 65,66 . Since our results clearly show that V O 0 is the most favorable form of charge neutral intrinsic  vacancy defects in CZO at both point C and D, in this section we explore the possibility of realizing oxygen vacancy clustering in CZO under the extreme reduction condition.
We have calculated Ω nX [ ] q for oxygen vacancy clustering (as explained in Method of Calculation section) using Equation 7 at the extreme reduction condition (Fig. 8(a) Table 3 and the strong hybridization of O-2p and Zr-4d orbitals just below E F shown in Fig. 4(c), it is clear that increasing oxygen vacancies clustering in the Zr 8 O 16 layer would lead to increase in charge delocalization which causes more attractive interaction among the Zr atoms in the layer. Moreover, a comparison of the defect formation energies presented in Figs 6 and 8(a) allows us to speculate that clustering of charged oxygen vacancies in the Zr 8 O 16 layer of CZO can be realized. The above findings are encouraging in view of using acceptor-doped CZO for high protonic conduction where proton trapping 67 causing major hindrance in long-range conduction can be effectively reduced by the formation of charged oxygen vacancy-acceptor clusters 62 .
For exploring the electronic properties of oxygen vacancy clustering in CZO, we have computed the electronic DOS for all the oxygen vacancy clattering cases. The calculated electronic band structure and DOS for V O 1 (not shown here) reveal deep donor-like levels of these charge neutral oxygen vacancy clustering cases similar to the ones shown in Fig. 5(c) and (d). In contrast, this donor-like level shifts into the CBM for the case of V O 8 2 ( Fig. 9(a)). The shifting of this defect level into the conduction band results from large charge delocalization from the zirconium atoms residing in defective Zr 8 O 16 layer (Zr 1 : residing in the Zr 8 O 16 layer containing the V O 8 2 vacancy) as shown in Fig. 9(a). The confinement of the delocalized charge in V O 8 2 containing Zr 8 O 16 layer of CZO ( Fig. 9(b)), therefore, allows the Zr 1 -4d states to move across E F . By using a 2 × 2 × 2 supercell we found that the presence of occupied states in CBM persists for charge neutral V O 8 2 clustering in CZO and is independent of the increasing number of defect free Ca 8   containing CZO which causes Zr 1 -4d states to move across E F . Although we have only considered neutral oxygen vacancy clustering in the present work, an easier incorporation of charged oxygen vacancy clustering in CZO can also be deduced from the formation energies presented in Fig. 6. These findings support the use of acceptor-doped CZO for high protonic conduction where proton trapping can be avoided by intentionally clustering charged oxygen vacancies (e.g. + V O 1 ) in the Zr 8 O 16 layer of CZO 68 . The clustering of charged oxygen vacancies would certainly eliminate the charge delocalization and n-type nature of CZO shown in Fig. 9, however, this reduction clearly supports the decreased electronic conductivity and increased ionic conduction in Y 2 O 3 (etc.)-doped CZO which makes calcium zirconate a potential solid-state electrolyte materials for solid oxide fuel cell technology 23,24 .

Conclusions
In summary, we have employed first-principles calculations for investigating the influence of isolated V Ca q , V Zr q and V O q vacancies and oxygen vacancy clustering on the electronic structure of orthorhombic CZO. Our results reveal that pristine as well as low concentration of charge neutral oxygen vacancy containing CZO are insulating. On the other hand, charge neutral Ca and Zr vacancy containing CZO are found to be hole-doped systems, where low formation energies of Ca vacancy confirms the contribution of V Ca 0 in experimentally observed mixed p-type and ionic conduction behaviour of CZO. For the case of neutral vacancy defects, calcium/oxygen vacancies are the most abundant form of vacancy defects in CZO under oxidation/reduction condition. We find that fully charged oxygen vacancies are the most favorable as compared to all other types and charge states of intrinsic vacancy defects in CZO. The calculated values of Schottky ζ permit us to predict that non-stoichiometric CZO could allow tunable p-type conductivity. It is shown that a high concentration of O vacancies can be experimentally realized in the Zr 8 O 16 layer of CZO which highlights potential utilization of acceptor-doped CZO for high protonic conduction where proton trapping can be avoided by means of introducing charged oxygen vacancies around the dopant site. The wide range of possibilities in tailoring the electronic structure of CZO by means of intrinsic vacancy defects makes it attractive for electronic, electrical and optical devices.