Effects on the Proton Conduction Limiting Barriers and Trajectories in BaZr0.875Y0.125O3 Due to the Presence of Other Protons

Article information

J. Korean Ceram. Soc.. 2016;53(5):521-528
Publication date (electronic) : 2016 September 30
doi : https://doi.org/10.4191/kcers.2016.53.5.521
Department of Chemistry, Mount Holyoke College 50 College Street, South Hadley, MA 01075
Corresponding author :Maria A. Gomez, E-mail : magomez@mtholyoke.edu, Tel : +1-1-413-538-3868 Fax : +1-1-413-538-2327
Received 2016 April 18; Revised 2016 July 25; Accepted 2016 July 26.

Abstract

Kinetic Monte Carlo (KMC) and graph searches show that proton conduction limiting barriers and trajectories in BaZr0.875 Y0.125 O3 are affected by the presence of other protons. At 1000 K, KMC limiting conduction barriers increase from 0.39 eV to 0.45 eV as the proton number is increased. The proton-proton radial distribution begins to rise at 2 Å and peaks at 4 Å, which is half the distance expected, based on the proton concentration. Density functional theory (DFT) calculations find proton/proton distances of 2.60 and 2.16 Å in the lowest energy two-proton configurations. A simple average of the limiting barriers for 7–10 step periodic long range paths found via graph theory at 1100 K shows an increase in activation barrier from 0.32 eV to 0.37 eV when a proton is added. Both KMC and graph theory show that protons can affect each other’s pathways and raise the overall conduction barriers.

1. Introduction

Perovskite oxides have attracted much interest due to their possible uses in proton conducting fuel cells.126) The traditional mechanism for proton incorporation in perovskites is directly tied to the dopant concentration.15,27) When an ABO3 perovskite is doped at the tetravalent B sites with trivalent ions, an oxygen vacancy results for every two doped sites. This vacancy dissolves water from the atmosphere into a hydroxyl, which fills the vacancy and a proton which quickly binds to an oxygen in the perovskite lattice. Both protons freely conduct through the oxygen ion network of the perovskite. Since the mechanism suggests a 1 : 1 ratio of protons to trivalent dopant ions, previous proton conduction simulations in perovskites have considered only one proton and one dopant simulations.413,1526) However, one study28) on a different perovskite system finds low energy minima in a perovskite with protons only 2.3 Å apart due to lattice relaxation as well as pathways for moving the two protons together. The same study uses a simplified model to show that the two proton motion may in fact be the dominant diffusing defect.

This paper uses KMC to explore if a system with multiple protons initially placed far from each other can result in protons in close proximity. Additionally, the paper discusses the effect of a second proton on the probable paths of the first using graph theory and DFT on a small system. One might expect that an excess proton could change conduction by occupying the lowest energy binding sites, causing lattice distortions, and/or providing alternative conduction mechanisms through correlated proton motions. These could potentially decrease or increase conduction. To assess the impact of an excess proton on the conduction of a probe proton, this paper calculates the binding and transition state sites of a probe proton using the same electronic structure methods as used in an earlier paper for a single proton.25)

In principle, both the excess proton and the probe proton move through the lattice. Since a proton can bind to an oxygen in a perovskite with four possible hydroxyl orientations, 2326) after one proton is added to the system, there are only N–4 possible sites for the second proton where N is the total number of binding sites. As a result, there are N (N–4) / 2 combinations of two protonated sites when two protons are in the system. Since the total number of sites in our earlier single proton studies25) is 96, there are 4,416 possible combinations of two protons in this system. Finding the minimum energy for all those combinations and the transitions between them is not practical. As an approximation of the effect of an excess proton on a probe proton’s motion, we will simply fix the position of our excess proton at the best binding site for a single proton and consider the pathways of the probe proton based on energetics alone.

Conduction pathways for the probe proton are then found using graph theory. Each binding site is a vertex in a graph. Vertices are connected by edges when a single transition state exists between them. Probabilities for the pathways are assigned using a Boltzmann probability for the initial binding site and conditional probabilities based on approximate rate constants for moving through the remaining vertices in the path. Dynamic programming is used to search for all pathways through the vertices as in our earlier work.25,26)

Section II A describes the KMC29,30) approach used while section III A describes the resulting trajectories in systems with multiple protons with an equal number of dopant ions. Sections II B and III B describe how the binding sites, transition states, and rate constants are calculated and how they change with the presence of an excess proton, respectively. Sections II C and III C describe how long range proton pathways are calculated through graph theory and how the excess proton affects the long range conduction paths of the probe proton, respectively. These results suggest that proton motion is correlated and the first proton may prepare the way for the next.

