The Impact of Nanoscale Percolation in Yttrium‐Doped BaZrO3 on the Oxygen Ion and Proton Conductivities: A Density Functional Theory and Kinetic Monte Carlo Study

Yttrium‐doped BaZrO3 is known for its high proton conductivity, making it an advanced energy material for various applications, like electrolyzers, fuel cells, or methane conversion cells. In the previous study, density functional theory (DFT) and kinetic Monte Carlo simulations (KMC) show that nanoscale percolation of Y ions enhances the proton mobility by providing fast migration pathways for protons. To investigate the general impact of nanoscale percolation on ionic conductivities, the same methods, DFT and KMC, are now used to calculate oxygen ion conductivity of Y‐doped BaZrO3. The results explain on a microscopic level why the macroscopic oxygen ion conductivity exhibits a higher activation energy than the proton conductivity, as known experimentally. They also show that nanoscale percolation pathways are not beneficial for oxygen ion conduction due to trapping of oxygen vacancies by single Y ions and blocking of their motion by nanoscale percolation of Y ions. Finally, the understanding and comparison of the microscopic jump processes of oxygen vacancies and of protons in acceptor‐doped BaZrO3 are used to propose a simple descriptor for the influence of a dopant on the ionic conductivities, of oxygen ions and of protons as well. This new descriptor allows an easy screening of various dopants.

