Conductivities in Yttrium-Doped Barium Zirconate: A First-Principles Study

: Yttrium-doped barium zirconate (BZY) has emerged as an attractive candidate for application as a proton (H + )-conducting solid electrolyte due to its high ionic conductivity and excellent chemical stability. In this study, the conductivities of BaZr (1 − x) Y x O 3 − δ (BZY, x = 0, 0.037, 0.074, 0.148, and 0.22) with different carriers were studied based on density functional theory (DFT) and experiments. The results revealed that yttrium doping can effectively reduce the energy barrier for the migration of protons and oxygen ions (O 2 − ). When comparing the energy barriers for protons and oxygen ions, the energy barriers for proton migration were found to be lower than those for oxygen ion migration, which indicates that a proton conductor can offer the advantages of lower activation energy and, possibly, higher conductivity. An analysis of the electronic structure of the BZYs found that the top of the valence band exceeded the Fermi energy level following yttrium doping. As a result, the electron conductivity increased as the yttrium content increased. Furthermore, this study also tested the total conductivity of BaZr (1 − x) Y x O 3 − δ (BZY, x = 0.1, 0.2, 0.3, and 0.4) and found the trend of the total conductivity to be consistent with the results of the DFT calculations.


Introduction
In recent years, solid oxide fuel cells (SOFCs) have been widely used in the conversion of chemical energy into electrical energy [1][2][3]].An SOFC is composed of three main components: the cathode, electrolyte, and anode [4,5].As a core component, the electrolyte is key to the performance of SOFCs [6].
Typically, zirconia-, ceria-, bismuth-oxide-, and perovskite-oxide-based electrolytes are used in SOFCs [7,8].The zirconia-based electrolyte is represented by 8 mol% Y 2 O 3stabilized ZrO 2 (8-YSZ), which is the earliest, most mature, and most widely used SOFC electrolyte [9].Its ionic conductivity is higher than 0.1 S/cm at 1000 • C. X.J. Chen et al. [10] reported a SOFC based on YSZ electrolyte, which showed a conductivity of 0.12 S/cm at 1000 • C.However, this conduction of ions is a thermal activation process, and high conductivity can only be obtained at high temperatures, which leads to a high cost with regard to the jointing material and hinders the long-term stability of SOFCs [5,11].A new solid proton-conducting electrolyte, such as yttrium-doped barium zirconate (BZY), exhibits lower working temperatures and higher conductivity than YSZ between 300 • C and 600 • C, meaning that it can facilitate the development of intermediate-to low-temperature fuel cell technology [12][13][14].We call it a proton ceramic fuel cell (PCFC).
As BZY is a protonic, oxygen ionic, and electronic mixed conductor, yttrium doping will directly affect these conductivities, thereby altering the performance of the fuel cell [15][16][17][18].To reveal the mechanism of BZY conductance, X.J. Li et al. [19] used reactive molecular dynamics (RMD) to simulate the O 2− transport properties and mechanisms of BZY coexisting oxygen vacancies, dopants, and edge dislocations.F.M. Draber et al. [20] combined DFT calculations and KMC techniques to simulate the macroscopic oxygen ion conductivity of BZY; the change in the average proton mobility with dopant fraction was also studied by the same method.The average proton mobility first decreases and then increases with the increase in dopant fraction [21].B. Merinov et al. [22] used density functional theory (DFT) to investigate the proton diffusion paths and rates in Y-doped BaZrO 3 .Y. Yamazali et al. [23] explored the effect of proton trapping on proton transport in BZY.Protons must overcome a certain association energy and general activation energy to achieve long-range transport, and it has been proposed that proton conductivity can be improved by determining a suitable alternative dopant to reduce the association energy.Although there are many theoretical studies to reveal the mechanism of BZY conductance, and while it is common for authors to report on some of these aspects in isolation, we consider all three aspects in complete research to achieve a greater understanding of this electrolyte material.V.P. Gorelov et al. [24] measured the ionic, proton, and oxygen conductivities as functions of air humidity (pH 2 O = 0.04 ~3.57kPa) in the BaZr (1-x) Y x O 3-δ (BZY, x = 0.02 ~0.15) over a temperature range of 600 ~900 • C.This represents a cumbersome way of revealing the mixed conductivity of BZY coexisting protons, oxygen ions, and electrons through experimental means, where the theoretical calculation is a reliable technology with which to understand the relevant mechanism.
In this work, we combine the computational and experimental methods to analyze the conductive mechanism of BaZr (1-x) Y x O 3-δ (BZY, x = 0 ~0.22) using DFT to clarify the contributions of the protons, oxygen ions, and electrons to the total conductivity under the condition of different yttrium contents.Furthermore, we also explain the dependence of the BZY conductivities on the Y content.First, we calculate the proton energy barriers for reorientation and intra-octahedral hopping in BZY with different yttrium contents.Second, we calculate the migration energy barrier for oxygen ions to investigate the effect on the oxygen ion conductivity in the different BZYs.Third, we calculate the energy band to reflect the electronic conductivities by analyzing the relationship between the valence band and Fermi energy level.Finally, we test the total conductivity of BZY by means of experiments.The influence of yttrium doping on the performance of BZY electrolytes is determined based on the experimental tests, which provide findings consistent with the calculation results.