2. Computational Procedure

2.1. Kinetic Monte Carlo including field generated by charged defects

The 12.5% Y doped BaZrO3 proton binding and transition state structures and energies found in our earlier study25) were used in a KMC study using the approach of Pornprasertsuk et al.30) In this approach, the effective barrier between binding sites is shifted by half the electric potential difference between the binding sites. The total electric potential at each site includes a contribution due to the fields created by charged defects as well as the applied field. Both of these contributions are calculated for planes perpendicular to the applied field. The new effective barrier is used to find the rate constants used in traditional KMC.29)

A 2 × 2 × 2 unit cell unit was replicated along the z-axis to form longer systems. Systems of one, two, four, and eight units were considered. This creates KMC simulation cells of length of 8.52, 17.04, 34.08, and 68.16 Å. The width and depth of the cells is 8.52 Å. Each original replicated unit contained one yttrium dopant ion. At the start a proton was placed in the lowest energy binding site in each replicated segment so that there are an equal number of protons and dopant ions.

The Vienna ab initio simulation package (VASP)31,32) was used in our earlier study25) to find the energies and normal mode frequencies at all the original unit binding and transition states. Specifically, the DFT method with the Perdew-Burke-Ernzerhof (PBE) functional was used to perform conjugate gradient optimization to find binding site energies and the Nudged Elastic Band (NEB)33) method was used to find minimum energy pathways between binding sites. The climbing NEB method33) was used to find the transition state. Finally, the dynamical matrix method34) was used to find normal mode frequencies at the binding and transition states. Details on these calculations are given in our earlier study.25) In this paper, we use the energies, barriers, and frequencies found for the original 2 × 2 × 2 unit cell system25) within the replicated longer systems.

The defect interactions between the z-replicated original system units are included in the electric potential contributions of the charged defects as described in the Pornprasertsuk KMC approach.30) The charged defect (proton and yttrium dopant) contribution includes the relative permittivity for barium zirconate which varies from 15 to 120 depending on experimental and calculation conditions.3540) References38,40) find the relative permittivity under conditions closest the system in this study (bulk and in Glazer (-,-,-) configuration) to be about 50, the value used in this simulation. The charge defect interactions damped by the relative permittivity along with the earlier DFT data were used to find the rate constants for proton motion between all binding sites. These rate constants were used in the KMC approach.

2.2. Binding and transition state calculations

As seen in our earlier study25) of 12.5% Y doped BaZrO3, the presence of a proton affects the structure of the lattice. Locally, oxygens get closer to the hydroxyl site making the interoctahedral distance closest to that site 3.4 Å rather than the more standard 3.8 Å in doped faces. This creates hydrogen bonds of about 2 Å. There are also some contractions of intraoctahedral OO distances near the excess proton. Further, as seen in the Fig. 1, interoctahedral transfer distances increase to 4.1 Å in the YZ plane making this plane’s distortions very similar to those in undoped planes. Earlier studies25) showed that proton conduction pathways without an excess proton have a significantly higher probability of occurring in the distorted doped planes. The increased interoctahedral distance in the YZ doped plane which effectively decreases the YZ plane’s distortions indicates that the presence of a proton could significantly impact the trajectory of another.

Fig. 1

The lattice distortions around the lowest energy single proton binding site are shown. This proton is referred to as the excess proton. While there are two interoctahedral distances in each cube face, only the shorter of these two distances is shown. The two planes containing the binding oxygen retain 3.8 Å or smaller distances between the closer interoctahedral oxygens while the remaining doped plane shows an increase in the shorter interoctahedral distance to 4.1 Å. Small contractions of intraoctahedral OO distances are also seen near the excess proton.

To find the effect of an excess proton at the lowest energy binding site on the energies of other binding sites, a probe proton was placed at each of the other possible binding sites and the energy was minimized. The same electronic structure methods as in our prior work25) which were briefly described in the previous section were used. However, with two protons and only one doped ion, a homogenous background charge of −1 is used to preserve charge neutrality. The NEB method33) was used to find transition states between neighboring binding sites using 2–6 intermediate images.

Using the binding and transition states found, rate constants are calculated using the standard harmonic transition state theory formula of kij(T)=Aijexp(-Eij-EjkBT).

