Incorporation of protons and hydroxide species in BaZrO₃ and BaCeO₃

Abstract

Barium zirconate (BaZrO₃ or BZO) and barium cerate (BaCeO₃ or BCO) are among the best-performing proton-conducting oxides used as electrolytes in all-solid-state fuel and/or electrolysis cells. During synthesis, they are seeded with oxygen vacancies (V²⁺_O), which charge-compensate with acceptor dopants such as yttrium (Y⁻_Zr) and, upon exposure to water vapor, are replaced by interstitial protons (H⁺_i). Here, we investigate this and alternative processes for protonation by calculating defect formation energies, concentrations, and migration barriers for several relevant species, including H⁺_i, V²⁺_O, interstitial oxygen (O²⁻_i), and interstitial hydroxide (OH⁻_i), using density functional theory. We confirm that V²⁺_O are favorable under typical operating conditions, although at lower partial pressures of H₂ gas and wet conditions, H⁺_i becomes the dominant donor species. Higher H⁺_i concentrations in BCO than in BZO under comparable conditions help to explain the higher conductivity measured in BCO. OH⁻_i species are present in low concentrations in the bulk (particularly in BZO; they may incorporate in BCO under wet conditions), and their migration is slow; however, they may form at surfaces and help seed materials with H⁺_i. Alloying BZO and BCO improves ionic conduction in general, although the presence of native defects tends to impede kinetics. Our results show that high ionic conductivity can be achieved through optimizing synthesis conditions to maximize the concentrations of H⁺_i, as well as reducing defect-rich regions such as grain boundaries.

---

1. Introduction

Proton-conducting oxides (PCOs) are widely studied as solid-state electrolytes in ceramic fuel and/or electrolysis cells, on account of their favorable proton kinetics and ease of hydrogen incorporation.^1–5 Two ABO₃ perovskite oxides, barium zirconate (BaZrO₃ or BZO) and barium cerate (BaCeO₃ or BCO), are among the most attractive of these materials, and they are often alloyed together to merge BZO's superior chemical stability with BCO's higher proton conductivity.^6–11

It is generally accepted that acceptor doping (e.g., with yttrium, Y⁻_Zr/Ce) is necessary to incorporate protons (H⁺_i) into PCOs, as they initially seed the materials with oxygen vacancies (V²⁺_O), which serve as intermediaries for protonation upon water exposure:¹

2Y⁻_Zr/Ce + V²⁺_O + H₂O → 2Y⁻_Zr/Ce + 2H⁺_i. (1)

In light of their shared status as electron donors, some have speculated that V²⁺_O may actually compete with the goal of high H⁺_i concentrations.¹² Indeed, defect formation energy calculations show that H⁺_i has lower formation energies than V²⁺_O,¹³ suggesting that the introduction of V²⁺_O may limit the achievable proton concentration.

It is therefore of interest to consider other routes for protonation that do not require oxygen vacancies. Direct exposure to dry H₂ during synthesis with an acceptor species is one option:

BaO(s) + {Zr/Ce}O₂(s) + H₂+ 2Y(s) → Ba{Zr/Ce}O₃ + 2H⁺_i+ 2Y⁻_Zr/Ce. (2)

However, the high partial pressures of H₂ necessary for this reaction may not be compatible with traditional synthesis approaches. Sol–gel synthesis, which requires organic chelating agents, is performed in air in order to remove carbon and nitrogen impurities. Conversely, solid-state synthesis requires high temperatures on the order of 1450 °C, and it is difficult to operate a furnace filled with H₂ at such temperatures.

Another alternative is that exposure to water could give rise to H⁺_i in addition to an interstitial hydroxide (OH⁻_i) via the following reaction:

H₂O → OH⁻_i + H⁺_i. (3)

In this case, a dopant is not necessary to achieve charge balance. This reaction is relevant in materials with intrinsically oxygen deficient regions, such as hexagonal perovskites, which can absorb water without the need for doping.¹⁴

Whether or not these pathways are relevant is unclear. For one thing, the formation energy of OH⁻_i in the bulk has not previously been calculated, to our knowledge. However, previous computational studies have suggested that hydroxide migration may be relevant for protonation, particularly when moving from the surface to the bulk. Polfus et al. found adsorbed OH⁻_i to be favorable on surfaces of BZO, even more so than adsorbed O²⁻_i.¹⁵ Jing et al. calculated hydroxide migration barriers to be one-third as large as those for protons at the surface of SrZrO₃, another perovskite proton conductor, which suggests that hydroxide migration is key for protons to migrate into the bulk.¹⁶ Halwidl et al. also showed that interstitial hydroxide species created through the dissociation of water are very mobile at the surface of perovskite-like Sr₂RuO₄ and Sr₃Ru₂O₇.¹⁷

In addition to considering migration of other proton-related species, it is important to consider the effect of local chemical environments on proton diffusion. Grain boundaries, for instance, exhibit environments characterized by different concentrations of defects (such as vacancies) than in the bulk.^18–22 Several studies have shown that ions are less mobile near grain boundaries than in the bulk, meaning that grain boundaries limit the overall proton conductivity in polycrystalline samples of BZO and BCO.^23–30 Furthermore, the importance of BZO–BCO alloys means that ionic migration should also be considered in the presence of alloy impurities, i.e., Ce on a Zr site in BZO (Ce_Zr) or Zr on a Ce site in BCO (Zr_Ce).

In this work, we use density functional theory (DFT) calculations to study the properties of protons, oxide ions, and hydroxide ions in BZO and BCO, as well as their mobility in the presence of cation vacancies and alloy impurities. To begin, we calculate defect formation energies for each of the relevant species, using a hybrid exchange–correlation functional to capture electronic properties accurately. Using these results, we compute concentrations of defects and H-related impurities in Y-doped BCO and BZO. V²⁺_O will be the most prevalent native point defect; however, H⁺_i can replace it very easily, particularly under wet conditions. Their presence reduces the concentration of free charge carriers, which implies that proton-rich samples will exhibit less electrical leakage. OH⁻_i species will not typically be present in large concentrations, although they will be slightly more common in BCO than in BZO.

Subsequently, we calculate migration barriers for these defects and impurities. We calculate barriers for H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i in the bulk, which we then compare with values calculated in the presence of an alloy impurity (i.e., Ce_Zr in BZO or Zr_Ce in BCO) or an intrinsic vacancy. Comparing these values gives us a sense for how migration pathways will be affected by alloying and the presence of defect-rich regions such as grain boundaries. Unsurprisingly, H⁺_i has the fastest kinetics, though other species will also be fairly mobile. Alloying generally reduces migration barriers, while the presence of cation vacancies increases barriers due to Coulombic binding. As such, defect-rich regions like grain boundaries should be limited in order to optimize ionic conductivity.

