First principles calculations of the thermodynamic stability of Ba, Zr, and O vacancies in BaZrO3

The temperature dependence of the stability of bulk BaZrO3 (BZO) and of the vacancies in this material are investigated by considering phonon contributions to the free energy. The stability diagram of BZO is determined for different chemical environments. With increasing temperature the stability region becomes smaller which is particularly caused by the strong temperature dependence of the chemical potential of gaseous oxygen. The free formation energy of Ba, Zr, and O vacancies in BZO is calculated for all possible charge states and for different atomic reservoirs. While the free formation energy of Zr vacancies is strongly influenced by temperature a weaker dependence is found for Ba and O vacancies. This also has an effect on the charge transition levels at different temperatures. The present results demonstrate that O poor reservoir conditions and a Fermi level close to the valence band maximum favour a high concentration of doubly positively charged O vacancies which is a prerequisite to get a large number of protonic defects and good proton conductivity. In such a chemical environment the number of Ba and Zr vacancies is low so that Ba and Zr deficiencies are not an important issue and BZO remains sufficiently stable.


Introduction
BaZrO 3 (BZO) is a prototypical perovskite-structured material with oxygen-ion and proton conducting properties. 1-3 BZO is a candidate material for various applications such as solid oxide fuel cells, gas sensors, and hydrogen pumps. [4][5][6][7][8] It shows good chemical and thermodynamic stability even at very high temperature. 9,10 BZO becomes proton conducting aer exposure to humid conditions. 10,11 The water from the gas phase dissociates into a hydroxide ion and a proton, the hydroxide ion lls a doubly positively charged oxygen vacancy, and the proton forms a covalent bond with an oxygen atom of the perovskite lattice. Finally, two positively charged protonic defects are formed. In Kröger-Vink notation this reaction is given as 11 Oxygen vacancy concentration in BZO can be increased by various techniques like strain engineering, doping, etc. Doping of BZO with lower valence cations such as Y can essentially improve the proton conductivity. [12][13][14] First principles Density Functional Theory (DFT) calculations can contribute to the understanding of the promising properties of BZO. Recently, the stability of bulk BZO and of Ba, Zr, and O vacancies at ground state has been extensively studied by DFT methods. 15 Since solid-state proton conductors such as BZO operate at relatively high temperatures, it is also very important to investigate the thermodynamic stability of the differently charged vacancy defects in dependence on the chemical environment. Sundell et al. 16 used DFT to study the neutral and fully charged Ba, Zr and O vacancies in the ground state and determined the temperature dependence of the formation of neutral and doubly positively charged O by thermodynamic modelling. In a previous DFT work we found that at 1000 K phonon contributions increase the migration barrier of the doubly positively charged oxygen vacancy by 0.35 eV. 17 Bjørheim et al. investigated the neutral and the doubly positively charged state of the O vacancy at nite temperatures and found that near the conduction band minimum the phonon contribution decreases the free formation energy of the neutral vacancy by around 0.8 eV at 1000 K, but does not alter the free formation energy of the doubly positively charged O vacancy. 18 However, in this study only the O poor chemical environment was considered. The absorption of CO 2 at high temperature was investigated in pristine as well as Y doped BZO by Polfus et al. and the role of Ba, Zr, and O vacancies was discussed. 19 Recent experimental studies clearly indicate that considerable Ba deciency may occur at high temperatures. 20 This leads to a drastic decreases of proton conductivity 21 and material stability. 22 Cation vacancy defects inuence the site occupancy of foreign atoms and dopants. 23,24 In this context, it is very important to determine the properties of the Ba vacancy at high temperatures. Also the Zr vacancies are responsible for the electrochemical properties and the thermodynamic stability of BZO. It is therefore important to study the temperature dependence of the free formation energy of all the possible vacancies in BZO, for different charge states and various atomic reservoir conditions. To the best of our knowledge, such comprehensive investigations have not been performed yet.
In this work, a systematic rst-principles study on the thermodynamic stability of Ba, Zr, and O vacancies in BZO is performed. In the rst part, the thermodynamic stability of BZO is determined in dependence on temperature and on the chemical potentials of atomic reservoirs, which are in equilibrium with BZO, i.e. on different conditions of BZO synthesis. The free formation energy of neutral and charged Ba, Zr, and O vacancies and the corresponding charge transition levels are calculated in the second part, with particular emphasis on the inuence of Fermi level, chemical environment, and temperature. In this manner, conditions for the preferential formation of doubly positively charged O vacancies and for the limitation of the number of Ba and Zr vacancies can be identied. These conditions are important prerequisites for achieving optimum proton conductivity in BZO.