Ai→j is the prefactor or the frequency of i → j attempts. Eij is the energy of the transition state between binding sites i and j. Ei is the energy of binding site i. As seen in our earlier studies,25) fixing the transition state energy pre-factor for all transitions at a constant value does not significantly affect the most probable pathways found and their average properties. For simplicity, we will fix the prefactor at 1 ps−1 for all transitions and effectively only consider the energetics of the pathways. The probability of moving from site i to site j given that the probe proton starts at site i (pi→j (T)) is the rate constant normalized by all the possible escape routes or kij(T)j=1Mkij(T). M is the number of binding sites connected to i via a single transition state. Notice that the normalization cancels out the constant prefactor making the specific constant choice irrelevant. These approximate rate constants are controlled fully by the transition energy barriers.

2.3. Finding long range proton pathways in the presence of an excess proton

Proton conduction pathways are found using a dynamic programming recursion on path weight to search the binding and transition state graph as in our earlier work.25,26) The binding sites are the vertices of the graph. These are connected by an edge only if there is a single transition state between them. In this and earlier papers,25,26) long range conduction pathways are defined as those that start at a binding site or vertex and end at a periodic image of that binding site traversing the entire simulation box. These pathways can be periodically replicated to achieve long range conduction. If the periodic pathway contains n steps, it is said to be n-step periodic. In another study, we find that this definition captures the essential features of long range conduction pathways.41) The overall probability of such an n-step periodic path (i1, i2,in, i1) is ρi1 (T)pi1i2(T)…pinin(T)pini1 (T) where ρi1(T) is the probability of being at site i1 at temperature T or the Boltzmann distribution at that site. The weight of the path, a more convenient quantity in the dynamic programming recursion, is minus the natural log of the probability of that path. The individual edge weights between sites i and j are wi → j (T) = − ln pi → j (T).

The n-step periodic path (i1, i2,in, i1) has n − 1 vertices between the starting and ending points and n distinct vertices. Dynamic programing is used to go through all subsets (S) of n distinct vertices and find the minimum weight of a path that starts at i1, visits the vertices in the subset once and ends at a vertex connected to the original vertex. Since all subsets are considered, all subsets with non-zero probability or non-infinite weights are kept and subsequently sorted.

To consider all subsets, we start paths at i1, define an array R with all vertices except for the first one, loop through all vertices j in R keeping only those that connect to the previous vertex. If vertices i and j are connected, pi → j (T ) is non-zero and wi→j is finite. All possible one step paths (i1, i2) are stored along with the current weight for that path. After only one step, the path weight would be ωi1i2(T). The next iteration considers all the connections to the second vertex in the currently stored paths. An array (R) with all vertices except for those already visited is defined. All possible next vertices from that array are considered and only those with edges connecting to the previous step are kept. Paths are duplicated when multiple next steps are found and terminated when there is no next step found. The latter is equivalent to finding only infinite weight edges. The path list is updated with all the new vertices and the path weights are increased by the new edge weight. The process continues until the paths have n vertices. In the final step, only the paths whose in vertex connects to the original vertex and that span the simulation box are kept. The path weights are augmented by − ln(ρi1 (T)) + wi1i (T)). The vertices removed from the array R for each path comprise the subset S for that path. All the paths found are sorted by probability or weight.

This algorithm is effectively what Scott et al.42) describes as dynamic programming using the overall path weight recurrence Wi1(S,T) = miniS−{i1} Wi(S−{i1}, T) + wi1i(T). The recurrence has been written in the notation of this paper. A final weight of −ln (ρi1 (T))+ wini1(T) is added to get the overall path probability in this physical setting. The final weight considers the probability of the initial vertex and the probability of the final connection from in back to the initial vertex. All possible initial vertices are considered. While the recurrence finds the minimum weight path, all possible n step paths are stored and sorted in order of increasing weight or decreasing probability. The stored paths are the ensemble of possible n-step periodic paths at temperature T. The ensemble contains the same paths in the forward and backward direction with equal probability.25)

Separate searches are done for each length path. To estimate the preferred periodic length, average estimated times for each length path can be found when rate constant prefactors are available.25,26) For this rough estimate of the effect of an excess proton on the conduction of a probe proton, estimated times are not calculated and a purely energetic assessment is made.

3. Results and Discussion

3.1. KMC trajectories show protons in very close proximity altering pathways