Calculation Details
The works were implemented on the basis of density functional theory (DFT) calculations with Vienna Ab initio simulation package (VASP), and the exchange-correlation functional was described using the Perdew-Burke-Ernzerhof (PBE) parametrization within the generalized gradient approximation (GGA) [25][26][27][28].The core electrons of each atom were described using the projector-augmented wave method (PAW) [29,30].The energy cutoff was chosen to be 520 eV.The valence electrons of the pseudopotentials were 2s 2 2p 4 for O, 5s 2 5p 6 6s 2 for Ba, 4s 2 4p 6 5s 2 4d 2 for Zr, 4s 2 4p 6 5s 2 4d 1 for Y, and 1s 1 for H.All structures were relaxed using the conjugate gradient method until the forces on each atom were less than 0.02 eV/Å.The on-site Coulomb-U correction was adopted for Zr 4d 2 and Y 4d 1 electrons [31].A 3 × 3 × 3 Monkhorst-Pack k-point mesh was used to optimize the supercells of BaZrO 3 , consisting of 3 × 3 × 3 cubic unit cells (135 atoms), and treat the diffusion behavior for protons in the compound.The minimum energy path connecting the initial and final configuration of proton rotation and migration was found using the climbing image nudged elastic band (CI-NEB) method [32,33].Three images were linearly interpolated between the initial and final states.The formation energy calculation Equation ( 1) is as follows: Crystals 2023, 13, 401 3 of 14 E substrate , N i , and E i−bulk are the total energies of the substrate, the number of the atoms of element i in the BZY structure, and the total energies of the bulk atom of element i, respectively.
BZO was cubic in the space group Pm3m [34].The optimized lattice parameter = 4.27 Å, which is essentially the same as the lattice parameter obtained in previous theoretical studies (a = 4.25 Å) [35,36] and in agreement with the experimental value (a = 4.19 Å).

Fabrication and Testing of Samples
BaZr (1-x) Y x O 3-δ (BZY, x = 0 ~0.4) powders were synthesized via the sol-gel method, and the chelating agent was citric acid.According to stoichiometric ratio, the starting commercial materials Y(NO 3 ) 3 •6H 2 O (Macklin, 99.9%), BaCO 3 (Macklin, 99.95%), and Zr(NO 3 ) 4 •5H 2 O (Macklin, 99.5%) were sequentially added to deionized water.The molar ratio of metal cations to citric acid was controlled at 1: 1.5, and the transparent solution was obtained.Then, the solution was continuously heated and stirred at 100 • C, and 25% NH 3 solution was dripped into the solution until it reached neutral (pH ~7).The water was continuously evaporated during the heating process, and opaline gel was obtained.After that, the gel was heated at 280 • C in air for 12 h.Whereafter, the obtained ash was ground in an agate mortar and calcined at 1100 • C for 2 h in air.The dry pressing method was used to fabricate the bulk electrolyte, and the binder solution was 3 wt % polyvinyl alcohol and mixed with the calcined powders.The disk-shaped samples were obtained by uniaxial pressing in a cylindrical mold at a pressure of 120 MPa.Then, the electrolyte samples were sintered at 1600 • C for 12 h in air.The platinum paste (Pt-7840, Sino-Platinum Metals, Guiyang, China) was screen-printed on both sides of the sample and calcinated at 900 • C for 2 h in air to prepare electrodes.Platinum wires were used as leads with a small amount of Pt paster as binder.Schematic diagrams of samples are shown in Figure 1.
final configuration of proton rotation and migration was found using the climbing image nudged elastic band (CI-NEB) method [32,33].Three images were linearly interpolated between the initial and final states.The formation energy calculation Equation ( 1) is as follows: E substrate , Ni, and E -bulk are the total energies of the substrate, the number of the atoms of element i in the BZY structure, and the total energies of the bulk atom of element i, respectively.
BZO was cubic in the space group Pm3 _ m [34].The optimized lattice parameter = 4.27 Å, which is essentially the same as the lattice parameter obtained in previous theoretical studies (a = 4.25 Å ) [35,36] and in agreement with the experimental value (a = 4.19 Å ).