Computational details
DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP) which employs plane wave basis sets. 25,26 The pseudopotentials are based on the projector augmented wave (PAW) method with exchange and correlation effects described by the generalized gradient approximation (GGA-RPBE). 27 Brillouin zone sampling is done using the Monkhorst-Pack scheme. 28 Ground state calculations are performed in a 40 atom supercell with 4 Â 4 Â 4 k-points. A plane wave basis cutoff of 500 eV is chosen for all calculations. Atomic positions as well as cell shape and volume are allowed to relax until the total stress/pressure on the supercell becomes zero. In order to obtain accurate nal equilibrium geometries, thresholds for total energy change and total force acting on each atom are chosen as 10 À8 eV and 10 À6 eVÅ À1 , respectively. To model charged defects, the appropriate number of electrons is removed from or added to the supercell and charge compensation is performed by a homogeneous jellium background charge. For example, in the case of the doubly positively charged oxygen vacancy two electrons are removed from the overall supercell. In the VASP code an approach similar to Makov and Payne is used in order to correct for the articial interactions between charged defects. 29 Table 1 demonstrates that the above mentioned choice of VASP input data yields reasonable values for ground state structural properties of solid Ba, Zr, BaO, ZrO 2 , and BZO, and the O 2 molecule.
The frozen phonon approach and the harmonic approximation 33,34 are applied to determine the vibrational frequencies which are used to investigate the effect of temperature on the materials properties. In these calculations, ground state atomic positions are slightly displaced (by 0.015Å) in positive and negative x-, y-and z-directions to obtain the respective force derivatives for the dynamical matrix. The symmetry of the underlying supercell is used to reduce the number of atomic displacements and thus computational costs. The diagonalization of the dynamical matrix yields all vibrational frequencies u i of the supercell. All the calculations are performed for the gamma point in the space of the phonon wave vectors. At a given temperature T the phonon or vibrational contribution G vib (T) to the Gibbs free energy of a system of N atoms is dened by where ħ and k B are the Planck's and Boltzmann's constants, respectively.