Kinetic Monte Carlo million picosecond trajectories at 1000 K were obtained for the original system replicated one, two, four and eight times in the absence of an external field. Initially, protons are placed in binding sites 8.52 Å apart along the longest side of the simulation box. The simulation box is periodic along the shorter two lengths. When a proton exits the replicated simulation box through the z-direction, a proton is added to the first layer on the z-axis based on the Boltzmann probability of having sites on that layer occupied. The process is repeated until the simulation time exceeds a million picoseconds. Limiting barriers were found by keeping a log of the barriers for protons traveling from entry to exit. The highest barrier seen in each proton trajectory was the limiting barrier for that proton trajectory. Over the million picoseconds of the simulation, many individual proton trajectories were found. To get the limiting conduction barrier, the individual limiting barriers were averaged weighted by the fraction of the time spent in their trajectory. The average limiting conduction barriers found at 1000 K for one, two, four, and eight proton systems were 0.39, 0.43, 0.44, and 0.45 eV, respectively. Intraoctahedral transfer limiting barriers represented 86%, 96%, 95% and 93% of the limiting barriers, for 1, 2, 4 and 8 proton systems, while rotational limiting barriers represent 12%, 2%, 0%, and 0% of limiting barriers, respectively. Finally, interoctahedral transfer limiting barriers represent 2%, 2%, 4%, and 7% of limiting barriers, respectively. This changing percentages of limiting barrier types for long range paths suggest that the presence of multiple protons alters chosen pathways.

When a field of 0.12 V/Å is added, the average limiting barriers change to 0.32, 0.35, 0.38, and 0.41 eV for 1, 2, 4, and 8 proton systems, respectively. This is a decrease relative to the zero field case as expected. More interestingly, the percent of each type of limiting barrier also changes. With the field, intraoctahedral transfer limiting barriers represent 56%, 69%, 81% and 89% of the limiting barriers, for 1, 2, 4 and 8 proton systems, while rotational limiting barriers represent 43%, 29%, 17%, and 9% of limiting barriers, respectively. Finally, interoctahedral transfer limiting barriers represent 1%, 1%, 2%, and 2% of limiting barriers, respectively. Overall, the field leads to more rotational limiting barriers and fewer intraoctahedral transfer limiting barriers. However, as protons are added, their motion along with the fixed lower valency dopant defects can create fields opposing the applied field. The competing effect of having more protons once again leads the intraoctahedral limiting barriers to become dominant. Once again, we see that protons affect the trajectories of other protons.

Let us consider the effect of the applied field in more detail. As in our earlier single proton work,45) the dopant defect acts as a trap for the proton. Occasionally, the trap is escaped through a long range limiting barrier which is most often an intraoctahedral transfer and less often a rotation. In the absence of field, the deepest potential energy wells are near the dopant though the trend is not monotonic. The field applied along the z-axis decreases the energy barrier to forward moves, increases the energy barrier to backward moves and leaves barriers to moves on planes perpendicular to the z-axis or field constant. The net effect is a tilting of the potential energy map lowering the region near the negatively charged electrode. This changes individual barriers differently. Using a flow of charge analogy and keeping in mind that the system has connections both in the field direction and in the planes perpendicular to it, the choke points in the charge transport can change as a result of the potential energy tilting. In fact, the proton residence time in some planes perpendicular to the field decreases while it increases in others due to multiple minima with different barriers and connections. The different planes have different connections to other planes leading to alternative routes and hence changes to the fraction of rotation and transfer limiting barriers. The 0.12 V/Å field tends to shift the residence time distribution and its choke points in the direction of the field moving the largest choke points from dopant planes to planes just after the dopant.

For all multiple numbers of protons, the proton-proton radial distribution shows a significant broad peak at about 4 Å as shown in Fig. 2. This broad peak shows that protons can get as close as 2 Å from each other. Since the energies used in this KMC simulation came from single proton DFT calculations with only the addition of proton/proton interactions through electric potential contributions of defects summed through planes of the simulation box, we next consider including the effects of additional lattice distortions by doing additional DFT calculations on the original small 2 × 2 × 2 unit cell system with the addition of a second proton.

Fig. 2

The proton proton radial distribution function shows a broad nearest neighbor peak at about 4 Å with some protons getting as close as about 2 Å to others.

3.2. Protons affect the binding sites, transition states, and rate constants of nearby protons

