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

Abstract

The temperature dependence of the stability of bulk BaZrO₃ (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₃ (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–8 It shows good chemical and thermodynamic stability even at very high temperature.^9,10 BZO becomes proton conducting after exposure to humid conditions.^10,11 The water from the gas phase dissociates into a hydroxide ion and a proton, the hydroxide ion fills 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¹¹

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–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.¹⁵ 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.¹⁶ 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.¹⁷ Bjørheim et al. investigated the neutral and the doubly positively charged state of the O vacancy at finite 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.¹⁸ However, in this study only the O poor chemical environment was considered. The absorption of CO₂ 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.¹⁹ Recent experimental studies clearly indicate that considerable Ba deficiency may occur at high temperatures.²⁰ This leads to a drastic decreases of proton conductivity²¹ and material stability.²² Cation vacancy defects influence 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 first-principles study on the thermodynamic stability of Ba, Zr, and O vacancies in BZO is performed. In the first 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 influence 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 identified. 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).²⁷ Brillouin zone sampling is done using the Monkhorst–Pack scheme.²⁸ 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 final equilibrium geometries, thresholds for total energy change and total force acting on each atom are chosen as 10⁻⁸ eV and 10⁻⁶ eV Å⁻¹, 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 artificial interactions between charged defects.²⁹

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₂, and BZO, and the O₂ molecule.

Table 1 Calculated values of the ground state structural properties of Ba, Zr, O, BaO, ZrO₂ and BZO, compared with theoretical and experimental data from literature. a: lattice constant, c/a: axial ratio and r₀: bond length

	This work	Theory	Experiment
Ba (BCC)
a (Å)	5.039	4.770 (ref. 30)	5.020 (ref. 31)
Zr (HCP)
a (Å)	3.220	3.229 (ref. 32)	3.232 (ref. 27)
c/a	1.732	1.600 (ref. 25)	1.592 (ref. 27)
O (dimer)
r0 (Å)	1.215	1.22 (ref. 25)	1.21 (ref. 26)
BaO (cesium chloride)
a (Å)	5.515	5.576 (ref. 15)	5.539 (ref. 28)
ZrO2 (monoclinic)
a (Å)	5.194	5.217 (ref. 15)	5.151 (ref. 29)
c/a	1.039	1.034 (ref. 15)	1.032 (ref. 29)
BaZrO3 (cubic)
a (Å)	4.269	4.208 (ref. 15)	4.192 (ref. 30)

---

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 ω_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 defined by

where ħ and k_B are the Planck's and Boltzmann's constants, respectively.

Gibbs free energy of formation of BZO, BaO and ZrO₂

At ground state the energy of formation of a compound is defined 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₂ the energy of formation per formula unit is given by:^35–38

Δ_fH_ZrO2 = E_ZrO2 − E_Zrhcp − E_O2 (5)

where, E_BaZrO3, E_BaO, E_ZrO2, E_Babcc, E_Zrhcp and E_O2 are the ground state energies of BZO, BaO, ZrO₂, Ba, Zr, and the O₂ molecule, respectively.

While DFT is able to determine rather precise values for the ground state energies of BZO, BaO, ZrO₂, Ba, and Zr, this method yields poor data for the ground state energy of the O₂ molecule. Wang et al.³⁹ attributed this to the over binding caused 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₂ 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_O2 which corresponds to the procedure proposed by Wang et al.³⁹ This correction was also used by Bjørheim et al. in the calculation of the free formation energy of the oxygen vacancy.¹⁸

The data obtained for the formation energy of BZO, BaO and ZrO₂ 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) and −17.82 eV (ref. 15) for BZO, −5.17 eV (ref. 16) and −5.55 eV (ref. 15) for BaO, and −10.25 eV (ref. 16) and −10.27 eV (ref. 15) for ZrO₂, since in these papers the O₂ correction was not taken into account.

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 (Δ_fG_BaZrO3), BaO (Δ_fG_BaO), and ZrO₂ (Δ_fG_ZrO2), 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–42

Fig. 1 Gibbs free energy of formation of BaO, ZrO₂ and BZO per formula unit, in comparison with reported data obtained from experiments and thermodynamic modelling.^40–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₂ 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₂), 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₂, respectively. This also demonstrates that the correction of the ground state energy E_O2 (see above) leads to consistent results.