2. Computational methods

We perform DFT calculations^31,32 using the Vienna ab initio simulation package (VASP).³³ We use the hybrid exchange–correlation functional of Heyd, Scuseria, and Ernzerhof (HSE06)³⁴ with 25% mixing of exact exchange in order to obtain accurate results for defect formation energies. For the calculation of migration barriers, we use the nudged elastic band (NEB) method with climbing images.^35,36 The generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) is used for NEB calculations for computational efficiency.³⁷ We apply projector augmented wave (PAW) potentials^38,39 with a plane-wave cutoff of 500 eV. The Ba 5s² 5p⁶ 6s², Zr 4s² 4p⁶ 5s² 4d², Ce 5s² 5p⁶ 6s² 5d¹ 4f¹, Y 4s² 4p⁶ 5s² 4d¹, and O 2s² 2p⁴ electrons are treated explicitly as valence. Spin polarization is included for all of our calculations. Supercells are constructed for defect formation and migration calculations; these consist of eight unit cells of BCO (2a × 2b × 2c in terms of the lattice vectors, 160 atoms) and 27 unit cells of BZO (3a × 3b × 3c, 135 atoms). The bulk properties we calculate for unit cells are described in our previous work.⁴⁰ For supercell calculations, a single special k point is used.

As a word of caution, we note that while our supercells are appropriately sized for defect calculations, they may be insufficient for studying global defect migration, particularly in highly disordered materials. As such, our results for defect migration provide relative trends, but our calculated migration barriers may not fully explain experimental measurements for these systems. Furthermore, as demonstrated in several review articles, experimental measurements of proton conductivity show significant variability, even in certain cases for materials with ostensibly the same composition.^3,41 As such, it is challenging to match computed barriers with experimental results.

2.1 Defect formation

We calculate the formation energy E^f(D^q) of a defect D in charge state q as:⁴²E(D^q) is the total energy of a supercell containing D^q; E_bulk is the total energy of the pristine supercell, containing no defects; n_A is the number of atoms of species A added (n_A < 0) or removed (n_A > 0) from the pristine system to create D^q; μ_A is the chemical potential of A; E_F is the position of the Fermi level; and Δ_corr is a finite-size correction term for charged defects.^43,44 The formation energy is exponentially related to the defect concentration c as:where N_sites is the site concentration of the defect, and k_B is Boltzmann's constant. As a result, lower formation energies correspond to exponentially higher defect concentrations. E_F is treated as an independent variable in calculating formation energies, meaning that plots of E^fvs. E_F will show each defect as a collection of lines with slopes q; where the slope changes, the preferred charge state changes as well.

For defect complexes, we can use formation energies to calculate a binding energy for a complex (AB) relative to its constituent, isolated defects (A and B):

E_bind(AB) = E^f(A) + E^f(B) − E^f(AB). (6)

A positive value of E_bind indicates that the complex is stable relative to the individual defects.

The chemical potential μ_A reflect the energetic preference for a specific element present in the system. We define it in terms of a deviation Δμ_A from a reference energy as:

μ_A = E_A + Δμ_A, (7)

where E_A is the total energies of the elemental reference structure, i.e., the ground-state structures of the Ba, Zr or Ce metals or an O atom in molecular O₂. In order to prevent formation of the elemental phases, we require that each Δμ_A ≤ 0. We assume conditions corresponding to thermodynamic equilibrium, which is captured by the following expression for the case of Ba{Zr,Ce}O₃:

Δμ_Ba + Δμ_{Zr,Ce} + 3Δμ_O = ΔH^f(Ba{Zr,Ce}O₃), (8)

where ΔH^f(Ba{Zr,Ce}O₃) is the enthalpy of formation of BZO or BCO. These and other enthalpies of formation have been calculated previously using a similar computational approach.^13,45

Additional bounds are placed on Δμ_Ba and Δμ_{Zr,Ce} to ensure that limiting phases such as BaO, ZrO₂, and CeO₂ do not precipitate. For Ba, this condition is expressed as:

Δμ_Ba + Δμ_O ≤ ΔH^f(BaO). (9)

And for {Zr,Ce}:

Δμ_{Zr,Ce} + 2Δμ_O ≤ ΔH^f({Zr,Ce}O₂). (10)

Based on eqn (8) and the upper limit of eqn (9), we can define “Ba-rich” conditions (equivalently, “Zr-poor”) as those for which Δμ_Ba is maximized while preventing precipitation of BaO or Ba metal. Similarly, eqn (8) and (10) allow us to define “{Zr,Ce}-rich” (“Ba-poor”) conditions, where Δμ_{Zr,Ce} is maximized while ensuring that {Zr,Ce}O₂ does not precipitate.

We treat impurity chemical potentials (e.g., for H and Y) in a similar fashion. As for eqn (9) and (10), we establish an upper limit on Δμ by ensuring that a limiting phase does not precipitate. For yttrium, that phase is Y₂O₃, leading to the condition:

2Δμ_Y + 3Δμ_O ≤ ΔH^f(Y₂O₃). (11)

For hydrogen, we choose H₂O as our limiting phase, as it is most relevant experimentally, although Ba(OH)₂ actually provides a more restrictive limit.¹³ The limiting condition is expressed as:

2Δμ_H + Δμ_O = Δμ_H2O ≤ ΔG^f(H₂O), (12)

where we introduce Δμ_H2O as the chemical potential of water vapor. By fixing Δμ_H to its upper limit for a given value of Δμ_O, we capture the “H-rich” limit that is likely to be most preferred for protonating BZO and BCO.

Note that ΔG^f(H₂O) is the upper limit here, rather than ΔH^f as with our other limiting phases. Strictly speaking, the free energies ΔG^f are the proper limits in each case and are reflective of conditions at finite temperatures. However, the entropic contribution is most pronounced for gaseous species such as water vapor, so while we can approximate the upper limits in eqn (8)–(11) with DFT-calculated ΔH^f values, we necessarily need to consider the temperature-dependent ΔG^f in this case. To that end, we use tabulated values for ΔG^f(T).⁴⁶

Finally, for elements with a gas-phase reference state, such as oxygen and hydrogen (and also water vapor), the Δμ values can be connected to experimentally measurable partial pressures using the expression:

where H⁰(T) and S⁰(T) are tabulated for O₂, H₂, and H₂O;⁴⁶P_A is the partial pressure; and P⁰ is the pressure at the standard conditions used in the tabulation.

