Y doping of BaZrO3 may lead to optimum conditions for proton conduction at operating temperature of solid oxide fuel cells: a first principles study

First-principle calculations are performed to investigate Y substitutional defects at ground state and at 1000 K, for Ba- and Zr-rich chemical environments. In dependence on the Fermi level, at ground state singly positively charged Y may be potentially stable on Ba site ( YBa1+ ) and neutral as well as singly negatively charged Y on Zr site ( YZr0 and YZr1− ). However, using recent results for the doubly positively charged oxygen vacancy ( VO2+ ) and taking account charge compensation, Fermi level pinning occurs, so that under Ba-rich conditions YZr1− and VO2+ are really stable. A similar consideration yields YBa1+ and YZr1− as stable defects in the Zr-rich case. Concerning VO2+, which occurrence is a prerequisite to obtain a good proton conductor, by Y doping, at ground state only in the Ba-rich case a moderate concentration can be formed. At 1000 K the situation is improved importantly. The consideration of vibrational contributions to the free formation energy of Y on Zr site shows an increase of the stability of YZr0 and YZr1−. Under Ba-rich conditions Fermi level pinning results in a free formation energy for the doubly positively charged O vacancy of 0.481 eV which corresponds to a high VO2+ concentration and optimum conditions for proton conduction. In Zr-rich case the respective value is 0.863 eV which leads also to relatively high VO2+ occurrence but the situation is less favourable than for the Ba-rich environment.

which occurrence is a prerequisite to obtain a good proton conductor, by Y doping, at ground state only in the Ba-rich case a moderate concentration can be formed. At 1000 K the situation is improved importantly. The consideration of vibrational contributions to the free formation energy of Y on Zr site shows an increase of the stability of Y Zr 0 and -Y . Zr 1 Under Ba-rich conditions Fermi level pinning results in a free formation energy for the doubly positively charged O vacancy of 0.481 eV which corresponds to a high + V O 2 concentration and optimum conditions for proton conduction. In Zr-rich case the respective value is 0.863 eV which leads also to relatively high + V O 2 occurrence but the situation isless favourable than for the Ba-rich environment.

Introduction
In solid oxide fuel cells (SOFC) electrochemical processes convert chemical to electrical energy [1][2][3]. SOFC are highly efficient in power generation, fuel adaptation and less pollutant compared to other energy generating technologies. The current challenges in the development of SOFC are to reduce the operating temperature, enhance efficiency and thermodynamic stability. The efficiency of SOFC mostly depends on the electrolyte material and the mechanism and type of ionic conduction. Zirconate based perovskite, especially BaZrO 3 (BZO), is a well-known electrolyte with high efficiency and high-temperature chemical stability [3,4]. BZO becomes proton conducting under humid conditions. Proton incorporation happens via dissociation of water from gas phase into a hydroxide ion and a proton. The hydroxide ion fills a doubly positively charged oxygen vacancy + V , O 2 and proton forms a covalent bond with an oxygen atom of the BZO lattice. Finally, two protonic defects + OH o are formed. This can be summarized in the following reaction [3,4]  conductor. Experimental and theoretical investigations led to the conclusion that Y doped BZO may be a very promising electrolyte [2][3][4][5][6][7][8][9][10][11][12][13]. However, it was found experimentally that the thermodynamic stability of Y doped BZO decreases with increasing Y concentration [14][15][16][17][18][19][20][21]. Moreover, substitution of Y on Ba-or Zr-site may vary the electronic properties by a large extent, because of the ionic radii mismatch and different electron valence state of Y compared to Ba and Zr. By using several experimental techniques, site occupation of Y and the corresponding synthesis conditions were discussed explicitly [15][16][17]. It was found that Y on Ba site increases the Ba deficiency [21,22]. Enhancing cation vacancies leads to suppression of chemical stability and Y on Ba site may act as a trap for proton, which results in the decrease of proton conductivity. Theoretical investigations were performed to understand the influence of Y substitutional defects in the ground state and explore the possibility of Y substitution on Ba and Zr sites of BZO [23][24][25][26][27][28][29]. However, temperature dependence of the properties of these defects is crucial since perovskite based SOFC are operated at elevated temperatures. Only few studies on this topic haven been performed yet [30][31][32][33][34][35]. For example, Björketun et al [30] and Løken [31] et al discussed the increase of proton mobility by the presence of dopants (amongst them Y) on Zr and Ba sites. Kang et al [32] studied the chemical stability and proton conductivity of BZO doped with a Y or other atoms at high temperatures. Bjørheim et al [33] determined the hydration entropy of BZO doped by Y, Ga, In, and Sc. Temperature dependence of formation energy of oxygen vacancies and hydration energy of protons in Zr-site Y-, Ga-, In-, Al-, Lu-and Sc-doped BZO was studied by Takahashi et al [34]. Using first principle calculations and Monte Carlo simulation, Eisele et al [35] determined the interaction free energy of defects in Y doped BZO and pointed out the importance of ionic interactions over the hydration process at elevated temperatures.
The present work is aimed at finding high-temperature properties of Y substitution on Ba and Zr sites and their influence on the concentration of the doubly positively charged oxygen vacancy + V . O 2 In this manner, it is possible to reveal optimum conditions for proton conduction in BZO (see above). At first, using density functional theory (DFT), ground-state structure and energetics as well as electronic properties of substitutional Y defects are determined by considering different chemical environments, i.e. Ba-or Zr-rich conditions, and different charge states of the defects. These results are compared with the corresponding data for + V O 2 [36] so that the concentration of all defects in the ground state can be estimated. In the following part, the dynamical stability of these substitutional defects is studied using phonon calculations and the determination of vibrational contributions to the free energy. Then, the free formation energy is calculated and the relationship to that of + V O 2 [36] is discussed leading to results for defect concentration at elevated temperatures. Finally, the results are summarized in a conclusion section.