nanoscale percolation of Y ions is beneficial or detrimental for the oxygen ion conductivity, despite being lower than the proton conductivity; and 3) By comparing the resulting detailed, microscopic picture for oxygen vacancies with our previous results for protons, we will propose for acceptor-doped BaZrO 3 a simple, new descriptor for the influence of a dopant on the ionic conductivities, of protons and of oxygen ions as well. This new descriptor will allow an easy screening of various new dopants by performing only a few DFT calculations.
There are only few studies that investigated oxygen vacancies [10,11] or even the oxygen ion conductivity in Y-doped BaZrO 3 . [2,[11][12][13][14] They reported activation energies of the oxygen ion conductivity between 0.93 and 1.2 eV, depending on the dopant fraction. In comparison, the activation energies of the proton conductivity in Y-doped BaZrO 3 are typically around 0.4 eV. [1,9,[14][15][16][17][18][19] In this article, we investigate the oxygen vacancy migration on a microscopic level and the resulting macroscopic oxygen ion conductivity in Y-doped BaZrO 3 under dry conditions using the same methods that we had applied successfully to investigate the microscopic proton migration and the resulting macroscopic proton conductivity. [9] For this purpose, we combine again two theoretical methods. We use density functional theory (DFT) to obtain the defect interaction energies, in particular the trapping energies of oxygen vacancies trapped by Y ions, and the migration energies of oxygen vacancies as well. Based on these energies we perform kinetic Monte Carlo (KMC) simulations to predict the resulting macroscopic oxygen ion conductivity. In our previous paper [9] on the proton conductivity of Y-doped BaZrO 3 , we considered the fully hydrated case. In this article, we will treat the other extreme, that is, oxygen ion conductivity under dry conditions. The general case of partially hydrated, Y-doped BaZrO 3containing mobile protons and mobile oxygen vacancies and considering all possible defect interactions will be treated in a forthcoming paper.
The article is organized as follows. In Section 2, we describe the defects in acceptor-doped BaZrO 3 , and in Section 3, we introduce the computational setup for DFT and KMC calculations. In Section 4, we present and discuss our results and compare oxygen ion and proton conductivities and their different dependencies on dopant fraction. Finally, in Section 5, we propose a simple, new descriptor for the influence of a dopant on the ionic conductivities of acceptor-doped BaZrO 3 .

Defects in Y-Doped BaZrO 3
Pure BaZrO 3 is doped at the B-site with Y ions, acting as acceptors, Y ' Zr , and introducing oxygen vacancies, V •• O , in the material (Kröger-Vink notation).
An Y ion may trap an oxygen vacancy, which is often described as a dopant-vacancy pair fY ' Zr ,V •• O g•. The clear definition of such a pair and the use of the corresponding law of mass action for the pair formation reaction Y ' Zr þV •• O ⇄ fY ' Zr ,V •• O g• is, however, only possible for rather small acceptor fractions. In the ABO 3 perovskite structure, the acceptor position (B-site) is coordinated by six oxygen sites, and one can easily estimate that isolated acceptor-vacancy pairs are only possible for defect fractions below 4%. For higher acceptor levels, the pairs start to overlap, and the usage of the pair concept leads to wrong results. As we want to describe the oxygen ion conductivity for the full range of acceptor levels, that is, from diluted acceptors to concentrated solid solutions, we do not use the concept of pair formation. Instead we will consider in DFT and KMC various local configurations of acceptor ions and oxygen vacancies and the corresponding interaction energies. As a consequence, we do not consider free and trapped vacancies separately, and the charge neutrality simply reads as total fractions of Y ions and oxygen vacancies. In the charge neutrality condition, we have assumed that electron holes, h • , that are polarons located on oxygen ions [20] and formed via ½O 2 Þ are only minority defects and can be neglected. At sufficiently high oxygen partial pressure, electron holes might, however, contribute to the experimental total conductivity due to their high mobility.
The transport mechanism of oxygen ions is believed to be a simple vacancy mechanism where the oxygen ion passes through a transition state which is energetically less favorable than the initial or final state. The migration energy for this process is given by E mig ¼ E ts À E ini , where E ts and E ini are the system energies in the transition state and the initial state, respectively. In the example shown in Figure 1, an oxygen vacancy jumps from a first-nearest-neighbor position (1NN) relative to a Y ion to a second-nearest-neighbor position (2NN). As will be shown later, trapping on the 2NN position is weaker than on the 1NN position, resulting in the asymmetry of the energy profile shown in Figure 1.

Ab Initio Calculations
All site and migration energies of oxygen vacancies were calculated by means of DFT within the generalized gradient approximation (GGA) according to Perdew, Burke, und Ernzerhof (PBE) [21] and the projector augmented wave (PAW) method, [22] as implemented in the Vienna Ab initio Simulation Package (VASP). [23,24] The climbing image-nudged elastic band method (CI-NEB) [25,26] was used to investigate the transition states and the minimum energy pathways. Plane waves with an energy cutoff of 500 eV and a 6 Â 6 Â 6 Monkhorst-Pack k-point mesh [27] were used for the Pm3 À m unit cell. The convergence parameters for the electronic and ionic relaxation were set to 10 À5 eV and 10 À2 eV Å À1 , respectively. The valence electrons for all elements can be seen in Table 1.
For all defective cells, the number of electrons was adapted to reproduce the real charge state of the defect. Thus, depending on the considered defects within the cells, they are positively charged, neutral, or negatively charged. Examples for all three cases are a cell containing the two defects V •• O and Y ' Zr , a cell containing one V •• O and two Y ' Zr , and a cell containing two Y ' Zr .

Size of the Supercell
Various seized supercells were used during the DFT calculations. The unit cell of BaZrO 3 contains five atoms. Supercells are built as multiples of this unit cell. In calculating energies of cells containing a defect, it is necessary to perform a finite size correction, to ensure that the defect is not influenced by its own image. [28] To do so, the energy is calculated in three supercells with different volumes V, and by plotting these energies against 1/V, it is possible to extrapolate to the energy corresponding to 1/V ¼ 0, meaning infinite dilution. Here, no interaction between a defect and its own image is possible. Two defects within the supercell exhibit long-range interactions (see later), and it is necessary to use supercells of at least the size of the 240 atom supercell. In addition, care has been taken to choose only those supercells that allow tilting of the oxygen octahedra. [29] Tilting is well known for BaZrO 3 containing protons. [7,9] Nevertheless, it appears also in supercells with Y dopants and oxygen vacancies (see later). In our calculations, we used 2 ffiffi ffi 2 p Â 3 ffiffi ffi 2 p Â 4 (240 atoms), 4 Â 4 Â 4 (320 atoms), and 3 ffiffi ffi 2 p Â 3 ffiffi ffi 2 p Â 4 (360 atoms) supercells, all of which fulfill the earlier requirements.

Kinetic Monte Carlo Simulations
KMC simulations were performed using the program MOCASSIN. [30] The Monte Carlo algorithm [31] was implemented in C# and random numbers were generated by the Minimal C Implementation of the permuted congruential generator (PCG). [32] All jumps were carried out in a supercell containing 16 Â 16 Â 16 unit cells with periodic boundary conditions. Following the method of Grope et al., [33] a weak electric field (E field ) was applied, which led to an average displacement Δx of the oxygen vacancies in field direction. Then, the average vacancy mobility u V •• O can be calculated as The elapsed time is given by with N att as the number of jump attempts, N poss as the number of jump possibilities, N V as the number of vacancies, and ν 0 as the attempt frequency of a jump. The oxygen ion conductivity σ O 2À is then given by where q V •• O and c V •• O are the charge and concentration of V •• O . KMC simulations were carried out in a simulation cell consisting of 16 Â 16 Â 16 unit cells of BaZrO 3 and using periodic boundary conditions. To achieve good accuracy, jump attempts were carried out until every oxygen site was visited by a vacancy 1000 times on average. This value is called Monte Carlo Steps per Particle (MCSP). The initial starting configuration for the (immobile) dopant ions and the oxygen vacancies was completely random. To equilibrate the oxygen vacancies before data collection for the ionic conductivity, initially 100 MCSP were carried out.

Defect Interaction Energies
We performed DFT calculations for various configurations of defect pairs, as shown schematically in Figure 2. In this pair interaction model, there are three basic types of pair interactions:  Zr -Y ' Zr (not shown). In each DFT calculation, the distance between the corresponding two defects was varied, and a structural relaxation was performed. For the finite size correction, supercells of different sizes were used (see Section 3.2), all of them containing only a single defect pair.
For the pair interaction between an oxygen vacancy and an yttrium dopant, an attractive interaction between these two oppositely charged species was found, which weakens with increasing distance. The finite size-corrected values (see Figure S1 in the Supporting Information) are shown in Figure 3. Because the supercell energy does not change between a distance of 7 and 10 Å, we set the interaction energy for the largest calculated distance to 0.0 eV, assuming that there is no more interaction between both species. All other interaction energies are calibrated relative to this reference point. Figure 3 shows that a single Y ion traps an oxygen vacancy on 1NN and 2NN sites and also weakly on 3NN sites. The trapping energy on 3NN is À0.07 eV. This value is so small, that this position is not counted as a trapped position, and we model the exact results in Figure 3 by considering trapping on 1NN and 2NN only, that is, up to a distance of about 6 Å around every yttrium dopant. The corresponding spherical volume with a radius of this distance is the trapping zone. These trapping zones start to over- Zr ] > 0.015. Only for smaller dopant fractions the concept of isolated dopant-vacancy pairs could be used. In these calculations, problems can emerge from finite size effects, meaning that a defect is affected by its own image in the next supercell. To avoid this, we used supercells that have a supercell size of at least twice the interaction radius between an yttrium ion and an oxygen vacancy.
We investigated the pair interaction between two oxygen vacancies analogously. This interaction was expected to be repulsive due to the same charge of both species, but Figure 4 does not show the expected dependence on distance due to coulombic pair interactions alone. The interaction energy does not simply decrease with increasing distance but depends also on lattice distortions. This effect can be visualized for the interaction between two vacancies at a distance of 4.23 Å (2NN). There are two different values, which differ considerably. This can be traced back to two different 2NN configurations that are possible with the same distance between both vacancies. In the first one (2aNN in Figure 2) with a zirconium ion between both vacancies, the interaction energy is rather small. In the other configuration (2bNN in Figure 2) with a zirconium ion between them, the interaction energy of about 2.3 eV is very high. A similar effect has been observed in doped ceria. [34] The same effect can be observed in alleviated terms, when one more unit cell lies between both Figure 2. Extract of the BaZrO 3 structure with yttrium in blue, zirconium in green, oxygen in red, and oxygen vacancies as black squares. The pair interaction energies that are most important for the KMC simulations are shown. Each pair interaction between defect 1 and defect 2 is described in terms of the neighborhood of defect 2 relative to defect 1. There are two pair interactions for V , which are denoted as 2aNN (with zirconium in between) and 2bNN (without zirconium in between). www.advancedsciencenews.com www.advenergysustres.com vacancies at a distance of 8.46 Å. Here, two interaction energies of 0.07 and 0.42 eV were found. Introducing one oxygen vacancy into the structure leads to distortion of the surrounding oxygen ions. This has a stabilizing and therefore energy-lowering effect on the structure. Introducing two oxygen vacancies into the structure leads to an intensification or attenuation of the distortion, depending on the positions of the vacancies relative to each other. In the structure with one zirconium ion between both vacancies, the distortion is very small due to symmetry reasons, which contributes to the very high site energy (see Figure S2, Supporting Information). The detailed pair interactions in Figure 3 and 4 are included in the energy model that was used in the subsequent KMC simulations. The interaction between two yttrium ions was also investigated analogously. It was found that only the interaction in 1NN is important with an interaction energy of þ0.2 eV, and interactions at longer distances could be neglected.
We checked our pair interaction model furthermore by analyzing a supercell containing one oxygen vacancy and two yttrium ions. We calculated the total energy of that supercell for several configurations of the three defects. This "exact" energy contains also many-body interactions due to local distortions of the structure. In comparing with our model that assumes pair interactions only (Table S1, Supporting Information), we found good agreement between exact and model results. Only for a cluster consisting of one oxygen vacancy trapped between two yttrium dopants on 1NN positions (linear defect triplet Zr )) we found a significant difference. In the pair interaction model, this vacancy is trapped with an energy of 2 · (À0.36) eV ¼ À0.72 eV (see Figure 3), and with an Y ' Zr -Y ' Zr interaction energy of þ0.2 eV the total supercell energy is À0.52 eV. In contrast, the exact total supercell energy is À0.74 eV. This means, that for the oxygen vacancy, the trap is 0.22 eV deeper than expected using the pure pair interaction model. We call this case "deep trapping" in the following, and it was considered explicitly in the energy model for the KMC simulations.
We also found less good agreement between the pair interaction model and the "exact" results with both yttrium ions in 2NN to the oxygen vacancy but far away from each other. This can be traced back to the missing tilting of the oxygen octahedra. If the defects are placed in a way that the tilting is hindered (e.g., in a row), then the pair interaction model systematically underestimates the energy of the cell. In the case of Y1-Vö in 2NN, Y2-Vö in 2NN, and Y1-Y2 in 5NN, the oxygen ions shift during the relaxation by only one-thirds of the corresponding case with Y1-Y2 in 3NN. In the perovskite structure, these configurations are rare and we stick also for them to the pair interaction model.  Figure 5. The migration energies of an oxygen vacancy jump were calculated for these eight cases. Furthermore, longer oxygen vacancy jumps are imaginable, but we found that the migration barriers of these jumps are typically >3 eV, allowing us to neglect them.