2.2 Defect concentrations

Using eqn (5), we can identify how defect concentrations change under different operating conditions. Specifically, by enforcing charge neutrality among all the defects we calculate, we can determine the concentrations for each species and the position of the Fermi level for certain conditions. This procedure requires two additional parameters, namely, the electron and hole concentrations, which we determine by integrating the calculated densities of states (DOS). We determine electron concentrations by integrating the calculated DOS near the CBM, using the expression:⁴⁷with g_C(E) being the conduction band DOS and f(E) being the Fermi–Dirac occupation function. Similarly, for holes, we determine the carrier concentration by integrating the DOS near the VBM via:with g_V(E) being the valence band DOS.

3. Results

3.1 Defect formation

3.1.1 Defect formation energies. We present results for the formation energies of relevant H- and O-related species (and Y⁻_{Zr,Ce} acceptors) in BZO and BCO in Fig. 1, assuming the Ba-rich limit and Y solubility limit, as discussed previously. As expected, V²⁺_O and H⁺_i are the lowest-energy electron donors for most of the bandgap, and they are particularly low in energy for Fermi level positions ∼1-to-2 eV above the VBM. Higher energy species include substitutional H⁺_O, and interstitials O_i and OH⁻_i. O_i preferentially adopts a split interstitial configuration, in agreement with a recent study,⁴⁸ while OH⁻_i resembles O²⁻_i with an additional H⁺_i attached to one of the oxygen atoms. It follows that we can treat OH⁻_i as a complex of O²⁻_i and H⁺_i and calculate its binding energy: 1.42 eV in BZO and 1.18 eV in BCO. Both of these binding energies indicate that OH⁻_i is highly stable with respect to the individual defects and is therefore likely to form if both are simultaneously present. In the absence of dopants, the formation energies of H⁺_i and other relevant species will be high; thus, a compensating acceptor species such as Y⁻_Zr/Ce will still be necessary.

Fig. 1 Formation energies of hydrogen- and oxygen-related defects and impurities in (a) BaZrO₃ and (b) BaCeO₃. Δμ_O is set at −2.42 eV, which corresponds to sintering conditions at 1650 °C,⁴⁹ and equilibrium with H₂O is assumed. Ba-rich conditions and the Y solubility limit are assumed for calculating the formation energies of Y_Zr and Y_Ce.

Our defect formation energies allow us to comment on the relative energetics of eqn (1) and (3), which describe possible routes toward protonation in BZO and BCO, by comparing the formation energies of products and reactants. For eqn (1), the enthalpy of the reaction will be related to 2E^f[H⁺_i] − E^f[V²⁺_O], while for eqn (3), the enthalpy is related to E^f[H⁺_i] + E^f[OH⁻_i]. The difference between these two quantities (2E^f[H⁺_i] − E^f[V²⁺_O] and E^f[H⁺_i] + E^f[OH⁻_i]) provides an estimate of the relative favorability of the two reactions; notably, this quantity is independent both of chemical potentials and the Fermi level position. In BZO, we find eqn (1) to be 4.14 eV lower in energy than eqn (3), while for BCO, we find eqn (1) to be 3.02 eV lower in energy. In both cases, these relative energies suggest that the widely proposed V²⁺_O-mediated reaction of eqn (1) will proceed much more favorably, although the process described by eqn (3), whereby V²⁺_O are not required and OH⁻_i species are created, will be more competitive in BCO. Most likely, this result extends qualitatively to surfaces and interfaces, where OH⁻_i are most likely to form initially.

3.1.2 Defect concentrations. Next, we calculate defect concentrations in bulk BZO and BCO under varying chemical potential conditions. We assume a fixed, typical Y doping level of 20 at%, as well as a temperature of 900 K, which is close to typical operating temperatures for BZO and BCO fuel and electrolysis cells.^4,9 To account for finite temperature effects, we include a harmonic correction to the formation energies of V_O and H⁺_i based on previous work;^50,51 these corrections increase the formation energies of V²⁺_O and H⁺_i by 0.13 eV and 0.21 eV, respectively.

In each of our concentration plots, we maintain a fixed value of Δμ_H and vary Δμ_O from −5 eV to 0 eV, which in turn yields a range of Δμ_H2O that we use as our x-axis. Using eqn (13), we can connect Δμ_H2O to partial pressures of water (P_H2O), which we show on the upper x-axis of our concentration plots. We consider two choices of P_H2 : 1 atm (corresponding to Δμ_H = −0.625 eV) and 10⁻⁵ atm (Δμ_H = −1.071 eV), to compare H-rich and H-poor conditions. Otherwise, we use intermediate chemical potentials halfway between the Ba-rich and Ba-poor limits. An upper limit on our x-axis is set by ΔG^f(H₂O), which at 900 K is −2.05 eV, although we extend slightly beyond this limit for illustrative purposes in the case of P_H2 = 1 atm. We include certain defect formation energies calculated previously^13,40,45,52 to generate these plots.

To start, we plot concentrations of defects and impurities in BZO in Fig. 2(a) and (b). In panel (a), for which P_H2 = 1 atm, we identify three distinct regimes based on which donor species compensates with Y⁻_Zr. At dry conditions, below approximately P_H2O = 10⁻²⁰ atm, the substitutional hydride H⁺_O will be the dominant donor. On the other hand, for wet conditions of P_H2O ≥ 10⁻¹⁰ atm, interstitial protons H⁺_i will dominate. In between those two extremes, V²⁺_O will have the highest concentration among donor species. The situation is similar in Fig. 2(b), for P_H2 = 10⁻⁵ atm, with the main difference being that H⁺_O concentrations will be significantly less competitive with V²⁺_O for dry conditions.

Fig. 2 Defect, impurity, and free carrier concentrations at 900 K as a function of Δμ_H2O in BaZrO₃ under (a) P_H2 = 1 atm, and (b) P_H2 = 10⁻⁵ atm; and for BaCeO₃ under (c) P_H2 = 1 atm, and (d) P_H2 = 10⁻⁵ atm. The Ba and Zr chemical potentials correspond to intermediate conditions between the Ba-rich and Ba-poor limits.