Fabrication and Testing of Samples
BaZr(1-x)YxO3-δ (BZY, x = 0 ~ 0.4) powders were synthesized via the sol-gel method, and the chelating agent was citric acid.According to stoichiometric ratio, the starting commercial materials Y(NO3)3•6H2O (Macklin, 99.9%), BaCO3 (Macklin, 99.95%), and Zr(NO3)4•5H2O (Macklin, 99.5%) were sequentially added to deionized water.The molar ratio of metal cations to citric acid was controlled at 1: 1.5, and the transparent solution was obtained.Then, the solution was continuously heated and stirred at 100 °C , and 25% NH3 solution was dripped into the solution until it reached neutral (pH ~ 7).The water was continuously evaporated during the heating process, and opaline gel was obtained.After that, the gel was heated at 280 °C in air for 12 h.Whereafter, the obtained ash was ground in an agate mortar and calcined at 1100 °C for 2 h in air.The dry pressing method was used to fabricate the bulk electrolyte, and the binder solution was 3 wt % polyvinyl alcohol and mixed with the calcined powders.The disk-shaped samples were obtained by uniaxial pressing in a cylindrical mold at a pressure of 120 MPa.Then, the electrolyte samples were sintered at 1600 °C for 12 h in air.The platinum paste (Pt-7840, Sino-Platinum Metals, Guiyang, China) was screen-printed on both sides of the sample and calcinated at 900 °C for 2 h in air to prepare electrodes.Platinum wires were used as leads with a small amount of Pt paster as binder.Schematic diagrams of samples are shown in Figure 1.X-ray diffraction (XRD) patterns of the BZY powders were analyzed on a D8 AD-VANCE X-ray diffractometer (Bruker, Billerica, MA, USA) using Cu Kα radiation.The testing system for conductivity testing is shown in Figure 2, which consisted of a tube X-ray diffraction (XRD) patterns of the BZY powders were analyzed on a D8 AD-VANCE X-ray diffractometer (Bruker, Billerica, MA, USA) using Cu Kα radiation.The testing system for conductivity testing is shown in Figure 2, which consisted of a tube furnace and mass flow controllers.The content of water vapor in the simulated gas was controlled by a peristaltic pump, and the water was vaporized with a heating coil around the stainless-steel tube.Flow rates of the gases were kept at 100 mL min −1 .The conductivities were measured using an electrochemical workstation (Interface 1010E, Gamry, Warminster, PA, USA) in the electrochemical impedance spectroscopy (EIS) model.During the tests, the frequency range was from 0.1 Hz to 1 MHz, and the perturbation voltage was 50 mV.
Testing temperature was fixed at 450 • C. The testing environment for the conductivity measurements of BZY was 3% H 2 O +air.
furnace and mass flow controllers.The content of water vapor in the simulated gas was controlled by a peristaltic pump, and the water was vaporized with a heating coil around the stainless-steel tube.Flow rates of the gases were kept at 100 mL min −1 .The conductivities were measured using an electrochemical workstation (Interface 1010E, Gamry, Warminster, PA, USA) in the electrochemical impedance spectroscopy (EIS) model.During the tests, the frequency range was from 0.1 Hz to 1 MHz, and the perturbation voltage was 50 mV.Testing temperature was fixed at 450 °C .The testing environment for the conductivity measurements of BZY was 3% H2O +air.

Supercell Structures
As shown in Figure 3a, the 3 × 3 × 3 supercell of BaZrO3 (BZO) consisted of 27 Ba atoms, 27 Zr atoms, and 81 oxygen atoms.The 27 Zr atoms constructed a 3 × 3 × 3 array, which was used for the next Y atom's replacement.Obviously, different yttrium doping contents have various crystal configurations.In this work, four models were employed with one, two, four, and six Y atoms to replace Zr atoms, corresponding to yttrium doping of 3.7%, 7.4%, 14.8%, and 22.2% in mole ratios.We refer to each of these systems as BZY3.7,BZY7.4,BZY14.8, and BZY22.2.In order to exhaust the most possible configurations of the above supercell, two situations were taken into account in the calculation models: Firstly, the types of B-site atoms were classified according to the distance between the atoms, as shown in Figure 3b.Y1 (yellow B-site atom) is fixed in the center of the supercell as the original point.Moreover, the distance between Y1 and Y2 (purple B-site atom) is a, the distance from Y3 (orange B-site atom) is √2a, and the distance from Y4 (brown B-site atom) is √3a, as shown in the amplified picture on the top left of Figure 3b.Secondly, according to the practical situations, a dispersed distribution of Y atoms in the Bravais lattice is most possible.When the number of doping atoms is more than 3, it can be ensured that B-site atoms are replaced in each (100) face.Typical configurations of BaZr(1x)YxO3-δ (BZY, x = 0, 0.037, 0.074, 0.148, and 0.22) are shown in Figure 4. BZY0 (BaZrO3) represents the BZY without Y atom replacement.BZY3.7 represents the BZY with one Y atom replacement, where the atom is set in the center of the supercell.In Figure 4c-e, the supercell with two Y atoms is illustrated.It can be seen that the distance between the two of them is a, √2a, and √3a.The next step is four Zr atoms, which are replaced by four Y atoms, as can be seen in Figure 4f-i.The total distances between the three surrounding Y atoms to the centered Y atom are 3 × a, 3 × √2a, and 3 × √3a, which cover all conditions of the four Y atoms' replacement.Furthermore, the six Y atoms' replacements are shown