The binding sites found look very similar to that of the original system. Only binding sites on oxygens without the excess proton are considered. Not surprisingly, binding sites on the same oxygen as the excess proton led to movement of one of the protons to a nearby oxygen. Protons on sites near the excess proton form hydroxyls which are tilted away from the excess proton. The probe proton distorts the lattice in a similar way to the excess proton but also is affected by the initial distortions made by the excess proton. For example, the lowest energy binding sites are on the same face as the excess proton as shown in Fig. 3. Comparison with Fig. 1 shows that the probe proton allows for continued contraction of OO distances. The optimized two proton positions allow for significant hydrogen bonding for the probe proton and the excess proton. The proton/proton distances in the lowest and second lowest binding sites for the probe proton shown in Fig. 3 are 2.68 Å and 2.16 Å, respectively.

Fig. 3

The binding configuration of the lowest and second lowest energy binding sites in the presence of the excess proton is shown in (a) and (b). The excess proton is at the right of the diagram. Comparison with Fig. 1 reveals that the probe proton enhances OO distance contractions. The proton/proton distances in the in (a) and (b) are 2.68 Å and 2.16 Å, respectively.

Interoctahedral transitions through OO distances larger than 4.0 Å often proceeded through two intraoctahedral transfers rather than one interoctahedral transfer even with the smallest number of images. As a result, only interoctahedral transitions between oxygens closer than 4.0 Å were considered when finding paths. This resulted in considering only interoctahedral transfers in planes containing the oxygen binding the excess proton, namely the XZ and XY planes in Fig. 1.

In our earlier work,25) without an excess proton, transitions to doped octahedra were preferred over transitions away from doped octahedra. While the proton is still repelled less by the Y3+ dopant ion than by the Zr4+ ion, there are now additional factors to consider. The probe proton is also repelled by the positive charge of the excess proton. At the same time, the lattice distortions made by the excess proton stabilize some sites near the excess proton by allowing increased hydrogen bonding. Most forward and backward barriers follow the same patterns seen before the addition of the excess proton. However, some are changed by the additional competing factors of repulsion by the excess proton and stabilization of some sites near the excess proton. This results in some barriers for moving away from the dopant being lower in energy than the reverse barriers.

3.3. The excess proton affects the proton conduction paths of the probe proton

Figure 4 shows several probable probe proton conduction pathways at 1100 K given that the excess proton is fixed at the lowest energy single proton binding site. These are a small sampling of the many probable pathways. For periodic paths with 7–10 steps, the average limiting barrier through the paths is in the range of 0.31 to 0.41 eV. For comparison, the range for these same path lengths when a single proton is considered is 0.31 to 0.33 eV. A simple average between all path lengths considered suggests activation barriers of 0.32 eV and 0.37 eV in the absence and the presence of an excess proton, respectively. Like the KMC trajectories, dynamic path searches for a probe proton in the presence of an excess proton show an increase in average limiting barrier to conduction. The distribution of limiting step types, however, is less clear for this method as it varies widely by length of periodic path. The artificial selection of path length to odd numbers forces more interoctahedral limiting barriers as seen in our earlier studies.25,26) The odd numbered step periodic pathways often take a single interoctahedral transfer step (I) in preference to two intraoctahedral transfer steps (TT) and as a result have a sizable fraction of interoctahedral transfer path rate limiting steps. This shortcut leads to higher limiting barriers for odd numbered paths than even numbered pathways. In contrast, even numbered step periodic pathways have a very small fraction of interoctahedral rate limiting steps with and without the excess proton. Paths are limited by both rotations and intraoctahedral transfers. The excess proton increases the fraction of intraoctahedral transfer limiting paths for the 8-step periodic pathways while decreasing it for the 10-step periodic pathways. The fraction of rotation limited pathways is shifted in the opposite direction. Even numbered step paths in all directions are seen in the presence and absence of an excess proton. The excess proton slightly biases 8-step periodic paths in favor of the x-direction and 10-step periodic paths in favor of the z-direction. Overall, the excess proton impacts the probable probe proton pathways. Future investigations are underway considering different excess proton locations and concerted motions.

Fig. 4

Probable 7, 8, 9, and 10 step periodic pathways for proton conduction are shown in (a), (b), (c), and (d). The perovskite lattice of the minimum excess proton binding site is shown. The probe proton makes additional changes. However, since these differ from site to site, they are not displayed in these whole path images. The excess proton is shown in black.

4. Conclusions