We note two more significant results for BZO. First, the Fermi level positions will be higher than those we calculated in previous work for BZO absent any hydrogen exposure,⁴⁰ spanning a range of about 1.6–3.1 eV above the VBM for P_H2 = 1 atm [Fig. 2(a)], and about 1.2–2.4 eV above the VBM for P_H2 = 10⁻⁵ atm [Fig. 2(b)]. Under the more H-poor conditions in panel (b), free hole (h⁺) concentrations are noticeably higher than under the H-rich conditions in panel (a); however, they plateau just above 10¹⁴ cm⁻³. Free electron (e⁻) concentrations do not even appear in the plotted range in either case. In general, these results suggest that the presence of H⁺_i will minimize harmful electrical leakage in these materials by limiting the presence of free carriers and/or polarons. This result is consistent with a previous report that proton conductivity will significantly exceed electrical conductivity in BZO at high P_H2O.⁵³ Second, we find that the hydroxide species, OH⁻_i, does not appear on our concentration plots; as such, we do not expect measurable concentrations of OH⁻_i in bulk BZO.

Next, in Fig. 2(c) and (d) we plot defect concentrations in bulk BCO. The trends are generally similar to those we observed in BZO. For P_H2 = 1 atm, shown in panel (c), the dominant compensating donor species varies from H⁺_O under dry conditions (P_H2O ≤ 10⁻²⁰ atm), to V²⁺_O for intermediate conditions (10⁻¹⁸ atm ≤ P_H2O ≤ 10⁻¹² atm), to H⁺_i under wet conditions (P_H2O ≥ 10⁻¹⁰ atm). In comparison with BZO, BCO has higher concentrations of cation antisites (Ce_Ba²⁺) and wrong-site Y donors (Y_Ba⁺). Ce_Ba²⁺ in particular will compete with V²⁺_O, with a more pronounced effect under more Ba-poor conditions than those shown here. In addition, OH⁻_i concentrations are noticeably higher than in BZO, particularly under wet conditions, where they approach concentrations of 10¹⁸ cm⁻³. For P_H2 = 10⁻⁵ atm, shown in panel (d), V²⁺_O is the dominant donor for P_H2O ≤ 10⁻¹² atm (again accompanied by Ce_Ba²⁺), and H⁺_i takes its place for wetter conditions. These H-poor conditions are less amenable to OH⁻_i formation; instead, O⁰_i concentrations will become increasingly sizeable with increasing Δμ_H2O.

Once again, as with BZO, we observe relatively high Fermi levels, ranging from 1.9–3.4 eV above the VBM for P_H2 = 1 atm, and 1.4–2.4 eV for P_H2 = 10⁻⁵ atm. Correspondingly, h⁺ concentrations are significantly lower than those we found in our previous study,⁴⁰ although e⁻ concentrations can be high under extremely dry/H₂-rich conditions.

The results of Fig. 2 suggest a possible explanation for BCO's superiority to BZO as a proton conductor. Specifically, the concentration of H⁺_i is higher in BCO than in BZO under comparable conditions, and H⁺_i is preferred over V²⁺_O over a wider range of chemical potentials in BCO. We demonstrate this finding schematically in Fig. 3, which plots the concentrations of H⁺_i and V²⁺_O side-by-side for analogous conditions in BZO and BCO. Specifically, we choose conditions that are typical of the two electrodes in an electrolysis cell under open-circuit conditions and change the chemical potentials linearly between these two extremes. We use H-rich conditions at the H₂ electrode, with P_H2 = 1 atm and P_O2 = 7 × 10⁻²⁰ atm, and O-rich conditions at the O₂/H₂O electrode, with P_O2= 0.7 atm and P_H2O = 0.3 atm. As our results indicate, H⁺_i concentrations are dominant over most of the electrolyte, with V²⁺_O only beginning to displace H⁺_i very close to the H₂ electrode. However, this crossover point between H⁺_i and V²⁺_O occurs more quickly for BZO than for BCO, implying a lower proton conductivity due to limited H⁺_i concentrations in BZO. Tuning the precise conditions at the two electrodes may be a viable strategy to keep H⁺_i concentrations high.

Fig. 3 Concentrations of H⁺_i and V²⁺_O in BaCeO₃ (solid lines) and BaZrO₃ (dashed lines) under a range of chemical potential conditions chosen to model operating conditions in an electrolysis cell. At the hydrogen electrode, we use P_H2 = 1 atm and P_O2 = 7 × 10⁻²⁰ atm; at the oxygen electrode, we use P_H2O = 0.3 atm and P_O2 = 0.7 atm. A temperature of 900 K is used.

Finally, we note that the concentrations we plot in Fig. 2 and 3 evolve slightly with increasing temperature. While our assumed temperature of 900 K is typical for electrolysis applications, other use cases for these materials (such as in oxide ion conduction) require higher temperatures. The most noticeable impact of higher temperatures is to increase the free carrier concentration, which can increase the degree of electrical leakage, as has been observed experimentally.^54,55 This rise in carrier concentrations is accompanied by a decrease in H⁺_i concentrations and an increase in V²⁺_O concentrations for analogous conditions. Increasing temperature appears to depress OH⁻_i concentrations, however. Thus, to the extent that water can be directly converted into H⁺_i and OH⁻_iviaeqn (3), we expect that process to be most favorable at the lower temperatures most characteristic of electrolysis.

3.2 Defect migration

3.2.1 Migration in bulk regions. Next, we use the NEB method to calculate the bulk migration barriers for H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i. Ranges encompassing our computed barriers are listed in Table 1; these capture the anisotropy present in many of the pathways, particularly those in BCO. In our present discussion, we do not explicitly consider the interactions between mobile species and dopants. The Coulombic interactions between Y⁻_Zr and H⁺_i, for example, are non-negligible,¹³ and they certainly influence defect migration in actual devices. However, at present, we are mostly concerned with the relative impacts of defects and alloy impurities on mobility, and as such, we neglect dopant interactions. For more accurate treatments of proton–dopant interactions in BZO, we refer the reader to other computational studies focused specifically on this problem, which employ larger supercells to capture varied dopant configurations.^56,57

Table 1 Calculated migration barriers for H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i in bulk BaZrO₃ and BaCeO₃, as well as in the presence of alloy impurities (Ce_Zr in BaZrO₃, BZO:Ce; and Zr_Ce in BaCeO₃, BCO:Zr) and vacancies (V²⁻_Ba, V⁴⁻_{Zr,Ce}, and V²⁺_O). The total charge state in the simulations accounts for the preferred charge states of both the vacancies and the mobile species

		BaZrO3 (eV)
Species	Pathway	Bulk	BZO:Ce	V^2−Ba	V^4−Zr	V^2+O
a This barrier is simply the change in energy between the initial (O^2−i separated from V^4−Zr) and final (O^2−i near V^4−Zr) states. The reverse process proceeds spontaneously.
H^+i	Hop	0.21	0.13–0.36	0.64	0.43	0.22–0.23
	Rotate	0.18	0.10–0.28	0.21–0.93	0.06–1.89	0.22–0.28