Supercell Structures
As shown in Figure 3a, the 3 × 3 × 3 supercell of BaZrO 3 (BZO) consisted of 27 Ba atoms, 27 Zr atoms, and 81 oxygen atoms.The 27 Zr atoms constructed a 3 × 3 × 3 array, which was used for the next Y atom's replacement.Obviously, different yttrium doping contents have various crystal configurations.In this work, four models were employed with one, two, four, and six Y atoms to replace Zr atoms, corresponding to yttrium doping of 3.7%, 7.4%, 14.8%, and 22.2% in mole ratios.We refer to each of these systems as BZY3.7,BZY7.4,BZY14.8, and BZY22.2.In order to exhaust the most possible configurations of the above supercell, two situations were taken into account in the calculation models: Firstly, the types of B-site atoms were classified according to the distance between the atoms, as shown in Figure 3b.Y1 (yellow B-site atom) is fixed in the center of the supercell as the original point.Moreover, the distance between Y1 and Y2 (purple B-site atom) is a, the distance from Y3 (orange B-site atom) is √ 2a, and the distance from Y4 (brown B-site atom) is √ 3a, as shown in the amplified picture on the top left of Figure 3b.Secondly, according to the practical situations, a dispersed distribution of Y atoms in the Bravais lattice is most possible.When the number of doping atoms is more than 3, it can be ensured that B-site atoms are replaced in each (100) face.Typical configurations of BaZr   On the basis of these considerations, the formation energy of each BZY structure was computed according to Equation (1).The calculated formation energies of those structures are tabulated in Table 1.Obviously, the BZY7.4_3,BZY14.8_2,BZY14.8_3, and BZY22.2_4configurations have the lowest formation energy in each replacement, which are −61685.09,−61075.28,−61075.28,and −60457.74kJ•mol −1 , respectively.It is indicated that these doped configurations are the most possible among our considered ones.Among them, BZY14.8_2 and BZY14.8_3have the same formation energy, which is caused by the periodicity and symmetry of the structure, indicating that these two doped structures can be regarded as one situation.In the following work, we used the BZY14.8_2structure for   On the basis of these considerations, the formation energy of each BZY structure was computed according to Equation (1).The calculated formation energies of those structures are tabulated in Table 1.Obviously, the BZY7.4_3,BZY14.8_2,BZY14.8_3, and BZY22.2_4configurations have the lowest formation energy in each replacement, which are −61685.09,−61075.28,−61075.28,and −60457.74kJ•mol −1 , respectively.It is indicated that these doped configurations are the most possible among our considered ones.Among them, BZY14.8_2 and BZY14.8_3have the same formation energy, which is caused by the periodicity and symmetry of the structure, indicating that these two doped structures can be regarded as one situation.In the following work, we used the BZY14.8_2structure for On the basis of these considerations, the formation energy of each BZY structure was computed according to Equation (1).The calculated formation energies of those structures are tabulated in Table 1.Obviously, the BZY7.4_3,BZY14.8_2,BZY14.8_3, and BZY22.2_4configurations have the lowest formation energy in each replacement, which are −61685.09,−61075.28,−61075.28,and −60457.74kJ•mol −1 , respectively.It is indicated that these doped configurations are the most possible among our considered ones.Among them, BZY14.8_2 and BZY14.8_3have the same formation energy, which is caused by the periodicity and symmetry of the structure, indicating that these two doped structures can be regarded as one situation.In the following work, we used the BZY14.8_2structure for analysis.The relationship between the doping ratio and formation energies is shown in Figure 5. analysis.The relationship between the doping ratio and formation energies is shown in Figure 5.