Migration Energies of Oxygen Vacancies
For the finite size correction after Makov and Payne, [28] we calculated the migration barriers for each of the eight configurations in Figure 5 in three differently sized supercells containing the same number of defects (Y ions and oxygen vacancies). The detailed results of the finite size correction are shown in Figure S3 (Supporting Information) and are listed in Table S2 (Supporting Information). The finite size-corrected barriers shown in Figure 5 vary between 0.46 and 1.16 eV depending on the dopant configuration.
In Figure 5, some results are remarkable. If the zirconium position at the top of the triangle is doped with yttrium (Figure 5b), the migration energy is 0.5 eV larger than in the corresponding configuration with zirconium at the top position (Figure 5a). One possible reason for this large value is the higher ionic radius of Y 3þ : R Y3þ ¼ 0.90 Å compared with R Zr4þ ¼ 0.72 Å. [35] The structures obtained from DFT show that in the transition state the distance between the moving oxygen ion and the Y ion at the top of the triangle in Figure 5b is about 0.13 Å higher compared with the occupation with a Zr ion. Accordingly, the "window" formed by the cation on the top B-position and both neighboring barium ions becomes smaller, which enlarges the migration energy.
The difference of the migration energies for a forward jump and the corresponding backward jump, ΔE mig , must be identical to the difference of both site energies ΔE site (microscopic balance). For the jump in Figure 5c and its reverse jump in Figure 5e, we obtain ΔE mig ¼ 0.08 eV. Because both configurations contain only a single Y ion, this value is in perfect agreement with the difference of the site energies for a vacancy on 1NN and 2NN positions relative to the Y ion calculated from    Figure 5. Possible (Zr,Y) configurations for oxygen vacancy jumps in BaZrO 3 and corresponding migration energies. Zirconium is shown in green, yttrium in blue, oxygen in red, and the oxygen vacancy as a black square. The shown values for the migration energies are the values extrapolated to infinite dilution. They were calculated with DFT and finite size corrected after Makov and Payne. [28] .
www.advancedsciencenews.com www.advenergysustres.com Y ions are involved (Figure 5d,f-h), the migration energies are higher than in a pure Zr environment (Figure 5a) but still smaller than for the 1NN!1NN jump around a single Y ion (Figure 5b). In literature, only the vacancy jump in Figure 5a was calculated by Stokes et al. using semiempirical methods. [36] Their calculated value 0.65 eV agrees well with our migration barrier in undoped BaZrO 3 . Figure 5 shows, however, that Y doping causes a complex microscopic behavior depending on the local environment.
All information from Figure 3-5 is brought together in Figure 6, showing a jump path of one oxygen vacancy through a possible dopant configuration in Y-doped BaZrO 3 .