V^2+O	Hop	0.72	0.42–0.96	0.59–1.00	0.75–2.25	0.21–0.99
O^2−i	Hop	0.58	0.16–0.55	0.48	1.37^a	—
	Rotate	1.28	0.78–1.23	0.78–0.95	0.33–1.02	—
OH^−i	Hop	0.55	0.24–0.85	0.57–1.69	0.85	—
	Rotate	1.67	1.54–1.68	0.24–2.17	0.61–1.52	—

---

		BaCeO3 (eV)
Species	Pathway	Bulk	BCO:Zr	V^2−Ba	VCe^4−	V^2+O
H^+i	Hop	0.21–0.39	0.13–0.55	0.26–1.00	0.21–1.57	0.21–0.38
	Rotate	0.07–0.14	0.06–0.22	0.11–0.58	0.13–1.60	0.12–0.20
V^2+O	Hop	0.45–0.50	0.32–0.81	0.30–1.08	0.30–1.51	0.05–0.72
O^2−i	Hop	0.22–0.27	0.20–0.64	0.30–0.50	0.85–0.93	—
	Rotate	0.58–0.64	0.53–0.87	0.66–0.74	0.21–0.22	—
OH^−i	Hop	0.37–0.46	0.09–0.73	0.24–1.10	0.78–0.84	—
	Rotate	0.80–0.95	0.72–1.10	0.53–1.31	0.34–0.98	—

---

To begin, we calculate barriers for H⁺_i migration. We consider two pathways, summarized visually in Fig. 4(a) and (b): hopping, whereby H⁺_i moves from its host oxygen atom to another nearby oxygen, and rotation, where H⁺_i simply changes its orientation while remaining attached to the same oxygen atom. Our results compare well to those of previous calculations on proton migration in BZO⁵⁸ and BCO,⁵⁹ exclusively considering regions far from dopants like Y⁻_Zr where additional interactions can be neglected. We compute nearly identical barriers for proton hopping in both materials, although due to its lower symmetry, BCO has more possible pathways, some of which are slightly higher in energy.

Fig. 4 Hopping and rotation migration pathways for hydrogen and oxygen diffusion in BaZrO₃ and BaCeO₃. Migration pathways of interstitial protons (H⁺_i) in (a) bulk BaZrO₃ and (b) bulk BaCeO₃, and hopping migration pathways of oxygen vacancies (V²⁺_O) in (c) bulk BaZrO₃ and (d) bulk BaCeO₃.

Next, we calculate barriers for V²⁺_O migration, which proceeds via a hopping mechanism as shown in Fig. 4(c) and (d). For a reaction such as that described by eqn (1) to drive protonation, V²⁺_O must necessarily be mobile in order to replace vacancies filled by water at the surface. Oxygen diffusion has been observed in BZO and BCO,^60–63 particularly at elevated temperatures.^64,65 We calculate a barrier of 0.72 eV in BZO, similar to the value of 0.66 eV recently calculated in another study.⁶⁶ For BCO, we calculate barriers ranging from 0.45–0.50 eV. We are unaware of previous NEB calculations on V²⁺_O migration in BCO, and experimental activation energies are higher (0.69–0.90 eV ^1,67,68). However, when adding our calculated migration barriers to the energy contribution from thermal activation of the mobile species (namely, their formation energies), our results are consistent with these values.

We also consider the migration of O²⁻_i, despite the higher formation energies of these species, as it is closely related to OH⁻_i. We identify pathways for hopping and rotation, shown for BZO in Fig. 5(a), and for BCO in Fig. 5(b). In the hopping pathway, migration proceeds via an interstitialcy (“kick-out”) process, meaning that the migrating species replaces a lattice O atom after “kicking it out” of its lattice site. Hopping is favored over rotation in both materials; in fact, in BCO, its barrier is roughly equal to that of H⁺_i hopping. These hopping barriers are lower than those of V²⁺_O; nonetheless, the markedly lower formation energy of V²⁺_O will ensure that any oxygen conductivity will be predominately vacancy mediated. We attribute the larger barriers for rotation to the mobile O atoms moving too close to large Ba atoms.

Fig. 5 Migration pathways for interstitial O-related species in BaZrO₃ and BaCeO₃. Hopping and rotation for O²⁻_i in (a) BaZrO₃, and (b) BaCeO₃. Hopping and rotation for OH⁻_i in (c) BaZrO₃, and (d) BaCeO₃. Mobile O and H atoms are highlighted in both the initial and final states for clarity, with yellow arrows used to indicate the directions of their motion.

Finally, we consider OH⁻_i migration, using the same approximate pathways as O²⁻_i, as shown in Fig. 5(c) for BZO and in Fig. 5(d) for BCO. As for O²⁻_i, hopping proceeds via an interstitialcy process; however, there is an additional step of rotation and transfer of hydrogen to complete the movement. As such, OH⁻_i migration can be considered to be a composite of H⁺_i and O²⁻_i migration. Hopping pathways for OH⁻_i generally proceed more favorably than for O²⁻_i, while rotation is less favorable, again due to the close proximity of Ba cations at the saddle point configuration. These barriers suggest that OH⁻_i will be mobile if it is present, and due to its large binding energy, it will be favored over O²⁻_i if hydrogen is present. OH⁻_i migration also proceeds more favorably in BCO than in BZO; as we showed earlier, OH⁻_i will also have much higher concentrations in BCO.

3.2.2 Migration in alloyed systems. As mentioned in the introduction, BZO and BCO are often alloyed together in devices. Therefore, we examine the effect of alloying on the migration barriers previously discussed by adding a single impurity atom (i.e., Zr_Ce in BCO, which we will refer to as BCO:Zr; or similarly, Ce_Zr in BZO, BZO:Ce) in the vicinity of the migrating species. We sample several possible sites for the impurity atom and compute the resultant migration barriers. In this way, we seek to determine whether alloying leads to a discernible trend in the migration of O- and H-related species. Given that our results capture the dilute limit of alloying, they are most representative of alloys at the compositional extremes, i.e., Ce-rich or Zr-rich. Of these, Ce-rich alloys in particular have attracted considerable attention.^4,9,11

In Table 1, we list the range of barriers computed for the hopping and rotation pathways discussed above. In each case, certain placements of the impurity reduce the energetic barrier for migration. This improvement is generally most pronounced in BZO:Ce, which follows from the lower barriers we identify in bulk BCO; however, BCO:Zr can also exhibit slightly improved kinetics for ion migration.