Computational details
DFT calculations were performed to investigate the properties of BZO with Y on Ba or Zr site. For this purpose the Vienna Ab Initio Simulation Package (VASP) [37][38][39] was used which employs plane wave basis sets, pseudopotentials obtained by the projector augmented wave (PAW) method [39,40], and the generalized gradient approximation (GGA-RPBE) [41] to treat exchange and correlation effects. For electronic structure calculations a supercell consisting of 3 × 3 × 3 BZO unit cells (135 atoms) was considered. The plane wave energy cut-off was set to 500 eV and the Brillouin zone sampling was done using the Monkhorst-pack scheme [42] with 4 × 4 × 4 k-points. The positions of the atoms as well as the volume and shape of the supercell were allowed to relax until convergence criteria are reached. In order to obtain high accuracy thresholds for total energy change and total force acting on each atom were chosen 10 −6 eV and 10 −4 eV Å −1 , respectively. Replacements of divalent Ba or tetravalent Zr by trivalent Y changes the chemical formula that describes the composition of the supercell, from BaZrO 3 to Ba 0.963 Y 0.037 ZrO 3 or BaZr 0.963 Y 0.037 O 3 , respectively. For convenience, hereafter Y on Ba site and Y on Zr site will be denoted as Y Ba and Y , Zr respectively. In order to model charged defects such as / Y Ba Zr q (Y Ba/Zr with a charge state q), the appropriate number of electrons is removed or added, and charge compensation is performed by a homogeneous jellium background charge. In the VASP code an approach similar to Makov and Payne is used in order to correct for the artificial interactions between charged defect [43,44] The temperature dependence of the Y site occupancy is investigated by considering the vibrational contributions to the free energy of Y Ba and Y Zr substitution (or the free formation energy). The vibrational frequencies were determined using the frozen phonon approach and the harmonic approximation. In these calculations the atoms were slightly displaced (by 0.015 Å) from their ground state positions in x-, y-and z-directions in order to obtain the force derivatives. Diagonalization of the dynamical matrix, which contains the force derivatives, provides the vibrational frequencies w ( ) i of the supercell. The calculations were performed for the gamma point in the space of the phonon wave vectors. The vibrational or phonon contribution to the Gibbs free energy of a supercell with N atoms is defined by Where  and k B denotes Planck's and Boltzmann's constants respectively, and T is the temperature. The free formation (or substitution) energy of the defects Y Ba q and Y Zr q is defined similarly to equation (2), but with the vibrational contributions (see equation (3)) to the corresponding total energies and chemical potentials Furthermore, the phonopy code [46] was used to determine the phonon dispersion relations along the high symmetrical Brillouin zone points.  [47,53]. Bader analysis reveals that in bulk BZO the Zr-O bond is less ionic (or more covalent) than the Ba-O bond. The replacement of a Ba or Zr atom by a Y atom at different lattice sites may alter bonding coordination and valence states. Therefore, to understand the properties of the Y substitutional defects in BZO, at first we analyzed the ground state properties of Y Ba 0 and Y Zr 0 configurations. Since a 3 × 3 × 3 supercell is used in the present work, the Y concentration on Ba or Zr sites is 3.7%. The relaxed supercell with Y Ba 0 shows a slightly smaller volume (change: −0.129 Å 3 ) than in the case of perfect BZO, but the cubic symmetry is preserved. Y Ba 0 causes a significant local distortion around the substitutional site. In particular, the nearest neighbor oxygen atoms relax towards the Y atom. The Y-O bond length is 2.888 Å, which is smaller than the Ba-O bond length in perfect BZO (3.019 Å). The bond length decrease is due to the smaller ionic radius of Y compared to that of Ba (Y 3+ : 0.9 Å, Ba 2+ : 1.49 Å), but also due to the stronger interaction between the Y atom and the nearest neighbor O atoms. As consequence, the nearest neighbor O-O bond length decreases by 0.141 Å and the central angle of Zr-O-Zr triple near the substitutional Y atom changes to 173.6°, compared to 180.0°in perfect BZO. This causes stress on the substitutional site. From Bader charge analysis [51,52], effective charge of Y and of neighboring O is found to be 74.86% and 68.50% of its nominal ionic charge (3+ and 2−), respectively. On the other hand, the effective ionic charge of Zr is equal to that in bulk BZO. Obviously, the nominally higher positive charge (3+) at the substitutional site is compensated by the negative charge of the neighboring oxygen. The strong charge redistribution between Y and O ions is an indication for the partially covalent character of this bond. This should be also the cause for the significant local distortion around the Y Ba 0 site. The relaxed Y Zr 0 configuration also preserves the cubic symmetry, and the volume of the supercell is by 0.625 Å 3 larger than that for perfect BZO. The Y-O bond length is found to be 2.242 Å, which is higher than the Zr-O bond length in bulk BZO (2.134 Å). Due to the ionic radii difference (Y 3+ : 0.9 Å, Zr 4+ : 0.72 Å), the nearest neighbor O atoms of Y move slightly outwards. The Zr-Y distance also increases to 4 a) and (b) illustrate the differential electron density of bulk BZO and Y Ba 0 configuration in the BaO plane, respectively. Similarly figures 1(c) and (d) represents the bulk BZO and Y Zr 0 configuration in the ZrO 2 plane, respectively. In the BaO plane of bulk BZO, there is a relatively constant but low electron density between Ba and O ions, which indicates that Ba-O bonding is ionic in nature (see figure 1(a)). In contrast, the Y Ba 0 configuration shown in figure 1(b) exhibits high electron density variations. The directional electron charge density variation can be attributed to the fact that Y possesses one electron more than Ba, and to the covalent character of the bond. The significant charge redistribution causes the above-mentioned local distortion near the substitutional atom and results in a decrease of the supercell volume. As already mentioned above, the Zr-O bond is partially covalent, see also figure 1(c) as well as figure S1 in the Supplementary Information, which shows a strong hybridization between Zr-d states and O-p-states. From the differential charge density plot for Y Zr 0 configuration ( figure 1(d)) one can understand that there is a high electron density close to Y, and a low electron density in the wider environment relative to that of bulk BZO. As a consequence, the Y-O bond length increases slightly in comparison with the Zr-O bond length in bulk BZO.