Proton Migration in BZYs
Using the ideal structures obtained above, we explore the effect of yttrium content on energy barriers of proton migration and Bader charge.Because of the lack of oxygen vacancy, proton transfer can only proceed in H-form, the proton hopping between adjacent oxygen atoms at normal lattice positions is attributed to the Grotius mechanism [37,38].On the basis of the analysis in Section 3.1.1,the same migration path around the central orthorhombic cell was selected to compare the proton migration energy barriers in every configuration, and the migration path (H1-H2-H3) was selected at the central octahedron of the BZY supercells, as shown in Figure 6b-f.In detail, as shown in Figure 6a, the energy barriers of reorientation (H1-H2) and intra-octahedral hopping (H2-H3) were calculated step by step.Note that it has been reported in previous works that Ba perovskites are unfavorable for proton inter-octahedral hopping, which will not be the focus in our work [39][40][41].
in every configuration, and the migration path (H1-H2-H3) was selected at the central octahedron of the BZY supercells, as shown in Figure 6b-f.In detail, as shown in Figure 6a, the energy barriers of reorientation (H1-H2) and intra-octahedral hopping (H2-H3) were calculated step by step.Note that it has been reported in previous works that Ba perovskites are unfavorable for proton inter-octahedral hopping, which will not be the focus in our work.[39][40][41].The energy curves of BZYs are shown in Figure 8a, which were calculated by the CI-NEB method along the above paths in Figure 6b-f.The calculated results show that the energy barriers of proton rotation and migration in BZY0 are 0.18 and 0.66 eV, respectively, which are consistent with other literature [42,43].The energy barriers of proton rotation and migration in BZY3.7,BZY7.4,BZY14.8, and BZY22.2 are 0.16 and 0.22 eV, 0.21 and 0.19 eV, 0.14 and 0.05 eV, and 0.23 and 0.05 eV, respectively.Similar work by F.M. Draber et al. [21] considered all eight and all three possible permutations of the occupation of the zirconium sites in the first nearest neighborhood of the initial and final positions of the translation jumps and the reorientation jumps, respectively.When all three sites are occupied by Zr, the proton migration barrier is 0.41 eV, which is slightly lower than our calculated value of 0.66 eV.When one site and two sites are occupied by Y atoms, the value of our energy barrier is roughly consistent with that reported in the reference.From the point of thermodynamics, the higher the energy barrier, the more difficult the proton migration.Therefore, when the proton migrates from the initial state H1 to the final state H3 in BZYs, the energy barrier is actually attributed to the higher local barrier in the reorientation (H1-H2) step and the intra-octahedral hopping (H2-H3) step.Similarly, the energy barriers of proton migration in BZY3.7,BZY7.4,BZY14.8, and BZY22.2 are 0.22, 0.21, 0.14, and 0.23 eV, respectively, as shown in the inset of Figure 8a.In terms of these energy barrier values, they are much smaller than those in the case of BZY0.This difference strongly suggests that proton migration is much easier for BZY3.7 ~ BZY22.2 than for BZY0.In other words, doping can significantly improve the migration of protons.More The energy curves of BZYs are shown in Figure 7a, which were calculated by the CI-NEB method along the above paths in Figure 6b-f.The calculated results show that the energy barriers of proton rotation and migration in BZY0 are 0.18 and 0.66 eV, respectively, which are consistent with other literature [42,43].The energy barriers of proton rotation and migration in BZY3.7,BZY7.4,BZY14.8, and BZY22.2 are 0.16 and 0.22 eV, 0.21 and 0.19 eV, 0.14 and 0.05 eV, and 0.23 and 0.05 eV, respectively.Similar work by F.M. Draber et al. [21] considered all eight and all three possible permutations of the occupation of the zirconium sites in the first nearest neighborhood of the initial and final positions of the translation jumps and the reorientation jumps, respectively.When all three sites are occupied by Zr, the proton migration barrier is 0.41 eV, which is slightly lower than our calculated value of 0.66 eV.When one site and two sites are occupied by Y atoms, the value of our energy barrier is roughly consistent with that reported in the reference.From the point of thermodynamics, the higher the energy barrier, the more difficult the proton migration.Therefore, when the proton migrates from the initial state H1 to the final state H3 in BZYs, the energy barrier is actually attributed to the higher local barrier in the reorientation (H1-H2) step and the intra-octahedral hopping (H2-H3) step.Similarly, the energy barriers of proton migration in BZY3.7,BZY7.4,BZY14.8, and BZY22.2 are 0.22, 0.21, 0.14, and 0.23 eV, respectively, as shown in the inset of Figure 7a.In terms of these energy barrier values, they are much smaller than those in the case of BZY0.This difference strongly suggests that proton migration is much easier for BZY3.7 ~BZY22.2than for BZY0.In other words, doping can significantly improve the migration of protons.More descriptions of several, such as defects proton-proton repulsion and vacancy-vacancy repulsion, can be referred to in the research of F.M. Draber et al. [20,21].Furthermore, we calculated the effective Bader charge of the atoms in BZYs.The Bader charge value of Zr in BZY0 is predicted to be +2.25e, and Y in BZY3.7,BZY7.4,BZY14.8, and BZY22.2 is predicted to be +2.10e,+2.11e, +2.13e, and +2.15e, respectively.These indicate that a smaller Coulomb exclusion barrier and a weaker electrostatic potential field for proton migration may exist, as the yttrium-doped Bader charge value was lower than Zr when Y atoms lost less electronic charge than the Zr atoms at the replacement site.
Similar work by F.M. Draber et al. [20] considered the possible occupations of the initial and final B-site positions during a vacancy jump, and eight configurations were obtained, and the migration energies of oxygen vacancies in these eight cases were calculated; among them, the energy barrier of oxygen ion migration in BZY7.4 is about the same as that reported in the reference.

Electronic Structure Analysis
For deep insight into the effect of yttrium content on the electronic conductivity of BZY, the density of states (DOS) and electronic band structure for BZY were explored, as shown in Figure 9.It is evident that pure BZO exhibited an indirect band gap at R-G.The band gap is 3.16 eV, which shows consistency with previously reported results [44,45].In Figure 9a, the TDOS of BZO is mainly contributed to by Zr_d and O_p electrons, as shown in the energy range from −2 eV to 6 eV.In the valence band, the DOS of Zr_d electrons is lower than that of O_p electrons.While in the conduction band, the DOS of Zr_d electrons From the results of the proton energy barrier, it can be seen that the energy barrier of the H2 to H3 step varies with the Y content, while the energy barrier of the H1 to H2 step changes slightly.Therefore, we analyzed the saddle point in transition states in the H2 to H3 step.The distance between the proton and the nearest oxygen atom (dOH) is given in Table 2.With the transition state changes, the dOH length also changes, and the longest dOH corresponds to the saddle point.The dOH length at the saddle point of BZY0 is 1.94 Å, and the dOH at the saddle points of each content is significantly reduced with the increase in Y content.Meanwhile, the lowest length of 1.19 Å is obtained at BZY14.8.As can be seen, the dOH length at the saddle point also exhibits the same trend as the energy barrier.One of the possible reasons is that yttrium doping leads to the change in dOH, which leads to the change in structural energy, thus, affecting the migration energy barrier.a The dOH at transition state of saddle point.