KMC simulations show an increased limiting barrier to proton conduction with increasing number of protons at 1000 K. The barrier starts at 0.39 eV with a single proton and reaches 0.45 eV with eight protons. The number of protons also affects the distribution of limiting barriers among the three possible step types: rotations, interoctahedral transfer, intraoctahedral transfer showing that protons affect the trajectories of other protons. Further, the proton-proton radial distribution function found for these trajectories shows a broad nearest neighbor peak at 4 Å with some proton-proton distances as small as 2 Å despite initial proton distances of 8.52 Å. The presence of multiple protons also increases barriers in the presence of an external field though the overall barriers are smaller than without the external field. Further, the field does not impact the radial distribution function significantly.

Since the KMC simulation treats proton-proton interactions approximately and does not consider additional lattice distortions added by additional protons, binding sites, transition states, and conduction pathways for a probe proton in the presence of an excess proton have been found using DFT with the PBE functional. The lower energy binding sites for the probe proton are near the excess proton. This allows the probe proton to make effective use of the already slightly distorted lattice and suggests that proton motion is correlated when there are multiple protons present. The probable conduction pathways reveal that the first proton raises the average periodic path-limiting barrier. The rise is not as great as with the KMC calculation because this path calculation only looks at paths of a probe proton around a single excess proton with fixed location and hence does not sample all the possible pathways. A simple average between all path lengths considered at 1100 K suggests activation barriers of 0.32 eV and 0.37 eV in the absence and the presence of an excess proton, respectively.

The experimental activation barrier in the 900 – 1300 K temperature range is 0.44 eV17,43) which is in excellent agreement with our 1000 K KMC average limiting barrier range for 2–8 proton systems of 0.43 to 0.45 eV and higher than the single proton KMC limiting barrier of 0.39 eV. Prior calculations with a single proton have found overall conduction barriers of about 0.3 eV.5,25) and some individual interoctahedral barriers of about 0.4 eV.6,25) The graph pathway search yields lower limiting barriers the KMC though still increasing in the presence of an excess proton. Graph pathway searches for periodic pathways lead to more direct pathways than KMC trajectories and often have lower limiting barriers.44,45) The advantage of having done DFT calculations in the presence of two protons was the revelation that the probe proton benefits from the lattice distortion that the excess proton introduced.

Acknowledgments

We would like to thank Frederick G. Haibach for his insightful comments. This research was supported by an NSF RUI Grant No. CHE-1111474, a Henry Dreyfus Teach-Scholar award, and a Cottrell Science Scholars Multi-investigator award. Computational resources were provided in part by the MERCURY supercomputer consortium http://mars.hamilton.edu under NSF MRI No. CHE-122954.

References