Results and discussion
where for the chemical potential of O the value given above is used and DH Y O 2 3 is the enthalpy (or energy) of formation of Y 2 O 3  [56]. Equation (7) yields m D Y = −5.44 eV as the upper bound, which is always used in the following considerations.
The equilibrium concentration of a particular defect i with charge q is determined by the formation energy via a Boltzmann distribution Where N and E f denote the number of possible sites available for the defect and the defect formation energy, respectively. Due to the exponential dependence in equation (9) only defects with low formation energy are of real relevance. We have considered the 0 and 1+ charge states for Y Ba q and 0, 1-and 2-charge states for Y .
Zr q Y 0 substitution for Ba is expected to be a donor since Y has one additional valence electron compared to Ba ion. Similarly, replacement of Y 0 for Zr is expected to be an acceptor. In figure 2 the ground state defect formation energies were calculated and plotted as a function of Fermi energy, for Ba-rich and Zr-rich conditions. In such metal rich regions the formation energy of O vacancy are comparatively low and therefore we consider these regions for our discussion [4]. It should be noticed that the real Fermi level in the system is determined by charge neutrality [45]. Due to equation (9) this condition is met approximately for charge compensation between the positively charged defect with the lowest formation energy and the negatively charged defect with the lowest formation energy (Fermi level pinning). Figure 2 shows the results of our defect formation energy calculation (using equation (2)) for Y q Ba and Y q Zr configurations along with the ground state defect formation energy of O vacancy data determined using data of [36]. It must be noticed that our DFT calculations yield a band gap value of 3.20 eV for BZO. It is well known that the GGA used in DFT underestimates the band gap [57]. In figure 2 the experimental value of 5.33 eV is used, also to be consistent with representations in our previous publication [36]. Under Ba-rich conditions ( figure 2(a)), at valence band maximum (VBM), the ground state formation energy of + Y Ba 1 (1.023 eV) is much lower than that of Y Ba 0 (6.243 eV, not shown in the figure). At conduction band minimum (CBM) the formation energy of Y Ba 0 is also higher than that of + Y Ba 1 (5.697 eV). Therefore, in the Ba rich case, Y on Ba site prefers the +1 charge state for all values of the Fermi energy. In the case of Y on Zr site, the formation energies of Y ,  Figure 2(a) shows that Y on Zr site prefers either the neutral or the −1 charge state in the region of the possible Fermi levels. Equilibrium charge transition from neutral to −1 state occurs at 0.091 eV above the VBM, indicates that the formation of -Y Zr 1 is dominant throughout the Fermi levels considered. The intersection of the lines for -Y Zr 1 and + V O 2 at a formation energy of 1.615 eV leads to Fermi level pinning at approximately 1.126 eV above the VBM. At this point, both defects may exist at moderate concentration. On the other hand, + Y Ba 1 is hardly formed because at this Fermi level + Y Ba 1 has relatively high formation energy (2.062 eV). Another charge compensation occurs between between + Y Ba 1 and -Y Zr 1 at 0.765 eV above the VBM. But this is of minor importance for a possible Fermi level pinning because the corresponding formation energy is 1.947 eV which is higher than that at the intersection of the lines for -Y Zr 1 and + V O 2 (see above). The obtained results mean that at Ba-rich conditions Y will preferentially form at acceptor site ( -Y Zr 1 ), and the formation of some oxygen vacancies ( is suppressed. In the Zr-rich case ( figure 2(b)), at VBM the formation energy of Y Ba 0 and + Y Ba 1 is 5.743 eV and −1.142 eV, respectively. Compared to figure 2(a) this means a reduction by 0.954 and 2.165 eV respectively. At CBM, the values for Y Ba 0 and + Y Ba 1 are 5.743 and 4.073 eV respectively, which indicates that no equilibrium charge transition occurs between them and Y Ba prefers to be in the +1 charge state. On the other hand, for Y , Therefore, in practice the intersection of the lines for -Y Zr 1 and + V O 2 is of minor relevance.

Phonon modes in BZO without and with Y Ba q and Y Zr q
Phonon modes of bulk BZO may depend on the supercell size as well as the choice of the pseudopotential [49]. Using VASP we calculated Γ-point phonon modes for 3 × 3 × 3 supercells. The use of a smaller cell size of 2 × 2 × 2 does not lead to a full set of real-number frequencies, due to the inadequate representation of forces [50,58]. In the case of a 3 × 3 × 3 supercell, no imaginary frequencies were found for bulk BZO, and figure 3 shows that their values range from 776 cm −1 to 45 cm −1 . This is consistent with the reported experimental data [59].
In the case of Y Ba 0 two imaginary phonon frequencies were found which indicates that the system is not stable. This is certainly due to the mismatch between the ionic size of Y and Ba, which creates a huge strain on the substitutional defect, leading to a dynamical instability. + Y Ba 1 substitution is also dynamically unstable with three imaginary modes. Additionally to the ionic radius mismatch, also the charge of the defect may be the cause of this instability. For Y Zr 0 substitution ( figure 3(a)), all frequencies are real numbers so that the structure is dynamically stable. The frequencies vary from 779 cm −1 to 39 cm −1 . An additional mode emerges at 523 cm −1 that corresponds to change in the Zr-O bonding due to the presence of Y atom. In addition, there is a red shift of the peaks indicating that the Y atom slightly decreases the material stability. Also in the case of -Y Zr 1 substitution ( figure 3(b)) no imaginary frequencies were found and the frequencies range from 776 cm −1 to 28 cm −1 . The new peak around 509 cm −1 is due to the modification of the Zr-O bonding by presence of the substitutional defect. Similarly tothe neutral case a slight red shift is observed as well.
Further, the phonon dispersion curves for bulk BZO and for BZO with Y Ba and Y Zr substitution were computed by the phonopy code [46]. The results for bulk BZO, Y , Ba 0 and + Y Ba 1 are presented and discussed in the Supplementary Information. Figure 4(a) shows the data for Y Zr 0 that are positive throughout the high symmetry points of the Brillouin zone, which confirms the material stability. In comparison with bulk BZO (see Supplementary Information, figure S2), noticeable decreases in the phonon frequencies were attained at the R point of the Brillouin zone, which corresponds to the ZrO 6 and YO 6 octahedra. Even though there were no negative frequencies found, softening of normal modes at R and M points indicates the decrease in the stability of Y .  4(b)). As like in Y Zr 0 a noticeable decrease in the phonon frequencies was attained at the R point of the Brillouin zone, in comparison with bulk BZO. Softening of normal modes at R and M points is somewhat higher than in the Y Zr 0 case.