Oxygen Ion Migration in BZYs
Similar to the strategy of studying proton migration, the ideal doping structures with the lowest formation energy were used to explore the energy barrier of oxygen ion migration in BZYs.Similarly, the number of oxygen vacancies created by doping is not considered.In order to form a system with the study of proton migration, the oxygen ion migration path in each doping content is still selected at the central octahedron of the BZY supercell.A schematic diagram of the oxygen ions migration path in BO 6 (B = Zr, Y) is shown in Figure 6a.The diffusion paths and the energy curves of oxygen ion migration from the initial state O1 to oxygen vacancy O2 are shown in Figures 7b-e and 8b, respectively, which were obtained using CI-NEB.The energy barriers for oxygen ion migration in BZY0, BZY7.4,BZY14.8, and BZY22.2 are 1.97, 1.08, 0.78, and 0.73 eV, respectively.This shows that Y doping can also reduce the energy barrier of oxygen ion migration.Initially, the energy barrier decreases rapidly with the increase in Y content.Specifically, from BZY0 to BZY7.4, the change in the energy barrier is 0.89 eV when the Y content reached 7.4 mol%.Then, the change in energy barrier decreases to 0.3 eV from BZY7.4 to BZY14.8.However, the change in the energy barrier between BZY14.8 and BZY22.2 is only 0.05 eV, which is without a significant difference.Therefore, the energy barrier of oxygen ion migration between BZY14.8 and BZY22.2 reached the lowest value, which is consistent with other literature [19].Moreover, the doping of Y will create lots of oxygen vacancies, which increase the concentration of charge carriers, and the increase can also improve ionic conductivity.Similar work by F.M. Draber et al. [20] considered the possible occupations of the initial and final B-site positions during a vacancy jump, and eight configurations were obtained, and the migration energies of oxygen vacancies in these eight cases were calculated; among them, the energy barrier of oxygen ion migration in BZY7.4 is about the same as that reported in the reference.
among them, the energy barrier of oxygen ion migration in BZY7.4 is about the same as that reported in the reference.

Electronic Structure Analysis
For deep insight into the effect of yttrium content on the electronic conductivity of BZY, the density of states (DOS) and electronic band structure for BZY were explored, as shown in Figure 9.It is evident that pure BZO exhibited an indirect band gap at R-G.The band gap is 3.16 eV, which shows consistency with previously reported results [44,45].In Figure 9a, the TDOS of BZO is mainly contributed to by Zr_d and O_p electrons, as shown in the energy range from −2 eV to 6 eV.In the valence band, the DOS of Zr_d electrons is lower than that of O_p electrons.While in the conduction band, the DOS of Zr_d electrons

Electronic Structure Analysis
For deep insight into the effect of yttrium content on the electronic conductivity of BZY, the density of states (DOS) and electronic band structure for BZY were explored, as shown in Figure 9.It is evident that pure BZO exhibited an indirect band gap at R-G.The band gap is 3.16 eV, which shows consistency with previously reported results [44,45].In Figure 9a, the TDOS of BZO is mainly contributed to by Zr_d and O_p electrons, as shown in the energy range from −2 eV to 6 eV.In the valence band, the DOS of Zr_d electrons is lower than that of O_p electrons.While in the conduction band, the DOS of Zr_d electrons is higher than that of O_p electrons, which indicates that some O_p electrons transfer to the conduction band and hybridize with the Zr_d electron orbitals.The DOS of Zr_d and O_p electrons indicates that Zr interacts covalently with O atoms.In contrast, the difference between the DOS peaks of the Ba and O atoms is relatively large, which indicates a weak bonding interaction between the two atoms, forming an ionic bond.Comparing Figure 9a with Figure 9b-e

Experimental Analysis
We also evaluated the total conductivities of BaZr(1-x)YxO3-δ (BZY, x = 0.1, 0.2, 0.3, and 0.4) using experiments.The influence of the yttrium doping content on the performance of BZY electrolytes was analyzed.The crystalline structure of the samples and XRD patterns of BZY powders with different yttrium contents are shown in Figure 10b.Using the standard PDF card, the peaks, which are marked with solid triangles, correspond to Ba-ZrO3, and some impurity peaks, including ZrO2 and BaCO3, were also identified and marked with open diamonds and open triangles, respectively.In addition, there was a slight shift in the main peak positions to lower angles with the increase in yttrium contents.The relationship between Goldschmidt tolerance factor (t) and ionic radius can explain this phenomenon, which can be defined by the following Equation (2):