Kinetic Monte Carlo Simulations
The oxygen ion conductivity of Y-doped BaZrO 3 was calculated as a function of the dopant fraction and the temperature using KMC simulations. As input parameters for the KMC simulations, the finite size-corrected migration energies in Figure 5 and the longrange defect interaction energies in Figure 3 and 4 were used. Most experimental studies of BaZrO 3 consider dopant fractions of x ¼ 0.1 or 0.2, [1,2,8,17,37] while there only few studies with dopant fractions up to x ¼ 0.3. [38,39] Although the solubility limit of Y is at x % 0.4, simulations for larger values of x up to x ¼ 1 are meaningful, as they yield deeper insight into the microscopic origin of the macroscopic conductivity and will help to identify descriptors for the influence of dopants and conductivity (see Section 5).
The simulated oxygen vacancy mobilities and oxygen ion conductivities are shown in Figure 7  At low dopant fractions, x, the oxygen vacancy mobility (Figure 7) decreases rapidly with increasing x due to trapping of oxygen vacancies on 1NN and 2 NN positions relative to Y ions and also due to blocking of 1NN ! 1NN vacancy jumps "around" Y ions (see Figure 5b). Both effects, trapping and blocking, become more pronounced with decreasing temperature. The corresponding oxygen ion conductivity ( Figure 8) increases with increasing x at high temperatures, due to the increasing fraction of oxygen vacancies, see Equation (5), while at low temperatures, it decreases with x, because the strong drop in the vacancy mobility cannot be compensated by the increase in the vacancy fraction. Li et al. [40] found with reactive molecular dynamics simulations also an increase of the oxygen ion conductivity with increasing x until x ¼ 0.2 for temperatures above 1000 K and a decrease for lower temperatures. They explained this result only with more effective trapping at low temperatures.
At intermediate dopant fractions, the random distribution of Y ions results in the formation of Y pairs and longer local  Figure 5 and the interaction energies in Figure 3, as described in the model of Grieshammer et al. [34]   www.advancedsciencenews.com www.advenergysustres.com structures (nanoscale percolation of Y ions, see Figure S4, Supporting Information). In these nanoscale structures, "deep traps" (an oxygen vacancy is trapped between two neighboring Y ions) become important. An oxygen vacancy in this deep trap has only a small probability of leaving it again during the simulation. The mobility decreases only slightly with increasing dopant fraction, and the oxygen ion conductivity reaches a plateau.
For very high dopant fractions (x > 0.7), the vacancy mobility and the oxygen ion conductivity increase again. Nearly every oxygen position is coordinated by Y ions, and these structures can form complete percolation pathways through the simulation cell. The migration energy for jumps in Y-Y-Y-configuration is 0.88 eV (see Figure 5h), which is rather high compared with Zr-Zr-Zr configuration (0.66 eV), but still lower than the energy for jumps out of a deep trap (1.04 eV, see Figure 6). At these high dopant fractions, the energy landscape and the resulting migration energies are, however, much more complicated than in Figure 5 due to the high vacancy concentration and the longrange vacancy-vacancy interaction (see Figure 4). The resulting large number of possible configurations renders a simple analysis impossible, but the influence of these configurations is of course correctly included in the energy model for the KMC simulations.
The activation energies of the macroscopic oxygen ion conductivity were obtained from an Arrhenius plot of the simulated oxygen ion conductivities ( Figure S5, see Supporting Information) and compared with experimental and theoretical studies from literature (see Table 2). Our simulated activation energies fit well to those data for x ¼ 0.1 and x ¼ 0.2 considering the experimental error bar [14] and the fact that no error bars were specified for the other literature values.
We report at this point another interesting finding. Assuming the usual temperature dependence for ion conduction, σ(T ) ¼ (A/T )exp(ÀE a /k B T ), the KMC simulations yield a linear correlation between the pre-exponential factor, ln(A), and the activation energy, E a (see Figure S6, Supporting Information). This relation is well known as the Meyer-Neldel rule. [41] In literature, there are several models explaining this empirical rule. A detailed investigation is beyond the scope of this article, but we note that the broad distribution of migration energies in this highly disordered oxide may be the origin of the observed correlation between ln(A) and E a .