Vibrational free energy and high temperature stability of Y Zr q
Since Y Ba 0 and -Y Ba 1 are not dynamically stable, in the following only the case of Y Zr q is studied. Figure 5 shows the quantity D Where DG vib Y O 2 3 is the phonon contribution to the free enthalpy (or energy) of formation given by In order to obtain the quantities G ,

Conclusions
Ground state DFT calculations on Y substitutional defects in BZO revealed that the formation of + Y , Ba 1 Y Zr 0 and -Y Zr 1 may occur. The dependence on the Fermi level was investigated in detail for two different the atomic reservoirs (Ba-rich or Zr-rich conditions). In the Ba-rich case the condition of charge neutrality between + V O 2 and -Y Zr 1 leads to Fermi level pinning at 1.126 eV and both defects form at moderate concentration. For a Zr-rich chemical environment self-compensation between + Y Ba 1 and -Y Zr 1 occurs at a Fermi level of 2.043 eV and results in a strong suppression of oxygen vacancy formation. The calculation of the vibrational modes for pure BZO and systems with Y substitutional defect on Ba and Zr sites leads to the conclusion that Y on Ba site is dynamically unstable due to the local strain and the interaction with neighbouring oxygen. Therefore, from the thermodynamic point of view this defect does not exist. This result indicates that there is no possibility for the Y Ba q to form alone at high temperature and it may form along with Ba deficiency as experimental results indicate [21,22]. Y on Zr site is dynamically stable. In particular the case of 1000 K was investigated. The results show that lattice vibrations increase the defects stability compared to the ground state. Under Ba-rich conditions charge neutrality between + V O 2 and -Y Zr 1 occurs at a Fermi level of 0.648 eV above the VBM and a free formation energy of 0.481 eV. The latter value means a significant reduction compared to the case of Fermi level pinning in the ground state where a formation energy of 1.615 eV was found. Therefore, in the Ba-rich case at 1000 K a rather high concentration of + V O 2 defects exists which means optimum conditions for proton conduction. On the other hand, in the Zr-rich case at 1000 K the charge neutrality between + V O 2 and -Y Zr 1 is found at about 0.848 eV above VBM. The corresponding free formation energy is 0.863 eV. That means that the concentration of + V O 2 islower than under Ba-rich conditions. However, also the Zr-rich case offers relatively good conditions for proton conductions.