1. Kreuer K. On the Development of Proton Conducting Materials for Technological Applications. Solid State Ionics 97(1):1–15. 1997;
2. Kreuer K. On the Complexity of Proton Conduction Phenomena. Solid State Ionics 136:149–60. 2000;
3. Kreuer KD. Proton-Conducting Oxides. Annu Rev Mater Res 33(1):333–59. 2003;
4. Björketun ME, Sundell PG, Wahnström G, Engberg D. A Kinetic Monte Carlo Study of Proton Diffusion in Disordered Perovskite Structured Lattices Based on First-Principles Calculations. Solid State Ionics 176(39):3035–40. 2005;
5. Björketun ME, Sundell PG, Wahnström G. Effect of Acceptor Dopants on the Proton Mobility in BaZrO3: A Density Functional Investigation. Phys Rev B 76(5):054307. 2007;
6. Merinov B, Goddard W III. Proton Diffusion Pathways and Rates in Y-doped BaZrO3 Solid Oxide Electrolyte from Quantum Mechanics. J Chem Phys 130(19):194707. 2009;
7. Bilić A, Gale JD. Ground State Structure of BaZrO3: A Comparative First-Principles Study. Phys Rev B 79(17):174107. 2009;
8. Bilić A, Gale JD. Proton Mobility in the In-doped CaZrO3 Perovskite Oxide. Chem Mater 19(11):2842–51. 2007;
9. Bilić A, Gale JD. Simulation of Proton Diffusion in In-doped CaZrO3. Solid State Ionics 179(21):871–74. 2008;
10. Islam MS, Davies RA, Gale JD. Proton Migration and Defect Interactions in the CaZrO3 Orthorhombic Perovskite: A Quantum Mechanical Study. Chem Mater 13(16):2049–55. 2001;
11. Liu Y, Yoshino M, Tatsumi K, Tanaka I, Morinaga M, Adachi H. Local Geometry and Energetics of Hydrogen in Orthorhombic SrZrO3. Mater Trans 46(6):1106–11. 2005;
12. Shi C, Yoshino M, Morinaga M. First-Principles Study of Protonic Conduction in In-doped AZrO3 (A=Ca, Sr, Ba). Solid State Ionics 176(11):1091–96. 2005;
13. Shi CS, Morinaga M. Doping Effects on Proton Incorporation and Conduction in SrZrO3. J Comput Chem 27(6):711–18. 2006;
14. Kreuer KD, Dipper T, Baikov Y, Maier J. Water Solubility, Proton and Oxygen Diffusion in Acceptor Doped BaCeO3: A Single Crystal Analysis. Solid State Ionics 86–88:613. 1996;
15. Kreuer KD, Munch W, Ise M, He T, Fuchs A, Traub U, Maier J. Defect Inter-Actions in Proton Conducting Perovskite-type Oxides. Ber Bunsen-Ges 101(9):1344–50. 1997;
16. Kreuer KD, Munch W, Traub U, Maier J. On Proton Transport in Perovskite-type Oxides and Plastic Hydroxides. Ber Bunsen-Ges 102(3):552–59. 1998;
17. Kreuer KD, Adams S, Munch W, Fuchs A, Klock U, Maier J. Proton Conducting Alkaline Earth Zirconates and Titanates for High Drain Electrochemical Applications. Solid Sate Ionics 145(1):295–306. 2001;
18. Munch W, Seifert G, Kreuer KD, Maier J. A Quantum Molecular Dynamics Study of Proton Conduction Phenomena in BaCeO3. Solid State Ionics 86–87:647. 1996;
19. Munch W, Seifert G, Kreuer KD, Maier J. A Quantum Molecular Dynamics Study of Cubic Phase of BaTiO3 and BaZrO3. Solid State Ionics 97(1):39–44. 1997;
20. Munch W, Kreuer KD, Seifert G, Maier J. A Quantum Molecular Dynamics Study of Proton Diffusion in SrTiO3 and CaTiO3. Solid State Ionics 125(1):39–45. 1999;
21. Munch W, Kreuer KD, Adams ST, Seifert G, Maier J. The Relation between Crystal Structure and the Formation and Mobility of Protonic Charge Carriers in Perovskite-type Oxides: A Case Study of Y-doped BaCeO3 and SrCeO3. Phase Transitions 68(3):567–86. 1999;
22. Munch W, Kreuer KD, Seifert G, Maier J. Proton Diffusion in Perovskites: Comparison between BaCeO3, BaZrO3, SrTiO3, and CaTiO3 Using Quantum Molecular Dynamics. Solid State Ionics 136:183–89. 2000;
23. Gomez MA, Griffin MA, Jindal S, Rule KD, Cooper VR. The Effect of Octahedral Tilting on Proton Binding Sites and in Pseudo-Cubic Perovskite Oxides. J Chem Phys 123(9):094703. 2005;
24. Gomez MA, Jindal S, Fletcher KM, Foster L, Addo ADA, Valentin D, Ghenoiu C, Hamilton A. Comparison of Proton Conduction in KTaO3 and SrZrO3. J Chem Phys 126(19):194701. 2007;
25. Gomez MA, Chunduru M, Chigweshe L, Foster L, Fensin SJ, Fletcher KM, Fernandez LE. The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO3, a Cubic Perovskite. J Chem Phys 132(21):214709. 2010;
26. Gomez MA, Chunduru M, Chigweshe L, Fletcher KM. The Effect of Dopant at the Zr Site on the Proton Conduction Pathways of SrZrO3: An Orthorhombic Perovskite. J Chem Phys 133(6):064701. 2010;
27. Jones CY, Wu J, LiPing L, Haile SM. Hydrogen Content in Doped and Undoped BaPrO3 and BaCeO3 by Cold Neutron Prompt-Gamma Activation Analysis. J Appl Phys 97(11):114908. 2005;
28. Bork N, Bonanos N, Rossmeist J, Vegge T. Thermodynamic and Kinetic Properties of Hydrogen Defect Pairs in SrTiO3 from Density Functional Theory. Phys Chem Chem Phys 13(33):15256–63. 2011;
29. Voter AF. Introduction to the Kinetic Monte Carlo Method. p. 1–23. Radiation Effects in Solids In : Sickafus KE, Kotomin EA, eds. Springer. Netherlands: 2007.
30. Pornprasertsuk R, Cheng J, Huang H, Prinz F. Electrochemical Impedance Analysis of Solid Oxide Fuel Cell Electrolyte Using Kinetic Monte Carlo Technique. Solid State Ionics 178(3):195–205. 2007;
31. Kresse G. Ab Initio Molekular Dynamik Für Flüssige Metalle. PhD thesis Technische Universittät; Wien: 1993.
32. Kresse G, Hafner J. Ab Initio Molecular-Dynamics for Liquid-Metals. Phys Rev B 47(1):558–561. 1993;
33. Henkelman G, Uberuaga BP, Jonsson H. A Climbing Image Nudge Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J Chem Phys 113(22):9901–4. 2000;
35. Azad A-M, Subramaniam S. Temperature Dependence of the Dielectric Response of BaZrO3 by Immittance Spectroscopy. Mater Res Bull 37(1):11–21. 2002;
36. Iguchi F, Tsurui T, Sata N, Nago Y, Yugami H. The Relationship between Chemical Composition Distributions and Specific Grain Boundary Conductivity in Y-doped BaZrO3 Proton Conductors. Solid State Ionics 180(6):563–68. 2009;
37. Yamazaki Y, Hernandez-Sanches R, Haile SM. High Total Proton Conductivity in Large-Grained Yttrium-doped Barium Zirconate. Chem Mater 21(13):2755–62. 2009;
38. Babilo P, Uda T, Haile S. Processing of Yttrium-doped Barium Zirconate for High Proton Conductivity. J Mater Res 22(05):1322–30. 2007;
39. Kjolseth C, Fjeld H, Prytz O, Dahl PI, Estournes C, Haugsrud R, Norby T. Space-Charge Theory Applied to Grain Boundary Impedance of Proton Conducting BaZr0.9Y0.1O3−δ. Solid State Ionics 181(5):268–75. 2010;
40. Bennett JW, Grinberg I, Rappe AM. Effect of Symmetry Lowering on the Dielectric Response of BaZrO3. Phys Rev B 73(18):180102. 2006;
41. Gomez MA, Shepardson D, Nguyen LT, Kehinde T. Periodic Long Range Proton Conduction Pathways in Pseudo-Cubic and Orthorhombic Perovskites. Solid State Ionics 213:8–13. 2011;
42. Scott J, Ideker T, Karp RM, Sharan R. Efficient Algorithms for Detecting Signaling Pathways in Protein Interaction Networks. p. 1–13. Research in Computational Molecular Biology 3500Springer Berlin Heidelberg. Heidelberg: 2005.
43. Bohn HG, Schober T. Electrical Conductivity of the High-Temperature Proton Conductor BaZr0.9Y0.1O2.95. J Am Ceram Soc 83(4):768–72. 2000;
44. Gomez MA, Liu F-J. Protons in Al doped BaZrO3 Escape Dopant Traps to Access Long Range Proton Conduction Highways. Solid State Ionics 252:40–7. 2013;
45. Krueger R, Haibach F, Fry D, Gomez M. Centrality Measures Highlight Proton Traps and Access Points to Proton Highways in Kinetic Monte Carlo Trajectories. J Chem Phys 142(15):154110. 2015;

Article information Continued

Fig. 1

The lattice distortions around the lowest energy single proton binding site are shown. This proton is referred to as the excess proton. While there are two interoctahedral distances in each cube face, only the shorter of these two distances is shown. The two planes containing the binding oxygen retain 3.8 Å or smaller distances between the closer interoctahedral oxygens while the remaining doped plane shows an increase in the shorter interoctahedral distance to 4.1 Å. Small contractions of intraoctahedral OO distances are also seen near the excess proton.

Fig. 2

The proton proton radial distribution function shows a broad nearest neighbor peak at about 4 Å with some protons getting as close as about 2 Å to others.

Fig. 3

The binding configuration of the lowest and second lowest energy binding sites in the presence of the excess proton is shown in (a) and (b). The excess proton is at the right of the diagram. Comparison with Fig. 1 reveals that the probe proton enhances OO distance contractions. The proton/proton distances in the in (a) and (b) are 2.68 Å and 2.16 Å, respectively.

Fig. 4

Probable 7, 8, 9, and 10 step periodic pathways for proton conduction are shown in (a), (b), (c), and (d). The perovskite lattice of the minimum excess proton binding site is shown. The probe proton makes additional changes. However, since these differ from site to site, they are not displayed in these whole path images. The excess proton is shown in black.