Gibbs free energy of formation of BZO, BaO and ZrO 2
At ground state the energy of formation of a compound is dened as the difference between the energy of the compound and the sum of energies of the elemental constituents. In the case of BZO, BaO and ZrO 2 the energy of formation per formula unit is given by  treatment of exchange and correlation effects by the GGA. These authors also found that this issue leads to a difference of about 1.36 eV (per O 2 molecule) between theoretical and experimental data for the oxidation energy of many metal oxides. Therefore, in the present work this value is added to the DFT value for E O 2 which corresponds to the procedure proposed by Wang et al. 39 This correction was also used by Bjørheim et al. in the calculation of the free formation energy of the oxygen vacancy. 18 The data obtained for the formation energy of BZO, BaO and ZrO 2 are À18.46 eV, À4.98 eV and À11.03 eV, respectively. These data are slightly different to previously reported ground state DFT results of À16.73 eV (ref. 16 The vibrational contribution to the Gibbs free energy of formation per formula unit is given by where G vib are the respective vibrational free energies at temperature T, see eqn (2). The total Gibbs free energy of formation of BZO (D f G BaZrO 3 ), BaO (D f G BaO ), and ZrO 2 (D f G ZrO 2 ), is then determined by the sum of the corresponding quantities in eqn (3)- (8). The temperature dependence of the (total) Gibbs free energy of formation is presented in Fig. 1 along with thermochemical reference data. [40][41][42] For all three compounds the absolute value of the Gibbs free energy of formation decreases if the temperature increases. There is a very good agreement between calculated and experimental data at lower temperatures. The discrepancy arising at high temperatures may be explained by anharmonic effects that are not taken into account in the calculation. The contribution of zero-point vibrations to the Gibbs free energy of formation is À0.16 eV, À0.67 eV and À0.53 eV for BZO, BaO and ZrO 2 respectively. Therefore, at 0 K the Gibbs free energy of formation is À18.62 eV (BZO), À5.65 eV (BaO), and À11.56 eV (ZrO 2 ), which is in good agreement with the thermochemical reference data of À18.27 eV (ref. 40), À5.68 eV (ref. 41), and À11.41 eV (ref. 42) for BZO, BaO and ZrO 2 , respectively. This also demonstrates that the correction of the ground state energy E O 2 (see above) leads to consistent results.
The Gibbs free energy of formation of BZO (per formula unit) from BaO and ZrO 2 is dened similarly to eqn (3)-(8) Fig. 2 shows that the absolute value of this quantity decreases with temperature, i.e. the probability of dissociation of BZO into BaO and ZrO 2 increases. This trend is in agreement with experimental results. 43 The quantitative difference at higher temperature may be due to anharmonic effects, which are not considered in the present work. The difference between the ground state value of D f G BaZrO 3 and that for 300 K is À1.34 eV, which is in good agreement with the experimental values of À1.33 eV and À1.19 eV at 300 K. 43,44

Chemical stability of BZO
The chemical stability of BZO can be determined from the constraints that give valid chemical potentials of atomic reservoirs in equilibrium with BZO. The chemical potential of an atom species i is dened as Here, the quantities m i 0 are the chemical potentials of Ba, Zr, and O atoms in their standard solid (Ba, Zr) and gaseous (oxygen) reference states. The stability of BZO is investigated by imposing constraints on Dm i . For example, the set of constraints that do not allow BZO to decompose into elements Ba, Zr and O 2 is given by Further, constraints that do not allow BZO to decompose into binary oxides BaO and ZrO 2 are governed by These constraints together with the following equation allow us to compute the stability diagram, where the suitable chemical environment for the formation of BZO can be determined. As the temperature increases, the region of chemical stability of BZO becomes smaller. This is due to the temperature dependence of the Gibbs free energies of formation D f G BaZrO 3 , D f G BaO , and D f G ZrO 2 . At 1000 K the stability region of BZO is limited by the points A 0 , B 0 , C 0 , 0, and D 0 . At points D 0 and C 0 the chemical potential of oxygen is À4.83 eV and À4.98 eV, respectively. This drastic change compared to 0 K is due to the high entropy of the gaseous O 2 molecule, which leads to a relatively fast decrease of the absolute value of the corresponding free energy if temperature increases. The strong temperature dependence of Dm O is also the cause for the loss of the stability of BZO as well as BaO and ZrO 2 by the release of oxygen, which occurs in regions that correspond to strong oxidation conditions at 0 K. Point X 0 in the intermediate stability region is obtained by shiing point X towards point 0, by a distance corresponding to the difference between A and A 0 or B and B 0 . On the other hand, the boundaries between the stability regions of BZO and BaO and those of BZO and ZrO 2 change relatively weakly if the temperature is increased.
It should be mentioned that the relation between the chemical potential m i and the atomic concentration C i is given by The partial oxygen pressure p O 2 can be obtained from the chemical potential m O using where p O denotes the standard or reference pressure. The quantity m O 0 is equal to the half of the free energy of the O 2 molecule at standard pressure (see e.g. ref. 18). In this manner, the stability of BZO may be also represented in dependence on concentration and partial pressure.

Gibbs free formation energy of Ba, Zr and O vacancies in BZO
In solid compounds like BZO, BaO, or ZrO 2 the vacancy formation is inuenced by the chemical potential of the individual atom species, which constitute the material. Furthermore, the vacancy may have different charge states. The ground state formation energy of a vacancy in a charge state q is dened as (see e.g. ref. 15) where E tot (V i q ) and E tot (bulk) are the total energy of the supercell with and without the vacancy, respectively. m i is the ground state chemical potential of the atom which was removed in order to form the vacancy. E VBM is the valence band maximum and E F is the Fermi level of the system. By denition E F is set to zero at VBM. The formation energy was calculated for different charge states of Ba, Zr and O vacancies in BZO, (Ba vacancy: 0, À1, À2; Zr vacancy: 0, À1, À2, À3, À4; O vacancy: 0, +1, +2). Our DFT calculations yield a band gap value of 3.20 eV. It is well known that the GGA used in DFT underestimates the band gap. 45 Therefore, in this work the experimental value of 5.33 eV is used in eqn (17). 46 It should be mentioned that the so-called over binding problem (oxygen energy see third section of this work) by the RPBE exchange-correlation functional as well as the underestimation of the band gap by this procedure might be overcome using the Heyd-Scuseria-Ernzerhof (HSE) functional 47-49 as well as other sophisticated approaches. 50 However, this would increase the computational costs and was therefore not considered in this work.
At temperature T the Gibbs free formation energy of a vacancy is dened similarly to eqn (17) Where G tot (V i q ,T) and G tot (bulk,T) are the total free energies of the system with and without the vacancy, respectively. These quantities correspond to the sum of the ground state energy and the respective vibrational contribution as given by eqn (2). The value of m i (T) is determined according to eqn (10). It should be mentioned the temperature dependence of the difference G tot (V i q ,T) À G tot (bulk,T) was also calculated by Bjørheim et al.
for the neutral and the doubly positively charged O vacancy. 18 Our results agree well with those of ref. 18 and 45. The solid lines in Fig. 4 depict the free formation energies of vacancies at 0 K for reservoir conditions corresponding to points A, B, C, D and X of the stability diagram (Fig. 3), whereas the dashed lines are the data for 1000 K and for points A 0 , B 0 , C 0 , D 0 , and X 0 . In dependence on the Fermi energy neutral or charged vacancies are energetically favored. In the case of points A and A 0 the charge transitions are illustrated in detail in Fig. 4b. Because of eqn (17), the Fermi level at which the charge transitions occur is independent of reservoir conditions, but the transition level depends on temperature. It should be noticed that for a given Fermi energy the diagrams show only values for the lowest free formation energy among all the possible charge states. Close to the VBM the O vacancy has the lowest free formation energy while the free formation energy of Ba and Zr vacancies is highest. The opposite behavior is found near the conduction band minimum (CBM). The values of the free formation energy strongly depend on the reservoir conditions. At points C and C 0 as well as at D and D 0 (oxygen poor conditions) the most stable defect is the O vacancy in the +2 charge state if the Fermi level is closer to VBM than to CBM. For other values of E F the negatively charged Ba and Zr vacancies are more stable. Under O rich conditions (points A and A 0 , B and B 0 ) Ba and Zr vacancies are more stable than the O vacancy, for all possible values of the Fermi energy. Near VBM these vacancies are preferentially neutral, towards the CBM they are negatively charged. In the case of the Ba and the O vacancy the difference between the data obtained for 0 and 1000 K is dependent on the value of E F , in contrast to the Zr vacancy where such a dependence is not signicant while the absolute value of the difference is higher.
For O poor/Zr rich conditions, i.e. at point C, the ground state formation energy obtained from our calculation is 5. respectively. Again, the differences are mainly due to the O 2 energy correction applied in this work (see above). Fig. 5 summarizes the data of the free formation energy of Ba, Zr, and O vacancies at 0 and 1000 K for the case that the Fermi level is at the VBM. The stability of the Ba vacancy is highest at points B and B 0 i.e. at O rich conditions, and the neutral charge state is most favored. On the other hand, at points D and D 0 , i.e. at Zr rich conditions, the free formation energy of the Ba vacancy is highest. The Zr vacancy shows a behavior similar to that of the Ba vacancy but has its minimum free formation energy at points A and A 0 . The stability of the Zr vacancy increases signicantly with temperature: while at 0 K the most probable state (points A and A 0 ) has a positive free formation energy (1.00 eV) this quantity is negative at 1000 K (À0.88 eV). Neither the Ba nor the O vacancy shows such a strong T dependence since in these cases the vibrational contributions to the free formation energy are lower. At and near the VBM the oxygen vacancy is most stable at points C and C 0 as well as at points D and D 0 whereby the doubly positively charged state is preferred. This condition favors the production of protonic defects according to reaction (1). The results depicted in Fig. 4 and 5 clearly show that the reservoir conditions at points C (C 0 ) and D (D 0 ) favor a high concentration of doubly positively charged O vacancies and a low concentration of Ba and Zr vacancies. This is an optimum prerequisite in order to achieve high proton conductivity in BZO while avoiding much Ba and Zr deciencies. As shown in a previous work 17 phonon contributions increase the migration barrier of the doubly positively charged O vacancy. This contributes to the high-temperature stability of BZO observed in experiments.  As already mentioned above, the charge state of Ba, Zr, and O vacancies changes with increasing Fermi energy. A diagram representing all possible charge transitions of Ba, Zr and O vacancies calculated at 0 K and 1000 K is given in Fig. 6. As the Fermi level increases the Ba vacancy undergoes a charge transition from 0 to À1 to À2, whereas the Zr vacancy exhibits a transition from 0 to À1 to À2 to À3 to À4, and the O vacancy undergoes transitions from +2 to +1 to 0. At 1000 K the charge transitions occur at lower values of E F than at 0 K. The inuence of temperature on the Fermi level of the charge transition is highest for the Zr vacancy. In this case at 1000 K the range of charge state À1 is extended compared to 0 K whereas the charge state À2 is not stable so that a direct transition between À1 and À3 takes place. For O vacancies, transition levels were found closer to CBM than for the other vacancies. At 0 K the (0/+) and the (+/2+) transitions occur at a distance of 0.46 and 1.43 eV from CBM, respectively. On the other hand, at 1000 K the rst and the second transition takes place at 1.15 eV and 1.77 eV, respectively, i.e. the range where the doubly positively charged O vacancy is stable is slightly smaller than at 0 K.

Conclusions
A systematic DFT study on the thermodynamic stability of BZO and of Ba, Zr, and O vacancies in this material was performed. In contrast to previous investigations, which were focused on the ground state, in the present work the stability diagram of BZO was calculated for a wide range of temperatures and for different atomic reservoir conditions. With increasing temperature, the stability region shrinks. This is mainly due to the strong temperature dependence of the chemical potential of gaseous oxygen. The temperature dependence of the free formation energy of Ba, Zr, and O vacancies in BZO was determined for all possible charge states and for the different chemical environments. To the best of our knowledge, such a comprehensive analysis has not been done yet. The free formation energy of Zr vacancies depends strongly on temperature whereas this dependence is weaker for Ba and O vacancies. Accordingly, the charge transition levels are differently inuenced by temperature. While at 0 K transitions between neutral, À1, À2, À3, and À4 states are found for the Zr vacancy, the À2 state disappears at 1000 K so that a direct transition between À1 and À3 takes place. Compared to 0 K at 1000 K the charge transition levels of the O vacancy are shied towards the valence band maximum. Present rst-principle investigations also indicate how an optimum proton conductivity in BZO can be achieved: the results clearly show that O poor reservoir conditions and a Fermi level close to the valence band maximum favor a high concentration of doubly positively charged O vacancies which is required to form many protonic defects. On the other hand, under these conditions the concentration of Ba and Zr vacancies is rather low so that Ba and Zr deciencies can be avoided.

Conflicts of interest
There are no conicts to declare.