However, certain pathways in alloyed systems also have larger barriers than in the bulk materials for almost every mobile species we consider. This mixed behavior likely arises from alloy impurities breaking the local symmetry in such a way as to create additional space for certain migration pathways, while at the same time constricting other pathways. As such, it is not entirely clear what the overall impact of alloying will be on large-scale ion migration. Presumably, in highly heterogeneous alloys, species will be able to migrate along more energetically favorable pathways, whereas in a phase-separated system, migration would be slower in Zr-rich (BZO-like) than in Ce-rich (BCO-like) regions. It is also possible that alloying may create higher-energy corridors that serve as bottlenecks for ion transport; however, a detailed examination of this possibility is beyond the scope of this work.

In Table 2, we list the binding energies (determined using eqn (6)) for complexes involving a mobile defect and an alloy impurity. These binding energies represent a coulombic barrier that must be overcome—in addition to the migration barrier—for a mobile species to break away from the complex. The values we list reflect the most energetically favored configurations of the complexes, which suggests a worst-case-scenario for migration. It is evident that alloy species exhibit some amount of Coulombic binding with V²⁺_O and O²⁻_i, thereby limiting their mobility, while the negative binding energies for complexes between H⁺_i and alloy species suggest that they will not inhibit proton mobility. The effect on OH⁻_i is more complicated, as our results suggest that Ce-alloying in BZO will bind OH⁻_i, while Zr-alloying in BCO will repel OH⁻_i. Of course, it is unlikely that alloy impurities will be present at the dilute limit, as they are in our simulation cells; for a near equal mix of Zr and Ce atoms, our results suggest very little net impact on OH⁻_i migration.

Table 2 Calculated binding energies for H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i in complexes with alloy impurities and vacancies, in BaZrO₃ and BaCeO₃. The total charge states account for the preferred charge states of constituent species (q = −4 for V⁴⁻_{Zr,Ce}, q = −2 for V²⁻_Ba and O²⁻_i, q = −1 for OH⁻_i, q = 0 for Ce_Zr and Zr_Ce, q = +1 for H⁺_i, and q = +2 for V²⁺_O)

Mobile species	Binding energy (eV)
	Material	CeZr/ZrCe	V^2−Ba	V^4−Zr/Ce	V^2+O
H^+i	BaZrO3	−0.09	0.77	1.53	−0.09
	BaCeO3	−0.02	0.75	1.29	−0.18
V^2+O	BaZrO3	0.28	0.36	1.74	−0.10
	BaCeO3	0.20	0.69	1.38	0.06
O^2−i	BaZrO3	0.43	−0.48	−0.45	—
	BaCeO3	0.21	−0.32	−0.39	—
OH^−i	BaZrO3	0.52	0.91	0.33	—
	BaCeO3	−0.29	0.33	0.16	—

---

3.2.3 Migration in defect-rich regions. Local atomic disorder is also highly relevant to transport in PCOs. Here, we consider a particularly small scale of disorder, namely, the presence of vacancies near mobile species. Previously, it has been shown that cation vacancies (V_Ba and V_{Zr,Ce}) are energetically unfavorable in bulk BZO and BCO;^13,52 however, they may be more prevalent in disordered regions, such as space charge layers around grain boundaries. Grain boundary cores exhibit excess positive charge,^69,70 which is accompanied by the depletion of positive charge in neighboring space charge layers;^71,72 negatively charged space charge layers could be consistent with increased cation vacancy concentrations. V_Ba concentrations, in particular, may be noticeable in these materials.⁷³ The positive charge at the grain boundary core, on the other hand, has been attributed to V²⁺_O accumulation.^18,74,75 Thus, their interactions with protons in bulk and grain boundary regions are likely to be relevant to ionic transport as well.

Tables 1 and 2 list migration barriers and binding energies, respectively, for mobile species in proximity to native vacancies. For each mobile species, we include both its preferred charge state (+1 for H⁺_i, +2 for V²⁺_O, −2 for O²⁻_i, and −1 for OH⁻_i), as well as that of the nearby vacancies (−2 for V²⁻_Ba, −4 for V⁴⁻_Zr/Ce, and +2 for V²⁺_O). Note that for O²⁻_i and OH⁻_i, we do not include results in the vicinity of V²⁺_O, as the interstitial O atom in both cases will spontaneously fill the vacancy. Several of the combinations have multiple pathways that we consider; in these cases, we list a representative range of migration barriers. Binding energies are computed using the most energetically stable complexes among those we identify.

The presence of cation vacancies has a significant impact on proton migration, with V²⁻_Ba and V⁴⁻_{Zr,Ce} giving rise to larger H⁺_i migration barriers in general. Additionally, protons have large, positive binding energies with cation vacancies (particularly with V⁴⁻_Zr/Ce), which will severely hamper proton migration in their vicinity. Furthermore, if cation vacancies dominate the defect chemistry in disordered regions, as may occur in space charge layers due to dopant segregation to grain boundary cores,^18,76 the concentration of protons will drop precipitously due to an increase in their formation energy. If we simply approximate the concentration of H⁺_i by assuming equilibrium with V²⁻_Ba, the concentrations will decrease by as much as six orders of magnitude compared with the results presented in Fig. 2. These large barriers and binding energies, and decreased proton concentrations, may help to explain the observation of higher activation energies for proton migration in grain boundaries compared to bulk regions.³⁰ V²⁺_O migration is similarly affected by the presence of cation vacancies, with V⁴⁻_Zr/Ce again leading to larger migration barriers and binding energies than V²⁻_Ba.

Bringing a V²⁺_O defect into proximity with H⁺_i does not significantly affect the migration barriers. However, having two V²⁺_O defects in close proximity does appear to enhance V²⁺_O diffusion along certain directions, with barriers as small as 0.21 eV in BZO and 0.05 eV in BCO. As both V²⁺_O and H⁺_i are positively charged, Coulombic trapping is not a concern, as evidenced by the small (and in most cases negative) binding energies.

Cation vacancies can have either a beneficial or harmful impact on O²⁻_i and OH⁻_i migration, depending on the placement of the vacancy relative to the mobile species. Rotation, as we have discussed, is limited in bulk systems by the presence of large Ba cations; removing these cations will therefore make rotation significantly easier along pathways including the vacancy. Hopping barriers, on the other hand, are generally increased slightly by the presence of vacancies (particularly V⁴⁻_Zr/Ce).