Experimental Analysis
We also evaluated the total conductivities of BaZr (1-x) Y x O 3-δ (BZY, x = 0.1, 0.2, 0.3, and 0.4) using experiments.The influence of the yttrium doping content on the performance of BZY electrolytes was analyzed.The crystalline structure of the samples and XRD patterns of BZY powders with different yttrium contents are shown in Figure 10b.Using the standard PDF card, the peaks, which are marked with solid triangles, correspond to BaZrO 3 , and some impurity peaks, including ZrO 2 and BaCO 3 , were also identified and marked with open diamonds and open triangles, respectively.In addition, there was a slight shift in the main peak positions to lower angles with the increase in yttrium contents.The relationship between Goldschmidt tolerance factor (t) and ionic radius can explain this phenomenon, which can be defined by the following Equation ( 2 In the perovskite structure (ABO 3 ), R A , R B , and R O are the radii of cation A, cation B, and oxygen, respectively.As shown in Figure 10a, the lattice parameter increases linearly with an increase in yttrium content due to the replacement of the Zr 4+ -sites by the bigger Y 3+ ions, which is consistent with other literature [46,47].The intercept in the y-axis is 4.19922 Å, which is consistent with the lattice parameter of BaZrO 3 (4.1987 Å).This suggests that yttrium was completely doped into the BaZrO 3 lattice and forms a pure perovskite oxide in all conditions.In the perovskite structure (ABO3), RA, RB, and RO are the radii of cation A, cation B, and oxygen, respectively.As shown in Figure 10a, the lattice parameter increases linearly with an increase in yttrium content due to the replacement of the Zr 4+ -sites by the bigger Y 3+ ions, which is consistent with other literature [46,47].The intercept in the y-axis is 4.19922 Å , which is consistent with the lattice parameter of BaZrO3 (4.1987 Å ).This suggests that yttrium was completely doped the BaZrO3 lattice and forms a pure perovskite oxide in all conditions.
The conductivities of the above-synthesized BZY electrolytes were measured using the EIS method in moist air (3% H2O + air), and the testing temperature was 450 °C .All spectra include a recessed arc and a tail; these spectra were fitted using ZView2 software and analyzed with an equivalent circuit, which consisted of a resistor (R) and two pairs of parallel resistor (R)-constant phase elements (CPE) connected in a series, as shown in the inset of Figure 11.The intercept point at high frequency (1 MHz) represents the grain resistance (Rg) of the electrolyte, with the first semicircular arc extended to intersect the xaxis and the value of the horizontal coordinate being the sum of the grain (Rg) and grain boundary (Rgb) resistances.The gas phase diffusion occurring at the electrode and the electrochemical reactions occurring at the TPBs can be revealed by the low-frequency tails.They were also described by the charge transfer resistance (Rct) and the constant phase angle element (CPEct).As the yttrium contents increased from 10 to 40 mol%, the conductivity of BZYs increased first and then decreased, and the maximum value was obtained The conductivities of the above-synthesized BZY electrolytes were measured using the EIS method in moist air (3% H 2 O + air), and the testing temperature was 450 • C. All spectra include a recessed arc and a tail; these spectra were fitted using ZView2 software and analyzed with an equivalent circuit, which consisted of a resistor (R) and two pairs of parallel resistor (R)-constant phase elements (CPE) connected in a series, as shown in the inset of Figure 11.The intercept point at high frequency (1 MHz) represents the grain resistance (R g ) of the electrolyte, with the first semicircular arc extended to intersect the x-axis and the value of the horizontal coordinate being the sum of the grain (R g ) and grain boundary (R gb ) resistances.The gas phase diffusion occurring at the electrode and the electrochemical reactions occurring at the TPBs can be revealed by the low-frequency tails.They were also described by the charge transfer resistance (R ct ) and the constant phase angle element (CPE ct ).As the yttrium contents increased from 10 to 40 mol%, the conductivity of BZYs increased first and then decreased, and the maximum value was obtained at BZY20, which can be seen in the inset table in Figure 11.These results showed the same regulation with the barrier energies in the calculation results and are consistent with other literature [47].
at BZY20, which can be seen in the inset table in Figure 11.These results showed the same regulation with the barrier energies in the calculation results and are consistent with other literature [47].

Conclusions
The evolution of protonic, oxygen ionic, and electronic conductivities in BZY electrolyte were measured using DFT calculations and experimental tests.Firstly, the energy barriers of proton and oxygen ion migration were effectively reduced by Yttrium doping, and the lowest values were 0.14 eV (Proton) for BZY14.8 and 0.73 eV (oxygen ion) for BZY22.2.The electronic band structure indicates that unfavorable electronic conductivity emerged with Y doping.Secondly, as the yttrium content increased from 10 to 40 mol%, the total conductivity of BZYs increased first and then decreased, and the maximum value was obtained in BZY20.Finally, the tendency of the tested total conductivity is consistent with the result from the DFT calculation.This work provides support to promote the wide application of BZY in PCFC.

Conclusions
The evolution of protonic, oxygen ionic, and electronic conductivities in BZY electrolyte were measured using DFT calculations and experimental tests.Firstly, the energy barriers of proton and oxygen ion migration were effectively reduced by Yttrium doping, and the lowest values were 0.14 eV (Proton) for BZY14.8 and 0.73 eV (oxygen ion) for BZY22.2.The electronic band structure indicates that unfavorable electronic conductivity emerged with the Y doping.Secondly, as the yttrium content increased from 10 to 40 mol%, the total conductivity of BZYs increased first and then decreased, and the maximum value was obtained in BZY20.Finally, the tendency of the tested total conductivity is consistent with the result from the DFT calculation.This work provides support to promote the wide application of BZY in PCFC.

Figure 1 .
Figure 1.(a) Sample photograph.(b) The front view of the sample.(c) The side view of the sample.

Figure 1 .
Figure 1.(a) Sample photograph.(b) The front view of the sample.(c) The side view of the sample.
(1-x) Y x O 3-δ (BZY, x = 0, 0.037, 0.074, 0.148, and 0.22) are shown in Figure 4. BZY0 (BaZrO 3 ) represents the BZY without Y atom replacement.BZY3.7 represents the BZY with one Y atom replacement, where the atom is set in the center of the supercell.In Figure 4c-e, the supercell with two Y atoms is illustrated.It can be seen that the distance between the two of them is a, √ 2a, and √ 3a.The next step is four Zr atoms, which are replaced by four Y atoms, as can be seen in Figure 4f-i.The total distances between the three surrounding Y atoms to the centered Y atom are 3 × a, 3 × √ 2a, and 3 × √ 3a, which cover all conditions of the four Y atoms' replacement.Furthermore, the six Y atoms' replacements are shown in Figure 4j-o.As well as the situations in the four Y atoms' replacement, we covered all arrangements by controlling the total atoms distance from 5 × a to 5 × √ 3a.Crystals 2023, 13, 401 5 of 14 Crystals 2023, 13, x FOR PEER REVIEW 5 of 14 in Figure 4j-o.As well as the situations in the four Y atoms' replacement, we covered all arrangements by controlling the total atoms distance from 5 × a to 5 × √3a.

Figure 3 .
Figure 3. (a) Schematic 3 × 3 × 3 supercell structure of BZO.(b) Color diagram of B-site atoms according to the doping site distance distribution.The Y1 (yellow B-site atom) is the center, with Y2 (purple B-site atoms) as the nearest neighbor, Y3 (orange B-site atoms) as the next closest neighbor, and Y4 (brown B-site atoms) as the farthest.

Figure 3 .
Figure 3. (a) Schematic 3 × 3 × 3 supercell structure of BZO.(b) Color diagram of B-site atoms according to the doping site distance distribution.The Y1 (yellow B-site atom) is the center, with Y2 (purple B-site atoms) as the nearest neighbor, Y3 (orange B-site atoms) as the next closest neighbor, and Y4 (brown B-site atoms) as the farthest.

Figure 3 .
Figure 3. (a) Schematic 3 × 3 × 3 supercell structure of BZO.(b) Color diagram of B-site atoms according to the doping site distance distribution.The Y1 (yellow B-site atom) is the center, with Y2 (purple B-site atoms) as the nearest neighbor, Y3 (orange B-site atoms) as the next closest neighbor, and Y4 (brown B-site atoms) as the farthest.

Figure 8 .
Figure 8. Energy curves of (a) proton and (b) oxygen ion migration in BZYs.

Figure 8 .
Figure 8. Energy curves of (a) proton and (b) oxygen ion migration in BZYs.

Figure 8 .
Figure 8. Energy curves of (a) proton and (b) oxygen ion migration in BZYs.
, the top of the valence band of BZO is below the Fermi energy level, and the DOS of BZO is changed by yttrium doping.Due to the yttrium-doped BZO, a donor level was formed, resulting in the Fermi level moving closer to the valence band in the BZY3.7 ~BZY22.2.This facilitates electron conduction and, in turn, decreases the activation energy of electron conduction.The DOS of d electrons for the B-site ions increased around the Fermi level and overlapped with the oxygen p electrons under the situation of yttrium doping, which may imply an increase in charge carrier concentration.Both of these conditions can enhance electron conduction.ence between the DOS peaks of the Ba and O atoms is relatively large, which indicates a weak bonding interaction between the two atoms, forming an ionic bond.Comparing Figure9awith Figure9b-e, the top of the valence band of BZO is below the Fermi energy level, and the DOS of BZO is changed by yttrium doping.Due to the yttrium-doped BZO, a donor level was formed, resulting in the Fermi level moving closer to the valence band in the BZY3.7 ~ BZY22.2.This facilitates electron conduction and, in turn, decreases the activation energy of electron conduction.The DOS of d electrons for the B-site ions increased around the Fermi level and overlapped with the oxygen p electrons under the situation of yttrium doping, which may imply an increase in charge carrier concentration.Both of these conditions can enhance electron conduction.

Figure 9 .
Figure 9. Electronic band structure and density of states for BZYs.

Figure 9 .
Figure 9. Electronic band structure and density of states for BZYs.

Figure 10 .
Figure 10.(a) Relationship between the calculated lattice parameter and the yttrium content for BZYs.(b) XRD patterns of calcined BZYs powders.

Figure 10 .
Figure 10.(a) Relationship between the calculated lattice parameter and the yttrium content for BZYs.(b) XRD patterns of calcined BZYs powders.

Figure 11 .
Figure 11.Nyquist plot of BZYs in a humid air atmosphere and at 450℃.The conductivity of BZYs with different yttrium doping contents is shown in the top right table.

Figure 11 .
Figure 11.Nyquist plot of BZYs in a humid air atmosphere and at 450°C.The conductivity of BZYs with different yttrium doping contents is shown in the top right table.

Table 1 .
Formation energy of BZY structures with different doping contents; the structures screened out by the behaviors marked in gray are used for calculation.

Table 1 .
Formation energy of BZY structures with different doping contents; the structures screened out by the behaviors marked in gray are used for calculation.
# of Y

Table 2 .
Distance between Hydrogen and the Nearest Oxygen Atom (dOH) at Transition States of Proton Transfer Step H2-H3 in BZY.