The Gibbs free energy of formation of BZO (per formula unit) from BaO and ZrO₂ is defined 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₂ increases. This trend is in agreement with experimental results.⁴³ 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 Δ_fG_BaZrO3 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

Fig. 2 Gibbs free energy of formation of BZO (per formula unit) from BaO and ZrO₂ as a function of temperature.

Fig. 1 and 2 clearly demonstrate the importance of vibrational effects on the thermodynamic stability of BZO, BaO and ZrO₂.

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 defined as

μ_i = μ_i⁰ + Δμ_i (10)

Here, the quantities μ_i⁰ 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 Δμ_i. For example, the set of constraints that do not allow BZO to decompose into elements Ba, Zr and O₂ is given by

Δμ_Ba ≤ 0, Δμ_Zr ≤ 0 and Δμ_O ≤ 0 (11)

Further, constraints that do not allow BZO to decompose into binary oxides BaO and ZrO₂ are governed by

Δ_fG_BaO ≥ Δμ_Ba + Δμ_O (12)

Δ_fG_ZrO2 ≥ Δμ_Zr + 2Δμ_O (13)

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.

Δ_fG_BaZrO3 = Δμ_Ba + Δμ_Zr + 3Δμ_O (14)

Fig. 3 represents the stability diagram of BZO with the competing binary oxides BaO and ZrO₂. Points A, B, C, D and X represent five different chemical potential conditions at 0 K. These conditions are obtained from eqn (10)–(14) and using the corresponding values for the Gibbs free energy of formation at 0 K. Points A and B correspond to O rich conditions with the chemical potential value Δμ_O = 0 eV. At point C the environment is Zr rich, at point D, Ba rich. The condition of extreme reduction can be obtained at point C, where the chemical potential of oxygen is most negative: Δμ_O = −5.78 eV. Point X is in the intermediate stability region of BZO between oxidation and reduction conditions.

Fig. 3 The area of chemical stability of BZO for different temperature, from 0 K to 1500 K. At 0 K the points A, B, C, 0 and D are the corners of the stability region, whereas at 1000 K BZO is stable within the area limited by A′, B′, C′, 0, and D′. The first, second, and third coordinates denote the values of Δμ_i for Ba, Zr, and O, respectively.

Fig. 3 can be compared with the results obtained by Abbas et al.¹⁵ They calculated the stability diagram for the ground state using all electron calculation (WIEN2k code) and did not include the O₂ correction (see above). Comparison of our ground state chemical potential values of Ba, Zr, and O with the literature data may helpful to understand the significance of O₂ correction. For example, at point C, O poor/Zr rich conditions, the atomic chemical potential values for the ground state obtained from our calculation is −1.96 (Δμ_Ba), 0.00 (Δμ_Zr) and −5.50 (Δμ_O), which is slightly different from the literature values of −2.423 (ref. 15), 0.000 (ref. 15) and −5.134 (ref. 15) for Δμ_Ba, Δμ_Zr, and Δμ_O respectively. The deviations are mainly due to the inclusion of the O₂ correction, but certainly also due to the use of a pseudopotential in the VASP calculations. At other points of the stability diagram similar differences between our ground state values of the chemical potential and the corresponding data of Abbas et al.¹⁵ occur.

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 Δ_fG_BaZrO3, Δ_fG_BaO, and Δ_fG_ZrO2. At 1000 K the stability region of BZO is limited by the points A′, B′, C′, 0, and D′. At points D′ and C′ 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₂ 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 Δμ_O is also the cause for the loss of the stability of BZO as well as BaO and ZrO₂ by the release of oxygen, which occurs in regions that correspond to strong oxidation conditions at 0 K. Point X′ in the intermediate stability region is obtained by shifting point X towards point 0, by a distance corresponding to the difference between A and A′ or B and B′. On the other hand, the boundaries between the stability regions of BZO and BaO and those of BZO and ZrO₂ change relatively weakly if the temperature is increased.

It should be mentioned that the relation between the chemical potential μ_i and the atomic concentration C_i is given by

The partial oxygen pressure p_O2 can be obtained from the chemical potential μ_O using