Interestingly, O²⁻_i has negative binding energies in complexes with cation vacancies, while OH⁻_i has positive binding energies. This discrepancy reflects the nature of OH⁻_i as a complex between H⁺_i and O²⁻_i; as such, the binding energies with cation vacancies are generally intermediate between the binding energies for isolated H⁺_i and O²⁻_i. The stronger Coulombic attraction between cation vacancies and H⁺_i may cause dissociation of OH⁻_i, as H⁺_i will be drawn to the vacancies while O²⁻_i are repelled. Such an effect might be observable near grain boundaries, where we expect cation vacancies to be most prevalent.

4. Conclusion

In conclusion, we have investigated the formation energies, concentrations, and migration barriers of H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i in BZO and BCO. Our formation energies suggest that it will be difficult to dissociate water and protonate BZO and BCO without oxygen vacancies (slightly more so in BZO than in BCO), lending credence to the widely accepted notion that V²⁺_O are needed for protonation. When doped with an element like yttrium, three donor species can be stabilized based on the partial pressure of water vapor: H⁺_O at extremely dry conditions, V²⁺_O at intermediate conditions, and H⁺_i under wet conditions. Decreasing the partial pressure of H₂ makes V²⁺_O favored over H⁺_O at extremely dry conditions, while leaving H⁺_i mostly unaffected. The presence of H⁺_i limits the concentration of free carriers and polarons, which has advantages for avoiding electrical leakage. O²⁻_i and OH⁻_i will have very low concentrations in bulk regions, although OH⁻_i should have measurable concentrations in BCO under wet conditions. Considering both kinetics and thermodynamics of proton incorporation in BZO and BCO, our results indicate that the tendency for higher conductivities in BCO is likely attributed to the greater relative solubility of protons under the same conditions.

The precise effect of alloying on ionic mobility in BZO and BCO is highly dependent on the position of alloy impurity states relative to the mobile species, meaning that certain ionic pathways may be more favored in alloys, and the ordering of Ce and Zr atoms may impact kinetics. The presence of vacancies has a mixed effect on migration for H⁺_i, V²⁺_O, O²⁻_i, and OH⁻_i, which has implications for understanding diffusion in disordered or off-stoichiometric regions (e.g. surfaces, interfaces, and grain boundaries). Cation vacancies hinder H⁺_i and V²⁺_O conduction due to Coulombic binding, while additional V²⁺_O defects have little impact. For O²⁻_i and OH⁻_i migration, cation vacancies reduce energetic barriers for rotation if they lie in the path of motion; otherwise, they tend to increase migration barriers. In addition, OH⁻_i will be bound to cation vacancies, perhaps leading to the dissociation of OH⁻_i into H⁺_i and O²⁻_i. Overall, these results suggest that in order to maximize the conductivity of the most mobile species in these systems (H⁺_i and V²⁺_O), defect-rich regions such as grain boundaries should be avoided as much as possible.

Author contributions

Andrew J. E. Rowberg: formal analysis, investigation, methodology, visualization, writing – original draft, writing – review & editing. Meng Li: resources, writing – review & editing. Tadashi Ogitsu: funding acquisition, project administration, writing – review & editing. Joel B. Varley: conceptualization, methodology, supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge support from the HydroGEN Advanced Water Splitting Materials Consortium, established as part of the Energy Materials Network under the U.S. Department of Energy (DOE), the Office of Energy Efficiency and Renewable Energy (EERE), the Hydrogen and Fuel Cell Technologies Office (HFTO). Part of this work was performed under the auspices of the DOE by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. M. Li additionally acknowledges support from the DOE Idaho Operations Office under contract no. DE-AC07-05ID14517. The research was performed using computational resources sponsored by the DOE's EERE and located at the National Renewable Energy Laboratory.

Notes and references