Comparison of Simulated Oxygen Ion and Proton Mobilities and Conductivities
Comparing the simulated oxygen ion conductivity in yttriumdoped BaZrO 3 under dry conditions with proton conductivity in the fully hydrated oxide, [9] we find some similarities but also important differences.
We start with the similarities. An oxygen vacancy is trapped in 1NN and 2NN of a Y ion, just like a proton (Figure 9a,b), both with only small differences in the depth of the trap on 1NN and 2NN positions. The fact that the trapping energies of oxygen vacancies and protons differ only slightly shows clearly that not only the direct Coulomb interaction but also lattice distortions are responsible for the trapping. In addition, both defect species, oxygen vacancies and protons, repel each other (Figure 9c,d), in particular when close together.
Differences exist in the interactions for configurations with more than one Y ion: For protons, the trapping phenomenon is broken down to a two-level system, that is, around yttrium dopants, the energy is lowered, but always approximately by the same amount. The situation is much more complicated for oxygen vacancies, where the trapping energies of several yttrium dopants are approximately additive, leading to a multilevel system. An important example is "deep trapping" of an oxygen vacancy between two Y ions, while there is no deep trapping for protons.
Oxygen vacancies show jumps between oxygen positions (translational), while protons perform rotational jumps around an oxygen ion and translational jumps between oxygen ions. For the oxygen ions, the migration energies corresponding to the different Y configurations in Figure 5 cover a range from 0.46 to 1.16 eV. In contrast, the migration energies of protons for the same Y configurations are always smaller and cover a range from 0.22 to 0.66 eV [9] (see Table 3). Due to their lower migration energies, protons move in general much faster in the material, resulting in higher ionic conductivities and mobilities.   [18] 0.1 1.1 exp. [13] 0.2 1.06 AE 0.01 DFT þ KMC (this work) 0.2 0.95 ReaxFF force field [12] 0.2 0.93 exp. [18] 0.2 1.2 AE 0.2 exp. [14] www.advancedsciencenews.com www.advenergysustres.com In summary, we emphasize at this point the most important, qualitative difference: Like protons, oxygen vacancies get trapped by yttrium dopants, but unlike them, they are not accelerated in these trapping zones. Instead, dopant ions block their way by means of high migration energies.
The Zr-Y-Zr jump (a in Table 3) has a very low migration energy for protons (0.22 eV or 54% of the jump in pure BaZrO 3 ) and a very high migration energy for oxygen vacancies (1.16 eV or 176% of the jump in pure BaZrO 3 ). This leads to an acceleration of protons around yttrium dopants while oxygen vacancies are blocked by them. The migration energies for jumps from 1NN to 2NN and backward (c: Y-Zr-Zr and e: Zr-Zr-Y) or (d: Y-Y-Zr and f: Zr-Y-Y), are the same for protons, supporting the two-level model, that was described above. For oxygen vacancies, there are, however, differences between these forward and backward jumps, and jumps with yttrium dopants on the central position (d: Y-Y-Zr and f: Zr-Y-Y-Zr) are also blocked, that is, they have migration energies higher than in pure BaZrO 3 .
All the aforementioned points are brought together in Figure 10. The higher (lower) migration energies for oxygen vacancies (protons) inside of the trapping zone are shown as well as the differences in the trapping behavior.
The differences in trapping and migration energies lead to the very different mobility curves of protons and oxygen vacancies: The mobility of the oxygen vacancies drops strongly with increasing dopant fraction for all temperatures (Figure 11a) due to trapping and blocking. In contrast, the proton mobility (Figure 11b) goes through a minimum around x ¼ 0.04 and increases again.  Table 3. Comparison of migration energies of proton jumps (E mig,H ) and oxygen vacancy jumps (E mig,Vo ). Proton migration energies are taken from Draber et al. [9] The nomenclature follows Figure 5. The defect species jumps from its position between the left two B-cations (e.g., Y and Zr in Y-Zr-Zr) to the position between the right two B-cations (Zr and Zr in Y-Zr-Zr). The defect species can either be an oxygen vacancy (this work) or a proton attached to an oxygen ion (Draber et al. [9] ) in which case only the proton jumps. Additional columns show the percentage of the respective migration energy compared with the undoped case Zr-Zr-Zr. The increase is caused by the lack of blocking of proton jumps around yttrium dopants and the low migration energies inside of overlapping trapping zones (nanoscale percolation [9] ). The resulting oxygen ion conductivity decreases for temperatures of 1000 K and lower and increases initially and then reaches a plateau for T > 1000 K. In contrast, the proton conductivity increases strongly with increasing dopant fraction for all temperatures.
In summary, nanoscale percolation pathways formed by the Y dopants are beneficial for proton mobility but blocking for the mobility of oxygen vacancies.