where p^O denotes the standard or reference pressure. The quantity μ_O⁰ is equal to the half of the free energy of the O₂ 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₂ the vacancy formation is influenced 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 defined as (see e.g. ref. 15)

ΔE_f(V_i^q) = E_tot(V_i^q) − E_tot(bulk) + q(E_VBM + E_f) + μ_i (17)

where E_tot(V_i^q) and E_tot(bulk) are the total energy of the supercell with and without the vacancy, respectively. μ_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 definition 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.⁴⁵ Therefore, in this work the experimental value of 5.33 eV is used in eqn (17).⁴⁶ 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.⁵⁰ 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 defined similarly to eqn (17)

ΔG_f(V_i^q,T) = G_tot(V_i^q,T) − G_tot(bulk,T) + q(E_VBM + E_F) + μ_i(T) (18)

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 μ_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.¹⁸ 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′, B′, C′, D′, and X′. In dependence on the Fermi energy neutral or charged vacancies are energetically favored. In the case of points A and A′ 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′ as well as at D and D′ (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′, B and B′) 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 significant while the absolute value of the difference is higher.

Fig. 4 Gibbs free formation energies of Ba (green), Zr (blue) and O (red) vacancies at 0 K (solid line) and 1000 K (dashed line), as function of Fermi energy, for different chemical environments (see Fig. 3).

For O poor/Zr rich conditions, i.e. at point C, the ground state formation energy obtained from our calculation is 5.39 eV (V_Ba⁰), 14.56 (V_Zr⁰), and −3.43 (V_O²⁺), whereas Abbas et al.¹⁵ obtained is 5.66 eV, 16.30 eV and −1.84 eV for V_Ba⁰, V_Zr⁰ and V_O²⁺ respectively. Again, the differences are mainly due to the O₂ 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′ i.e. at O rich conditions, and the neutral charge state is most favored. On the other hand, at points D and D′, 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′. The stability of the Zr vacancy increases significantly with temperature: while at 0 K the most probable state (points A and A′) 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′ as well as at points D and D′ 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′) and D (D′) 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 deficiencies. As shown in a previous work¹⁷ 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.

Fig. 5 Free formation energy of Ba, Zr and O vacancies at the valence band maximum (VBM), for different points of the stability diagram (see Fig. 3). Data for 0 K and 1000 K data are shown by solid and dashed lines, respectively.

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 influence 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 first 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.

Fig. 6 Charge transitions of Ba, Zr and O vacancies at 0 K and 1000 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 influenced 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 shifted towards the valence band maximum. Present first-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 deficiencies can be avoided.

Conflicts of interest

There are no conflicts to declare.

References