1. K.-D. Kreuer, Annu. Rev. Mater. Res., 2003, 33, 333–359.
2. T. Norby, Solid State Ionics, 1999, 125, 1–11.
3. N. Kochetova, I. Animitsa, D. Medvedev, A. Demin and P. Tsiakaras, RSC Adv., 2016, 6, 73222–73268.
4. S. Rajendran, N. K. Thangavel, H. Ding, Y. Ding, D. Ding and L. M. Reddy Arava, ACS Appl. Mater. Interfaces, 2020, 12, 38275–38284.
5. C. Duan, J. Huang, N. Sullivan and R. O’Hayre, Appl. Phys. Rev., 2020, 7, 011314.
6. K. H. Ryu and S. M. Haile, Solid State Ionics, 1999, 125, 355.
7. K. Katahira, Y. Kohchi, T. Shimura and H. Iwahara, Solid State Ionics, 2000, 138, 91–98.
8. J. Lü, L. Wang, L. Fan, Y. Li, L. Dai and H. Guo, J. Rare Earths, 2008, 26, 505–510.
9. W. Bian, W. Wu, B. Wang, W. Tang, M. Zhou, C. Jin, H. Ding, W. Fan, Y. Dong, J. Li and D. Ding, Nature, 2022, 604, 479–485.
10. S. Choi, C. J. Kucharczyk, Y. Liang, X. Zhang, I. Takeuchi, H.-I. Ji and S. M. Haile, Nat. Energy, 2018, 3, 202–210.
11. L. Yang, S. Wang, K. Blinn, M. Liu, Z. Liu, Z. Cheng and M. Liu, Science, 2009, 326, 126–129.
12. T. Norby, M. Widerøe, R. Glöckner and Y. Larring, Dalton Trans., 2004, 3012–3018.
13. A. J. E. Rowberg, L. Weston and C. G. Van de Walle, ACS Appl. Ener. Mater., 2019, 2, 2611–2619.
14. S. Fop, K. S. McCombie, E. J. Wildman, J. M. Skakle and A. C. Mclaughlin, Chem. Commun., 2019, 55, 2127–2137.
15. J. M. Polfus, T. S. Bjørheim, T. Norby and R. Bredesen, J. Mater. Chem. A, 2016, 4, 7437–7444.
16. Y. Jing, H. Matsumoto and N. R. Aluru, Chem. Mater., 2018, 30, 138–144.
17. D. Halwidl, B. Stöger, W. Mayr-Schmölzer, J. Pavelec, D. Fobes, J. Peng, Z. Mao, G. S. Parkinson, M. Schmid, F. Mittendorfer, J. Redinger and U. Diebold, Nat. Mater., 2016, 15, 450–455.
18. J. M. Polfus, K. Toyoura, F. Oba, I. Tanaka and R. Haugsrud, Phys. Chem. Chem. Phys., 2012, 14, 12339–12346.
19. M. Vollmann, R. Hagenbeck and R. Waser, J. Am. Ceram. Soc., 1997, 80, 2301–2314.
20. R. A. De Souza, Z. A. Munir, S. Kim and M. Martin, Solid State Ionics, 2011, 196, 1–8.
21. G. L. Burton, S. Ricote, B. J. Foran, D. R. Diercks and B. P. Gorman, J. Am. Ceram. Soc., 2020, 103, 3217–3230.
22. S. Haile, G. Staneff and K. Ryu, J. Mater. Sci., 2001, 36, 1149–1160.
23. F. Iguchi, N. Sata, T. Tsurui and H. Yugami, Solid State Ionics, 2007, 178, 691–695.
24. J.-H. Yang, B.-K. Kim and Y.-C. Kim, Solid State Ionics, 2015, 279, 60–65.
25. B. Joakim Nyman, E. E. Helgee and G. Wahnström, Appl. Phys. Lett., 2012, 100, 4355.
26. A. Lindman, E. E. Helgee, B. Joakim Nyman and G. Wahnström, Solid State Ionics, 2013, 230, 27–31.
27. A. Lindman, E. E. Helgee and G. Wahnström, Chem. Mater., 2017, 29, 7931–7941.
28. S. Yue, Y. Jing, Y. Sun, J. Zhao and N. Aluru, Ceram. Int., 2022, 48, 2097–2104.
29. M. Shirpour, R. Merkle, C. Lin and J. Maier, Phys. Chem. Chem. Phys., 2012, 14, 730–740.
30. R. Pornprasertsuk, O. Kosasang, K. Somroop, S. Jinawath and F. B. Prinz, ECS Trans., 2010, 25, 367.
31. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
32. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
33. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
34. G. E. Heyd, J. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
35. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
36. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
39. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
40. A. J. E. Rowberg, M. Li, T. Ogitsu and J. B. Varley, Phys. Rev. Mater., 2023, 7, 015402.
41. C. Y. R. Vera, H. Ding, D. Peterson, W. T. Gibbons, M. Zhou and D. Ding, J. Phys.: Energy, 2021, 3, 032019.
42. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van de Walle, Rev. Mod. Phys., 2014, 86, 253.
43. C. Freysoldt, J. Neugebauer and C. G. Van de Walle, Phys. Rev. Lett., 2009, 102, 016402.
44. C. Freysoldt, J. Neugebauer and C. G. Van de Walle, Phys. Status Solidi B, 2011, 248, 1067–1076.
45. A. J. E. Rowberg, M. W. Swift and C. G. Van de Walle, Phys. Chem. Chem. Phys., 2021, 23, 14205–14211.
46. M. W. Chase Jr, J. Phys. Chem. Ref. Data, Monograph, 1998, 9, 1.
47. B. J. Van Zeghbroeck, Principles of Semiconductor Devices, University of Colorado, Boulder, 2011.
48. K. Hoang, C. Latouche and S. Jobic, J. Am. Ceram. Soc., 2022, 105, 4242–4249.
49. T. Yajima, H. Suzuki, T. Yogo and H. Iwahara, Solid State Ionics, 1992, 51, 101–107.
50. P. G. Sundell, M. E. Björketun and G. Wahnström, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104112.
51. M. E. Björketun, P. G. Sundell and G. Wahnström, Faraday Discuss., 2007, 134, 247–265.
52. M. W. Swift, A. Janotti and C. G. Van de Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 214114.
53. D. Han, K. Toyoura and T. Uda, ACS Appl. Energy. Mater., 2021, 4, 1666–1676.
54. J. Exner, T. Nazarenus, J. Kita and R. Moos, Int. J. Hydrog. Energy, 2020, 45, 10000–10016.
55. F. J. Loureiro, D. Pérez-Coll, V. C. Graça, S. M. Mikhalev, A. F. Ribeiro, A. Mendes and D. P. Fagg, J. Mater. Chem. A, 2019, 7, 18135–18142.
56. K. Toyoura, W. Meng, D. Han and T. Uda, J. Mater. Chem. A, 2018, 6, 22721–22730.
57. T. Fujii, K. Toyoura, T. Uda and S. Kasamatsu, Phys. Chem. Chem. Phys., 2021, 23, 5908–5918.
58. M. E. Björketun, P. G. Sundell and G. Wahnström, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 054307.
59. J. Hermet, M. Torrent, F. Bottin, G. Dezanneau and G. Geneste, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 104303.
60. R. Glöckner, M. Islam and T. Norby, Solid State Ionics, 1999, 122, 145–156.
61. K. Kreuer, E. Schönherr and J. Maier, Solid State Ionics, 1994, 70, 278–284.
62. S. J. Stokes and M. S. Islam, J. Mater. Chem., 2010, 20, 6258–6264.
63. H. Zhu, S. Ricote, C. Duan, R. P. O'Hayre, D. S. Tsvetkov and R. J. Kee, J. Electrochem. Soc., 2018, 165, F581–F588.
64. V. Gorelov, V. Balakireva and A. Kuz’min, Russian J. Inorg. Chem., 2018, 63, 930–937.
65. F. Zhao, Q. Liu, S. Wang, K. Brinkman and F. Chen, Int. J. Hydrog. Energy., 2010, 35, 4258–4263.
66. F. M. Draber, J. R. Denninger, P. C. Müller, I. K. Sommerfeld and M. Martin, Adv. Energy. Sustain. Res., 2022, 3, 2200007.
67. R. De Souza, J. Kilner and C. Jeynes, Solid State Ionics, 1997, 97, 409–419.
68. T. He, K. Kreuer, Y. M. Baikov and J. Maier, Solid State Ionics, 1997, 95, 301–308.
69. X. Xu, Y. Liu, J. Wang, D. Isheim, V. P. Dravid, C. Phatak and S. M. Haile, Nat. Mater., 2020, 19, 887–893.
70. S. M. Haile, D. L. West and J. Campbell, J. Mater. Res., 1998, 13, 1576–1595.
71. M. Shirpour, R. Merkle and J. Maier, Solid State Ionics, 2012, 216, 1–5.
72. T. Bondevik, J. M. Polfus and T. Norby, Solid State Ionics, 2020, 353, 115369.
73. M. Shirpour, R. Merkle and J. Maier, Solid State Ionics, 2012, 225, 304–307.
74. F. Iguchi, N. Sata and H. Yugami, J. Mater. Chem., 2010, 20, 6265–6270.
75. C. Kjølseth, H. Fjeld, Ø. Prytz, P. I. Dahl, C. Estournès, R. Haugsrud and T. Norby, Solid State Ionics, 2010, 181, 268–275.
76. M. Shirpour, B. Rahmati, W. Sigle, P. A. van Aken, R. Merkle and J. Maier, J. Phys. Chem. C, 2012, 116, 2453–2461.

---