New Descriptor for the Influence of a Dopant on Ionic Conductivities
Finally, we use our detailed results for the microscopic migration energies and the resulting macroscopic ionic conductivities in Y-doped BaZrO 3 to propose a simple, new descriptor for predicting the influence of other dopants on the ionic conductivity of acceptor-doped BaZrO 3 . On first sight, the optimum acceptor dopant would not act as a trap for the mobile defects (oxygen vacancies or protons) and would produce only chargecompensating defects. Then, the defect mobility would not depend on the dopant fraction (except for site-blocking effects at very high defect fractions) and would correspond to the mobility at infinite dilution. The lack of trapping is, however, rather improbable, because the different charge of the acceptor dopant and its different ionic radius will most probably result in the trapping of oxygen vacancies or protons. Accepting that trapping is unavoidable, Figure 11b shows, however, that even with trapping the defect mobility can be increased compared with its value at infinite dilution. There are two necessary conditions for this behavior. 1) The defect migration energies within the trapping zone must be smaller than outside of the trapping zone (see Table 3); 2) The trapping zones must overlap to form nanoscale percolation pathways with fast defect transport within the connected trapping zones. Condition (1) is determined by the nature of the dopant, while condition (2) can be achieved by adjusting the dopant fraction to obtain overlapping trapping zones (see Figure 10). If trapping occurs on 1NN and 2NN positions (as we found for the dopant yttrium), the minimum of the defect mobility would occur at lower dopant fractions compared with trapping only on 1NN positions. As a DFT-based descriptor of the behavior of a new dopant, we propose to calculate by means of DFT only two defect migration energies: 1) in undoped BaZrO 3 (see Figure 5a); and 2) "around" the dopant, that is, for the 1NN!1NN jump (see Figure 5b). If the second value is smaller than the first one, the dopant is a good candidate for Zirconium ions are shown in green, oxygen atoms in red, and the oxygen vacancy as black square. The migration energies E mig of each jump are given below the jump. These values correspond to the combination of the migration energies in Figure 6 and the interaction energies in Figure 3 as described in the model of Grieshammer et al. [34] b) Proton jumps are shown outside and inside of the trapping zone. The colors are the same as above. Taken from Draber et al. Nature Materials, 2020, 19, 338, Springer Nature. [9] www.advancedsciencenews.com www.advenergysustres.com an improved defect mobility, that is larger than the defect mobility at infinite dilution, once the dopant fraction is high enough for the onset of nanoscale percolation. We believe that this simple descriptor that needs only two migration energies represents the necessary condition for the beneficial impact of nanoscale percolation of dopant ions on the ionic conductivity, for oxygen ions and for protons as well. This new descriptor will allow an easy screening of various new dopants.