1. H. Iwahara, T. Esaka, H. Uchida and N. Maeda, Solid State Ionics, 1981, 3–4, 359–363.
2. K. D. Kreuer, Chem. Mater., 1996, 8, 610–641.
3. K. D. Kreuer, Solid State Ionics, 1997, 97, 1–15.
4. H. Iwahara, Proton Conductors: Solids, Membranes and Gels-Materials and Devices, Cambridge University Press, Cambridge, England, 1992.
5. S. M. Haile, G. Staneff and K. H. Ryu, J. Mater. Sci., 2001, 36, 1149–1160.
6. D. Chuancheng, T. Jianhua, S. Meng, N. Stefan, S. Michael, R. Sandrine, A. Ali and O. H. Ryan, Science, 2015, 349, 1321–1326.
7. H. Malerød-Fjeld, D. Clark, I. Yuste-Tirados, R. Zanón, D. Catalán-Martinez, D. Beeaff, S. H. Morejudo, P. K. Vestre, T. Norby, R. Haugsrud, J. M. Serra and C. Kjølseth, Nat. Energy, 2017, 2, 923–931.
8. E. Vøllestad, R. Strandbakke, M. Tarach, D. Catalán-Martínez, M. L. Fontaine, D. Beeaff, D. R. Clark, J. M. Serra and T. Norby, Nature Materials, 2019, 18, 752–759.
9. M. E. Björketun, P. G. Sundell and G. Wahnström, Faraday Discuss., 2007, 134, 247–265.
10. W. Munch, G. Seifert, K. D. Kreuer and J. Maier, Solid State Ionics, 1997, 97, 39–44.
11. K. Kreuer, Annu. Rev. Mater. Res., 2003, 33, 333–359.
12. P. Babilo, T. Uda and S. M. Haile, J. Mater. Res. Technol., 2007, 22, 1322–1330.
13. S. Tao and J. T. S. Irvine, J. Solid State Chem., 2007, 180, 3493–3503.
14. D. Han, N. Hatada, T. Uda and R. Koc, J. Am. Chem. Soc., 2016, 99, 3745–3753.
15. S. M. Alay-e-Abbas, S. Nazir and A. Shaukat, Phys. Chem. Chem. Phys., 2016, 18, 23737–23745.
16. P. G. Sundell, M. E. Björketun and G. Wahnström, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104112.
17. R. Nadarajan, D. Murali, M. Posselt and S. V. M. Satyanarayana, Phys. Status Solidi B, 2018, 255, 1700398.
18. T. S. Bjørheim, M. Arrigoni, D. Gryaznov, E. Kotomin and J. Maier, Phys. Chem. Chem. Phys., 2015, 17, 20765–20774.
19. J. M. Polfus, B. Yildiz, H. L. Tuller and R. Bredesen, J. Phys. Chem. C, 2018, 122, 307–314.
20. S. Imashuku, T. Uda and Y. Awakura, Electrochem. Solid-State Lett., 2007, 10, B175–B178.
21. D. Han, K. Kishida, H. Inui and T. Uda, RSC Adv., 2014, 4, 31589–31593.
22. R. Sažinas, C. Bernuy-López, M.-A. Einarsrud and T. Grande, J. Am. Ceram. Soc., 2016, 99, 3685–3695.
23. D. Han, K. Kishida, K. Shinoda, H. Inui and T. Uda, J. Mater. Chem. A, 2013, 1, 3027.
24. R. Sažinas, M. F. Sunding, A. Thøgersen, I. Sakaguchi, T. Norby, T. Grande and J. M. Polfus, J. Mater. Chem. A, 2019, 7, 3848–3856.
25. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11–169.
26. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
27. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
28. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
29. G. Makov and M. C. Payne, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4014.
30. P. Erhart and K. Albe, J. Appl. Phys., 2007, 102, 084111.
31. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 8 edn, 2004.
32. Y. Nie and Y. Xie, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 174117.
33. C. Wolverton, V. Ozolin and M. Asta, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 144109.
34. S. Huang, D. L. Worthington, M. Asta, V. Ozolins, G. Ghosh and P. K. Liaw, Acta Mater., 2010, 58, 1982–1993.
35. C. G. V. d. Walle and J. Neugebauer, J. Appl. Phys., 2004, 95, 3851–3879.
36. P. Erhart and K. Albe, J. Appl. Phys., 2007, 102, 084111.
37. T. Tanaka, K. Matsunaga, Y. Ikuhara and T. Yamamoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 205213.
38. H. Raebiger, S. Lany and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 045209.
39. L. Wang, T. Maxisch and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 195107.
40. M. Chase, J. Phys. Chem. Ref. Data, 1996, 25, 551–603.
41. J. A. Dean, Lange's handbook of chemistry, McGraw-Hill, New York, London, 1999.
42. E. Cordfunke and R. Konings, Thermochim. Acta, 1989, 156, 45–51.
43. K. T. Jacob and Y. Waseda, Metall. Mater. Trans. B, 1995, 26, 775.
44. M. D. Gonçalves, P. S. Maram, R. Muccillo and A. Navrotsky, J. Mater. Chem. A, 2014, 2, 17840–17847.
45. J. M. Polfus, M.-L. Fontaine, A. Thøgersen, M. Riktor, T. Norby and R. Bredesen, J. Mater. Chem. A, 2016, 4, 8105–8112.
46. J. Robertson, J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.--Process., Meas., Phenom., 2000, 18, 1785.
47. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
48. J. Heyd and G. E. Scuseria, J. Chem. Phys., 2004, 121, 1187–1192.
49. J. Heyd and G. E. Scuseria, J. Chem. Phys., 2006, 124, 219906.
50. A. Lindman, P. Erhart and G. Wahnström, Metall. Mater. Trans. B, 2015, 91, 245114.

---