Conclusion
In this article, we investigated the impact of nanoscale percolation in yttrium-doped BaZrO 3 on its oxygen ion and proton conductivities. For this purpose, we investigated the oxygen vacancy migration on a microscopic level and the resulting macroscopic oxygen ion conductivity using the same methods that we applied successfully to investigate the microscopic proton migration and the resulting macroscopic proton conductivity. [9] We combined DFT calculations and KMC techniques to simulate the macroscopic oxygen ion conductivity of yttrium-doped BaZrO 3 . We used DFT to calculate long-range interaction energies between one yttrium dopant and one oxygen vacancy, between two oxygen vacancies, as well as between two yttrium ions. We found a trapping zone around the yttrium dopants where oxygen vacancies are trapped on 1NN and 2NN positions and one "deep trap" between two yttrium dopants. The interaction between two oxygen vacancies turned out to be more complicated than the expected Coulomb behavior due to structural distortions that are induced by the oxygen vacancies.
Migration energies of oxygen vacancies were calculated explicitly for eight different configurations, where the zirconia positions neighboring the initial and final vacancy position during its jump were doped with yttrium ions. The migration energy for undoped BaZrO 3 is 0.66 eV, while those for the doped configurations vary between 0.46 and 1.16 eV. These migration energies were used in a pair-interaction-based KMC model to simulate the mobilities of oxygen vacancies and the oxygen ion conductivities.
A detailed comparison of our present results for the oxygen ion conductivity in Y-doped BaZrO 3 with our previously published results on proton conductivity shows the impact of nanoscale percolation of Y ions on both ionic conductivities: Both mobile defects, protons and oxygen vacancies, get trapped by single yttrium dopants and also in nanoscale configurations of Y ions. However unlike protons, oxygen vacancies are not accelerated in these trapping zones. Instead, Y ions block their way by means of high migration energies. In summary, configurations of Y ions with overlapping trapping zones form nanoscale   [9] www.advancedsciencenews.com www.advenergysustres.com percolation pathways with fast proton transport, while they block the transport of oxygen vacancies. Based on the insights from the above comparison, we propose a descriptor of the impact of a dopant on the oxygen ion or proton conductivity of acceptor-doped BaZrO 3 . This descriptor consists of only two migration energies of oxygen vacancies or protons: 1) in undoped BaZrO 3 ; and 2) "around" the dopant, that is, for the 1NN!1NN jump. If the second value is smaller than the first one, the dopant is a good candidate for improved defect mobility, that is larger than the defect mobility at infinite dilution, once the dopant fraction is high enough for the onset of nanoscale percolation of dopant ions. We believe that this simple descriptor will work for oxygen ions and for protons as well and will allow an easy screening of various new dopants.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.