- Email: [email protected]
Ordered yttrium concentration effects on barium zirconate structure, proton binding sites and transition states
Solid State Ionics 304 (2016) 126–134
Contents lists available at ScienceDirect
Solid State Ionics journal homepage: www.elsevier.com/locate/ssi
Ordered yttrium concentration effects on barium zirconate structure, proton binding sites and transition states Maria A. Gomez a, * , Gillian Kwan a , Wanshu Zhu a , Monica Chelliah a , Xintong Zuo b , Audrey Eshun a , Virginia Blackmer a , Truc Huynh a , Mai Huynh a a b
Department of Chemistry, Mount Holyoke College, South Hadley, MA 01075,United States Medical School, Duke-NUS, 169857, Singapore
A R T I C L E
I N F O
Article history: Received 30 October 2016 Received in revised form 22 March 2017 Accepted 27 March 2017 Available online xxxx
A B S T R A C T Acceptor doped perovskites show promise as fuel cell proton conducting ceramic membranes and continue to be rapidly developed. Proton conduction in barium zirconate with 6.25% ordered yttrium doping at the zirconium site is explored. The dopant increases octahedral distortions in the lattice similar to our earlier studied 12.5% dopant system though the angle distribution is slightly broader. Normal modes promoting vibrations about a non-transient octahedral tilt are found for both concentrations at frequencies thermally accessible at typical proton conduction temperatures. The energy of proton binding sites and transition states roughly increases with distance away from the dopant though there are many oscillations in the trend due to local proton induced distortions. Protons can both locally enhance and decrease the octahedral distortions induced by the dopant. Transition state energies roughly increase with proton-dopant distance. Rotational and intraoctahedral transfer transition state energies vary through a similar magnitude range though the intraoctahedral transition states are often higher in energy. Interoctahedral transfer transition states tend to have the highest energies. The diversity of minima and transition states suggests that proton conduction in 6.25% yttrium doped barium zirconate is more complex than conduction in 12.5% doped barium zirconate. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Acceptor doped perovskite oxides, especially yttrium doped barium zirconate, show increasing promise as fuel cell high to intermediate temperature proton conductors [1–3]. Central to improving performance of these materials is understanding how the energy landscape for proton conduction changes with dopant concentration. High proton concentration and fast conduction paths for individual protons lead to high conductivity. When barium zirconate is doped with yttrium, Y3+ occupies Zr4+ sites through a maximum load of 16% [3,4], forming charge-compensating oxygen vacancies. In humid conditions, the water can fill oxygen vacancies. Electron density shifts allow an oxide ion to remain in the former vacancy while two protons conduct through the perovskite moving from one electron rich oxide ion to another. Thus, the greater the yttrium content, the greater the expected proton concentration. However, since the yttrium dopant distorts barium zirconate effectively changing the
* Corresponding author. E-mail address: [email protected] (M.A. Gomez).
http://dx.doi.org/10.1016/j.ssi.2017.03.027 0167-2738/© 2016 Elsevier B.V. All rights reserved.
landscape in which protons conduct, increasing yttrium content may not lead to consistent increases in proton conduction. Many experiments have elucidated specific details of the proton conduction landscape. Braun et al. used X-ray and neutron diffraction studies to show that yttrium is arranged in an ordered superstructure in the barium zirconate lattice and noted that proton conduction occurs on a landscape set by the dopant [5]. Giannici et al. [4] proposed a landscape including octahedral tilting in light of the appearance of additional Raman peaks in the yttrium doped barium zirconate system and noted that these led to a permanent lowering from the cubic symmetry seen in previous studies [6]. Yamazaki et al. showed that to conduct, protons need to overcome a barrier to escape from the dopant and also a barrier to conduct once free from the dopant trap giving an overall migration barrier of 0.47 eV below 470 K [7]. Several other experiments have found migration barriers in this range for temperatures in the range 350–900 K [3,8,9]. Both Braun et al. and Yamazaki et al. have found that migration barriers decrease significantly at high temperatures suggesting a change in mechanism which lowers the activation energy significantly though they differ on the specific temperature critical to the change [7,10]. Characterizing structure flexibility and energetics may shed light on potential mechanism change at higher temperatures.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
Simulations see similar proton conduction landscape patterns and can shed further light on the specific proton conduction mechanisms in barium zirconate. Octahedral tilting has been seen to increase upon yttrium doping using Density Functional Theory (DFT) [11–13]. Very different proton conduction pathways have been found for different acceptor dopants that varied the extent of octahedral tilting [14,15]. Björketun et al. have found that larger dopants result in a larger trapping region in the conduction landscape decreasing energies of proton binding sites with high hydrogen bonding in second coordination shells. Using kinetic modelling, they also found that larger dopants had decreased activation barriers to conduction. A migration barrier of 0.30 eV was found for 12.5% yttrium doped barium zirconate in the 300–600 K temperature range [16,17]. In 12.5% Y doped BaZrO3 , Merinov and Goddard [18] used DFT to find activation energies of 0.48 to 0.49 eV and 0.41 eV for a variety of intraoctahedral and interoctahedral transfers, respectively. van Duin et al. [19] found an activation energy of 0.45 eV with a yttrium dopant concentration of 12.5% using molecular dynamics with a ReaxFF reactive force field. Using a combination of DFT, kinetic Monte Carlo, and graph theory, we have found that long range periodic pathways in the same system have limiting barriers of 0.3 eV [13] while kinetic Monte Carlo pathways which allow for the indirect meandering of diffusion have barriers of 0.39 eV [14] for a single proton system and 0.45 eV in an eight proton system at 1000 K [20]. These findings are in general agreement with experimental proton landscape details and low temperature migration barriers. Several experiments have highlighted changes in the conduction landscape with dopant concentration. A study by Kreuer et al. shows that proton mobility increases with increasing yttrium dopant concentration through a range of 2% to 20% though the mobility increase saturates by 10% doping. At the same time, the activation barrier remains 0.43–0.44 eV for 2–10 % and goes up to 0.48 eV by 20% in the temperature range of 313–625 K [3]. Hydration shows a similar pattern to mobility but starts to decrease from 20% to 25% dopant concentration [3]. Another study by Fabbri et al. [9] probed the yttrium concentration range of 20–60% in barium zirconate which is beyond the 16% concentration where yttrium occupies exclusively zirconium sites [4]. Proton concentration increased with dopant concentration and decreased with temperature [9] as expected and as seen in the earlier Kreuer et al. study [3]. Migration barriers [9] at the lower concentration range were also consistent with the earlier study [3]. This paper considers the effect of decreased yttrium concentration from 12.5% [13] to 6.25% in cases where there is a single acceptor dopant ordered through periodic boundary conditions. There are other groups [21] currently enumerating many possible multiple dopant and oxygen vacancy configurations. Such studies are critical as at least one experiment suggests that yttrium is arranged in an ordered super structure in barium zirconate [5]. In our prior work with 12.5% Y doped BaZrO3 [13], one in eight Zr sites contained yttrium in a periodic simulation box. As a result, every other octahedron in each direction contained a yttrium ion, a highly ordered acceptor dopant situation. In this study, the system is twice as large in one direction and again one Zr site is occupied by a yttrium ion which is then periodically replicated. Increasing size in all directions would lead to needing to consider multiple relative yttrium configurations and an overall significantly higher computational cost. Section 2 describes all the methods used in this contribution. Section 3.1 describes how lowering the dopant concentration from 12.5% to 6.25% affects the overall perovskite structure, tilting angle distribution, and low energy normal modes. Sections 3.2 and 3.3 describe the coupled effect of yttrium induced octahedral tilting and proton induced local distortions of the structure and energy of the binding sites and transition states, respectively. Finally, Sections 4 and 5 contain the discussion and conclusion.
127
2. Methods Bilic´ and Gale [12] show that even the structure of the undoped barium zirconate global energy minimum has alternating layers with opposite octahedral tilting in each direction corresponding to (−, −, −) in Glazer notation [22]. Bennet, Grinberg and Rappe also found octahedral tilting in the undoped barium zirconate structure [11]. In an earlier study [13], we found that the structure of 12.5% yttrium doped BaZrO3 showed greater octahedral distortions than the undoped system. In this study, we lower the yttrium concentration to 6.25%. To obtain the initial structure, we double the size of our earlier studied 2 × 2 × 2 unit cell system with 12.5% yttrium at the zirconium site and replaced one dopant ion with a zirconium ion thereby setting the doping level to 6.25%. The initial structure was in the lowest energy Glazer distortions [22] namely (−, −, −). Density functional theory (DFT) calculations with the PBE functional and a generalized gradient approximation (GGA) have been performed to find energies using the ab-initio total-energy and molecular-dynamics program VASP (Vienna ab-initio simulation program) developed at the Fakultät für Physik of the Universität Wien [23–28]. The projector augment wave (PAW) method [29] with the valence states of Ba(5s, 5p, 6s, 5d), Zr(4s, 4p, 5s, 4d), Y(4s, 4p, 5s, 4d), and O(2s, 2p) [30] is used. Partial occupancies are calculated with Gaussian smearing starting at 0.2 eV and automatically extrapolated to zero using VASP. Bands are optimized iteratively using the preconditioned residual minimization scheme. The first 10 iterations are non-self consistent. A 2×2×2 Monkhorst-Pack k-point mesh with zero shift is used. Accurate level precision and real space projections are done. The initial dimensions of the system were set to 2a×2a×4a where the a is the lattice constant of the undoped barium zirconate system of 4.26 Å found in our earlier study with the PBE functional [13]. This lattice constant is larger than the experimental lattice size of 4.19 Å [4] but similar to other DFT GGA VASP studies as seen in Table 1 of reference [12]. The LDA functional does match the experimental value better and gives similar structure patterns to the GGA functions [31] though GGA functionals are expected to do better with transition states. The conjugate gradient minimization method was used to find the optimum configuration of the nuclei in the 2a×2a×4a system with one yttrium dopant ion at a zirconium site. A uniform positively charged cloud is used to neutralize the charge. This uniform cloud naturally misses the non-uniform effect of localized charge compensating oxygen vacancies but is a common approximation in periodic systems when one does not want to enumerate all the possible configurations of two to one yttrium ion to oxygen vacancy ratios at each yttrium concentration. This structure was further allowed to relax through changing of both ionic positions and cell shape and volume. While the symmetry group was not explicitly set, VASP by default identifies the symmetry of the original structure and maintains it through the optimization. In this case, the minimal C1 symmetry is identified likely because earlier work [13] used random angle octahedral tilting in different directions to explore the different possible Glazer configurations [22]. Glazer does identify the (−, −, −) tilt system as having F 1¯ symmetry but we did not explicitly fit to or constrain for those symmetry elements. Further, as seen in Section 3.1, the length and angles of the simulation box suggest nearly cubic average unit cells. Once the lowest energy structure is found, normal modes are calculated. Our earlier 0% and 12.5% yttrium concentration system had fixed the lattice parameter at 4.26 Å [13] and box angles at 90◦ while the eight octahedra were allowed to tilt. To see the effect of dopant concentration on lattice size, the simulation box for these two concentrations is optimized with the same conditions as our earlier work [13] except that we allow lattice shape and size to relax in addition. Note, however, that all other comparisons with earlier data are comparing the earlier work with the 4.26 Å lattice constant [13].
128
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
Proton binding sites at 6.25% yttrium concentration are found using the conjugate gradient minimization method to relax ion positions but not cell shape or volume. Four initial proton locations are considered about each oxygen with the OH bond pointing into each cubic plane as in our earlier work [13]. Optimization is halted when forces are less than 0.01 eV/Å. Energies are expected to be accurate to about 0.02 eV. Relaxed transition state energies are found using a climbing nudged elastic band method [32]. Initially, a standard nudged elastic band method is used followed by a climbing nudged elastic band method where the highest image is set to climb. In an effort to find as many transition states as possible, the looser convergence criteria of converging forces at the transition state to 0.1 eV/Å was set. For some transition states, convergence in force was set to 0.01 eV/Å. The energy changes to the barriers are expected to have errors around 0.04 eV. 3. Results 3.1. Octahedral tilting in the optimized doped perovskite without protonation The final optimized simulation box dimensions for the 6.25% yttrium concentration case were 8.46 Å by 8.47 Å by 16.97 Å with angles 90.0◦ , 90.0◦ , and 90.1◦ for this 2×2×4 unit cell system. Table 1 shows the average lattice constants and angles for 0%, 6.25%, and 12.5% yttrium simulations where the box size and shape is allowed to change. The experimental lattice constant column shows individual lattice constants for cubic systems and two lattice constants for tetragonal systems. The calculated lattice constant values overestimate experimental values [3,4] though they do show the same increasing trend with dopant concentration. Notice that the 12.5% yttrium system shows a larger lattice constant when the unit cell is allowed to vary in shape and size than in our earlier work which held it at 4.26 Å [13]. The remainder of the 12.5% comparisons in this paper use the data from reference [13] which kept a box size that is closer to the experimental value at the higher concentration. Including oxygen vacancies and starting with different Glazer octahedral tilting can impact calculated lattice constants. Visual inspection of the structure shows a Glazer octahedral tilting pattern of (−, −, −). Fig. 1 compares one side of the optimized 6.25% doped structure with the earlier 12.5% doped structure [13] suggesting that octahedral tilting roughly decreases with decreasing dopant concentration. Octahedral tilting brings two oxide ions closer together while bringing two other oxide ions further apart. This results in a variation of O-O distances across cubes as shown in Fig. 1. The greater range in O-O distances in the higher doped system is a result of the larger octahedral tilting. For a more quantitative comparison, the distribution of ZrOZr and ZrOY angles is explored. In a perfectly cubic system without octahedral tilting these angles would be 180◦ . Fig. 2 shows the distribution of ZrOX angles with X=Zr and X=Y for both 12.5% and 6.25% yttrium doped systems. All angles are less than 180◦ highlighting that both Table 1 The average unit cell lengths are shown for optimized structures allowing the cell shape and size to relax. Experimental lattice constants are shown in the last column for comparison. The experiments with a single lattice constant found cubic symmetry while those with two lattice constants found tetragonal symmetry. %Y
Average lattice constants (Å)
Angles (◦ )
Experiment lattice constants (Å)
0 5 6.25 12.5 15
4.23, 4.22, 4.22
90.0, 90.1, 90.0
4.19 [4] 4.21 [3]
4.23, 4.24, 4.24 4.29, 4.29, 4.29
90.0, 90.0, 90.1 90.1, 90.1, 90.1 a=4.23, c=4.21 [3]
systems exhibit some octahedral tilting. The average ZrOZr angles are 170◦ and 172◦ with range widths of 6◦ and 10◦ for 12.5% and 6.25% yttrium, respectively. The average ZrOY angles are 164◦ and 165◦ with range widths of 1◦ and 8◦ for 12.5% and 6.25% yttrium, respectively. Both averages show a decrease in octahedral tilting or an increase in angle as the yttrium concentration decreases. The details of the distributions reveal more significant subtleties than the averages. The ZrOZr angle distribution in the lower yttrium concentration system has more larger angles i.e. angles showing less tilt or distortion from the cubic case. However, there is a small number of angles that show increased tilting. The lower yttrium concentration ZrOY angle distribution shows two smaller angles and two larger angles with an overall higher average angle showing overall decreased tilting. It is important to note that this 2 × 2 × 4 system is only lower in concentration in one direction. The contrasting higher concentration in the perpendicular plane may be leading to the increased tilting peaks. An earlier experimental study [4] noted that EXAFS data could be described by a ZrOY angle of about 160◦ . Both concentrations considered here have most ZrOY angles in the range of 160–165◦ . O-Zr distances are 2.13(3) Å and 2.12(2) Å at 12.5% and 6.25% doping, respectively. O-Y distances are 2.23(1) Å and 2.22(1) Å at 12.5% and 6.25% doping, respectively. The value in parenthesis is one standard deviation of the distribution at the last significant figure. XRD and EXAFS data by Giannici [4] in Fig. 10 shows distances within a standard deviation of these values. It is worthwhile noting that the tilt system used to fit the data in Giannici [4] is (0, 0, +) while our earlier work [13] predicted a (−, −, −) tilt system for 12.5% doping and had the (0, 0, +) case as 0.17 eV higher in energy. The (−, −, −) tilt system was also found in the undoped case in another study [12]. The zirconium and yttrium ions seem to stay very close to the center of the octahedra in the minimum energy structures. Normal mode analysis of the minimum energy structures shows many low energy vibrations involving zirconium and barium oscillations and higher energy octahedral tilting vibrations. One imaginary mode with a frequency corresponding to an energy less that 0.01 eV was found for the low and high dopant concentrations. Since the energy of this mode is less than the error expected in these calculations and moving in the mode direction did not yield a lower energy, the structure was considered to be a true minimum. The lowest energy normal mode at 6.25% and 12.5% yttrium concentrations has frequencies of 38 cm −1 and 40 cm −1 , respectively. These frequencies are again lower than the accuracy of the energy calculations and correspond to 55 K and 58 K. Low frequency modes involve vibrations of the barium ions and some zirconium ions within their octahedra. Modes at 200 cm −1 and 199 cm −1 for 6.25% and 12.5% yttrium concentrations, respectively, involve ZrOX angle changes and are indicative of octahedral tilting. These higher frequencies correspond to temperatures of 288 K and 286 K. Other higher frequency normal modes also involve tilting. This suggests that at typical proton conduction temperatures that low frequency vibrational modes will involve moving of cations in their sites and higher frequency though still thermally accessible modes will show octahedral tilting oscillations. Higher frequency modes (up to 787 cm −1 and 763 cm −1 in low and high dopant systems) show even more octahedral tilting vibrations but are less thermally accessible. Overall, even at moderate temperatures, there is enough thermal energy that the structure is flexible and small normal mode oscillations around the lowest energy structure are possible. 3.2. Structure and trends for proton binding sites As shown in Fig. 3, the relative energy for proton binding sites (black squares) roughly increases as a function of distance away from the dopant up to about 5.5 Å. The energy shown is relative to the lowest energy proton binding site. The trend is not monotonic due to proton induced local distortions which can change the proximity of
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
129
Fig. 1. (a) shows a section of the 6.25% yttrium doped barium zirconate system while (b) shows the 12.5% doped system. Oxide ions are shown in red. Barium ions are shown in pink. Octahedra around zirconium ions are in green and octahedra around yttrium ions are in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the proton and non-bonded oxide ions beyond the dopant induced tilting. We will refer to the decreasing of the distance between the proton and non-bonded oxide ions as increasing hydrogen bonding. Fig. 4 shows the corresponding data from our earlier work [13] on 12.5% doping to be clustered. The trend is for increasing relative energy as a function of proton distance from the dopant. The one exception to the increase in energy occurs at slightly under 4.0 Å where local distortions bring some binding sites on oxide ions that are second nearest neighbors to the dopant to a lower energy. The low energy for protons on some first and some second nearest neighbor oxide ions are more spread out in the 6.25% yttrium case than the 12.5% case. The increased octahedral tilting in the 12.5% yttrium system results in a reduction of binding site diversity. Fig. 5 shows three binding sites around the same oxide ion for the 6.25% yttrium case. In each case, the proton induces significant local distortions. In binding site (a), the proton brings two oxide ions which were 4.35 Å apart in the unprotonated system to a distance of 3.57 Å apart. This allows both oxide ions to be closer to the proton.
Since all octahedra in the system are connected by vertices and the long side of this simulation box is just four octahedra long, this proton induced distortion reverses the dopant induced distortion in this configuration. This is an artificial effect of the boundary conditions. Releasing the periodic constraints by freezing the plane of octahedra in the center of this simulation box would have prohibited long range tilting changes. However, the artificial system result highlights that when the dopant concentration is low, local proton induced tilting can be significant. In binding site (b), the proton increases the original dopant induced tilting bringing two oxide ions from 4.07 Å apart to 3.51 Å apart. Again, this brings both oxide ions closer to the proton. Binding site (c) shows that the proton can also bring the oxide ion in the direction of the OH bond closer. The distance between the oxygen in the hydroxide and the opposite oxide ion started at 4.05 Å and contracted to 3.55 Å. In each of these cases, the proton brings oxide ions closer enhancing hydrogen bonding. Overall, protons can induce oxide ions to approach each other or move apart through local tilting of octahedra around zirconium
12 12.5% X = Zr 12.5 % X = Y 6.25% X = Zr 6.25% X = Y
Number of Observations
10
8
6
4
2
0 160
165
170
175
180
Fig. 2. The distribution of ZrOX angles is shown for both 12.5% and 6.25% yttrium concentrations. Traditional histogram rectangles are used for the higher concentration while points are used for the lower concentration. ZrOZr angle distributions are shown in black while ZrOY angle distributions are shown in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
130
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
0.90 Inter TS Intra TS Rot TS Minima
0.80
0.70
Relative Energy (eV)
0.60
0.50
0.40
0.30
0.20
0.10
0.00 2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Distance from dopant (Å) Fig. 3. The relative energy of minima and transition states is shown as a function of proton-dopant distance for 6.25% dopant concentration. All energies are relative to the global minimum for this concentration. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
and yttrium ions. The reversal in the tilting induced by the dopant in some cases is likely just an artifact of small system size and periodic boundary conditions. Figs. 6 and 7 highlight that the ZrOX distribution with X=Zr and Y can widen significantly upon addition of a proton in both the 6.25% and 12.5% yttrium doped systems. The three binding sites considered in the higher dopant concentration case have a hydroxide in a plane without dopant. While both distributions broaden relative to the no proton cases, the lower concentration ZrOX distribution has lower octahedral tilting or angles closer to 180◦ for the binding sites far from the dopant.
3.3. Transition states connecting proton binding sites Fig. 3 shows the energy of all the transition states between proton binding sites relative to the binding site global minimum. In general, all transition state energies roughly increase as the proton-dopant distance increases until about 5.5 Å. Transition states energies for rotations (green plus symbols) and intraoctahedral transfers (red triangles) are similar with the latter more often being higher in energy. The highest transition state energies are for interoctahedral transfers (blue circles). As we saw in Section 3.2, the backbone changes due to different local distortions in the different minima.
0.90 Inter TS Intra TS Rot TS Minima
0.80
0.70
Relative Energy (eV)
0.60
0.50
0.40
0.30
0.20
0.10
0.00 2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
Distance from dopant (Å) Fig. 4. The relative energy of minima and transition states is shown as a function of proton-dopant distance for 12.5% dopant concentration. All energies are relative to the global minimum for this concentration.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
131
Fig. 5. Three different binding site structures on the same oxide ion are shown in (a), (b), and (c).
Not all possible backbone distortions were considered to find the transition states so our data represents only a subset of possible minimum energy configurations and transition states. Fig. 4 shows the corresponding transition states for the 12.5% dopant concentration system. The overall pattern is the same but the energies for binding and transition states are more localized in the higher concentration system. A variety of rotational, intraoctahedral transfer, and interoctahedral transfer transitions are seen in the 6.25% yttrium concentration system. A small sampling is presented in Figs. 8, 9, and 10. Fig. 8 highlights a rotation including a backbone distortion. The distance between two oxide ions is shown to highlight the face in which the backbone change is most prominently seen. Notice that the two oxide ions tracked go from being the closer two oxide ions on that face to being the further two oxide ions on that cube face. At the same time, the dopant-proton distance increases from 2.32 Å to 3.06 Å to 3.38 Å while the proton-barium distance goes from 2.54 Å to 2.42 Å to 2.50 Å. Rotations tend to bring the proton in close proximity to barium at the transition state. To try to overcome this repulsion, the system distorts allowing both hydroxide and barium to move away as in other studies [33]. These two oxide ions near the proton in (a), (b), and (c) with average distances to the proton of 2.17 Å, 2.58 Å, and 2.15 Å, respectively. The array of competing factors makes it difficult to predict whether the forward or backward barrier is higher for this transition. Both forward direction increases in dopant-proton distance as well as sharper decreases in barium-proton distance
suggest a higher forward barrier. The average hydrogen bond length, however, has a larger increase in the backward direction suggesting a higher backward barrier. No doubt, there are other factors not yet considered as well. The actual forward barrier of 0.21 eV is slightly lower than the reverse barrier of 0.26 eV. Fig. 9 shows an intraoctahedral transfer that retains the same basic backbone with smaller proton induced changes. This transfer occurs on a face with higher dopant concentration and has a larger barrier to leave the dopant than to reach it namely 0.32 and 0.25 eV, respectively. Fig. 10 shows an interoctahedral transfer in the higher concentration plane and near the dopant. Forward and backward barriers are 0.29 and 0.23 eV, respectively, with the barrier for leaving the dopant (forward) being slightly larger. The broad change in relative energy of both minima and transition states as the proton/dopant distance increases seen in Fig. 3 for the 6.25% yttrium concentration case leads to a very broad pattern for barriers to moving closer and further from the dopant. Individual cases highlighted show that while sometimes the energy decreases as the dopant is approached, there are local proton induced changes that may counter this effect. This can be a result of increased hydrogen bonding or decreased proton-barium repulsion. At 12.5% yttrium concentration, Fig. 4 shows significantly more clustering of binding sites and transition state energies leading to clearer patterns in activation barriers as the proton-dopant distance increases. The patterns are still prone to oscillations due to local distortions caused by hydrogen bond enhancement and proton-barium repulsion.
12 no protons X = Zr no protons X = Y (a) X = Zr (b) X = Zr (c) X = Zr (a) X=Y (b) X = Y (c) X = Y
Number of Observations
10
8
6
4
2
0 155
160
165
170
175
180
Fig. 6. The distribution of ZrOX with X=Zr and Y for Fig. 5 (a), (b), and (c) proton binding site structures are shown.
132
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134 12 no protons X = Zr no protons X = Y (a) X = Zr (b) X = Zr (c) X = Zr (a) X = Y (b) X = Y (c) X = Y
Number of Observations
10
8
6
4
2
0 155
160
165
170
175
180
Fig. 7. The distribution of ZrOX with X=Zr and Y for three binding sites (a), (b) and (c) on an oxide ion on a plane between the dopant planes in the 12.5% yttrium doped system.
4. Discussion The optimized 6.25% yttrium doped barium zirconate structure shows slight octahedral tilting while still retaining roughly cubic average unit cells. Visual comparison with the 12.5% yttrium system shows that decreasing the yttrium concentration decreases octahedral tilting. Comparison of distributions of ZrOX angle with X=Zr and Y between 6.25% and 12.5% yttrium doped barium zirconate shows roughly increased tilting with increasing yttrium concentration. Low energy normal modes with zirconium and barium ions oscillating in their lattice positions are seen in addition to higher energy normal modes with some octahedral tilting. Some of these octahedral tilting normal modes have thermally accessible energies at proton conduction temperatures suggesting that small octahedral tilting oscillations maintaining overall tilt direction are possible during proton conduction. Octahedral tilting vibrations at higher frequencies have been calculated before for undoped barium zirconate by the Rappe group [11]. Since there is more octahedral tilting in the doped systems, a lowering in the frequency of these vibrations is not unexpected. Transient octahedral distortions were suggested for the undoped system in addition to non-transient distortions for doped systems by Giannici et al. [4]. ZrOY angles, and O-Zr and O-Y distances are in also in good agreement with experiment [4].
The energy of proton minima roughly increases as the proton is moved away from the dopant. The increased diversity of backbone shifts leads to an increase in the spread of minimum energies in the 6.25% dopant case in comparison to the more clustered 12.5% dopant energies. Both show low energies for protons bound to oxide ions in both the first and second neighbor shells of the yttrium as seen in other studies [16,34,35]. Overall, energy trends are not only affected by the dopant induced octahedral tilting but also by proton induced distortions. Proton induced distortions include decreasing proton and non-bonded oxide ion distances enhancing hydrogen bonding as seen in other studies [16,34]. Occasionally, the periodicity in this small system artificially allows octahedral tilting to be reversed across a plane. While this is not expected in the real system, this highlights some of the deficiencies of this small periodic model and at the same time shows that proton induced distortions are more significant when the dopant induced octahedral tilting is less. Considering larger systems and different dopant distributions would increase physical reality and allow further testing of current conclusions with an increase in computational effort. Rotational, intraoctahedral transfer and interoctahedral transfer transition states also roughly increase in energy with increasing proton-dopant distance for the 6.25% yttrium concentration system.
Fig. 8. The initial minimum (a), rotational transition state (b), and the final minimum (c) for a rotation that involves a change in backbone is shown.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
133
Fig. 9. The initial minimum (a), intraoctahedral transfer transition state (b), and the final minimum (c) for an intraoctahedral transfer with a smaller backbone change is shown.
However, as with the binding sites, the trend is not monotonic. Competition between dopant induced octahedral tilting and proton induced local distortions leads to energy diversity in these different states. Proton induced distortions can enhance hydrogen bonding for transition states as they did for binding sites. However, in addition, because the rotational transition state places the proton closer to the central barium, distortions moving the barium and hydroxide are also seen as in other perovskites [33]. The 12.5% concentration transition state energies are more clustered and show a clearer increasing energy trend as a function of increasing proton-dopant distance though oscillations still occur. A full map of all possible binding sites and transition states was not found at the lower concentration as the many proton induced distortions makes the size of such a map beyond our computational resources. Without the full map, it is not possible to find the best pathways or overall migration barriers. However, it is clear that the smaller dopant concentration leads to decreased octahedral tilting and can make proton induced distortions more significant. This leads to a greater diversity in binding and transition states though the range of barriers is only slightly expanded. As a result, it is likely that the limiting barrier is still of similar magnitude for the 6.25% yttrium doped system as for the 12.5% yttrium doped system though the distribution of paths in the lower concentration system should be more diverse. Our earlier work on the higher concentration system where we could make a full map gave an overall proton migration barrier of 0.39 eV [14] for a single proton system and 0.45 eV in an eight proton system [20]. Based on the similar range in activation barriers for the lower concentration system, we would expect similar overall migration barriers. In fact, experiments show that the activation barrier remains 0.43 − 0.44 eV for 2–10% and goes up to 0.48 eV by 20% yttrium doping [3]. Since the kinetic Monte Carlo calculations we have used in the past use zero temperature energies and harmonic vibrational frequencies in the pre-exponential factors,
we expect that the limiting migration barriers found are only valid for the lower temperature regime. However, at higher temperatures, more of the normal modes including significant octahedral tilting should become accessible and these may increase proton conduction through enhancing trap escape and maintaining protons in a trap free region lowering activation barriers at high temperatures as shown in experiments [7,10]. More detailed calculations of phonon modes including the proton and avoiding low temperature approximations could help understand these effects. 5. Conclusion Comparison of barium zirconate structures with 6.25% and 12.5% ordered yttrium doping at the zirconium site show increased octahedral tilting with increasing dopant concentration. The energy of proton binding sites and transition states roughly increases with proton-dopant distance. The greater diversity of possible proton binding sites in the lower concentration system leads to a greater spread in energies. The energy trends are affected by both octahedral tilting and proton induced distortions. Overall, the greater diversity of minima and transition states in the lower doped system suggests a more complex conduction mechanism in this system than in our earlier studies of the 12.5% yttrium doped system. Acknowledgments We would like to thank Frederick G. Haibach as well as anonymous reviewers for their insightful comments. This research was supported by an NSF RUI Grant No. CHE-1111474. Computational resources were provided in part by the MERCURY supercomputer consortium http://mars.hamilton.edu under NSF MRI No. CHE1229354.
Fig. 10. The initial minimum (a), interoctahedral transfer transition state (b), and the final minimum (c).
134
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
References [1] C. Duan, J. Tong, M. Shang, S. Nikodemski, M. Sanders, S. Ricote, A. Almansoori, R. O’Hayre, Readily processed protonic ceramic fuel cells with high performance at low temperatures, Science 349 (6354) (2017) 1321–1326. [2] K.D. Kreuer, Proton-conducting oxides, Annu. Rev. Mater. Res. 33 (2003) 333–359. [3] K.D. Kreuer, S. Adams, W. Münch, A. Fuchs, U. Klock, J. Maier, Proton conducting alkaline earth zirconates and titanates for high drain electrochemical applications, Solid Sate Ionics 145 (2001) 295. [4] F. Giannici, M. Shirpour, A. Longo, A. Martorana, R. Merkle, J. Maier, Long-range and short-range structure of proton-conducting Y:BaZrO3 , Chem. Mater. 23 (2011) 2994–3002. [5] A. Braun, A. Ovalle, S. Erat, V. Pomjakushin, A. Cervellino, W. Stolte, T. Graule, Yttrium and hydrogen superstructure and correlation of lattice expansion and proton conductivity in the BaZr0.9 Y0.1 O2.95 proton conductor, Appl. Phys. Lett. 95 (2009) 224103. [6] T. Schober, H.G. Bohn, Water vapor solubility and electrochemical characterization of the high temperature proton conductor BaZr0.9 Y0.1 O2.95 , Solid State Ionics 127 (3-4) (2000) 351–360. [7] Y. Yamazaki, F. Blanc, Y. Okuyama, L. Buannic, J.C. Lucio-Vega, C.P. Grey, S.M. Haile, Proton trapping in yttrium-doped barium zirconate, Nature Materials 12 (2013) 647–651. [8] H.G. Bohn, T. Schober, Electrical conductivity of high-temperature proton conductor BaZr0.9 Y0.1 O2.95 , J. Am. Ceram. Soc. 83 (4) (2000) 768–772. [9] E. Fabbri, D. Pergolesi, S. Licoccia, E. Traversa, Does the increase in Y-dopant concentration improve the proton conductivity of BaZr1 −x Yx O3 −d fuel cell electrolytes? Solid State Ionics 181 (2010) 1043–1051. [10] A. Braun, S. Duval, J. Embs, F. Juranyi, P. Ried, P. Holtappels, R. Hempelmann, U. Stimming, T. Graule, Proton diffusivity in the BaZr0.9 Y0.1 O3 −d proton conductor, J. Appl. Electr. Chem. 39 (4) (2009) 471–475. [11] J.W. Bennett, I. Grinberg, A.M. Rappe, Effect of symmetry lowering on the dielectric response of BaZrO3 , Phys. Rev. B 73 (2006) 180102. ´ J.D. Gale, Ground state structure of BaZrO3 : a comparative first[12] A. Bilic, principles study, Phys. Rev. B 79 (2009) 174107. [13] M.A. Gomez, M. Chunduru, L. Chigweshe, L. Foster, S.J. Fensin, K.M. Fletcher, L.E. Fernandez, The effect of yttrium dopant on the proton conduction pathways of BaZrO3 , a cubic perovskite, J. Chem. Phys. 132 (2010) 214709. [14] R. Krueger, F. Haibach, D. Fry, M. Gomez, Centrality measures highlight proton traps and access points to proton highways in kinetic Monte Carlo trajectories, J. Chem. Phys. 142 (2015) 154110. [15] M.A. Gomez, F.-J. Liu, Protons in Al doped BaZrO3 escape dopant traps to access long range proton conduction highways, Solid State Ionics 252 (2013) 40–47. [16] M.E. Björketun, P.G. Sundell, G. Wahnström, Effect of acceptor dopants on the proton mobility in BaZrO3 : a density functional investigation, Phys. Rev. B 76 (2007) 054307. [17] M.E. Björketun, P.G. Sundell, G. Wahnström, D. Engberg, A kinetic Monte Carlo study of proton diffusion in disordered perovskite structured lattices based on first-principles calculations, Solid State Ionics 176 (2005) 3035. [18] B. Merinov, W. Goddard, III, Proton diffusion pathways and rates in Y-doped BaZrO3 solid oxide electrolyte from quantum mechanics, J. Chem. Phys. 130 (2009) 194707.
[19] A.C.T. van Duin, B.V. Merinov, S.S. Han, C.O. Dorso, W.A. Goddard, III, ReaxFF reactive force field for Y-doped BaZrO3 proton conductor with applications to diffusion rates for multigranular systems, J. Chem. Phys. 112 (2008) 11414–11422. [20] M.A. Gomez, D. Fry, M. Sweet, Effects on the proton conduction limiting barriers and trajectories in BaZr0.875 Y0.125 O3 due to the presence of other protons, J. Korean Ceram. Soc. 53 (2016) 521–528. [21] A. Kuwabara, M. Shiga, C.A.J. Fisher, H. Moriwake, First principles calculations of solid solution states in Y-doped BaZrO3 , Solid State Proton Conductors Meeting, 2016. [22] A.M. Glazer, The classification of tilted octahedra in perovskites, Acta Crystallogr. Sect. B: Struct. Sci. 28 (1972) 3384. [23] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186. http://link.aps.org/doi/10.1103/PhysRevB.54.11169. http://dx.doi.org/10.1103/ PhysRevB.54.11169. [24] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (1999) 1758–1775. http://link.aps.org/doi/10.1103/PhysRevB.59.1758. http://dx.doi.org/10.1103/ PhysRevB.59.1758. [25] G. Kresse, Ab initio molekular dynamik für flüssige metalle, Ph.D. thesis. Technische Universität at Wien. 1993. [26] G. Kresse, J. Hafner, Ab initio molecular-dynamics for liquid-metals, Phys. Rev. B 47 (1993) RC558. [27] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1996) 15–50. [28] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169. [29] G. Kresse, J. Joubert, From ultrasoft pseudopotentials to the projector augmented wave method, Phys. Rev. B 59 (1999) 1758. [30] The labels on these PAW method files within the VASP library are Ba_sv, Zr_sv, Y_sv, Al, O, and H. [31] M.A. Gomez, M.A. Griffin, S. Jindal, K.D. Rule, V.R. Cooper, The effect of octahedral tilting on proton binding sites and in pseudo-cubic perovskite oxides, J. Chem. Phys. 123 (2005) 094703. [32] G. Henkelman, B.P. Uberuaga, H. Jonsson, A climbing image nudge elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (2000) 9901. [33] J.M. Polfus, T. Norby, R. Bredesen, Protons in oxysulfides, oxysulfates, and sulfides: a first-principles study of La2 O2 S, La2 O2 SO4 , SrZrS3 , and BaZrS3 , J. Phys. Chem. C 119 (42) (2015) 23875–23882. http://dx.doi.org/10.1021/acs. jpcc.5b08278. [34] F. Blanc, L. Sperrin, D. Lee, R. Dervisoglu, Y. Yamazaki, S.M. Haile, G. De Paëpe, C.P. Grey, Dynamic nuclear polarization NMR of low-gamma nuclei: structural insights into hydrated yttrium-doped BaZrO3 , J. Phys. Chem. Lett. 5 (2014) 2431–2436. [35] Y. Yamazaki, R. Hernandez-Sanches, S.M. Haile, High total proton conductivity in large-grained yttrium-doped barium zirconate, Chem. Mater. 21 (2009) 2755–2762.
Contents lists available at ScienceDirect
Solid State Ionics journal homepage: www.elsevier.com/locate/ssi
Ordered yttrium concentration effects on barium zirconate structure, proton binding sites and transition states Maria A. Gomez a, * , Gillian Kwan a , Wanshu Zhu a , Monica Chelliah a , Xintong Zuo b , Audrey Eshun a , Virginia Blackmer a , Truc Huynh a , Mai Huynh a a b
Department of Chemistry, Mount Holyoke College, South Hadley, MA 01075,United States Medical School, Duke-NUS, 169857, Singapore
A R T I C L E
I N F O
Article history: Received 30 October 2016 Received in revised form 22 March 2017 Accepted 27 March 2017 Available online xxxx
A B S T R A C T Acceptor doped perovskites show promise as fuel cell proton conducting ceramic membranes and continue to be rapidly developed. Proton conduction in barium zirconate with 6.25% ordered yttrium doping at the zirconium site is explored. The dopant increases octahedral distortions in the lattice similar to our earlier studied 12.5% dopant system though the angle distribution is slightly broader. Normal modes promoting vibrations about a non-transient octahedral tilt are found for both concentrations at frequencies thermally accessible at typical proton conduction temperatures. The energy of proton binding sites and transition states roughly increases with distance away from the dopant though there are many oscillations in the trend due to local proton induced distortions. Protons can both locally enhance and decrease the octahedral distortions induced by the dopant. Transition state energies roughly increase with proton-dopant distance. Rotational and intraoctahedral transfer transition state energies vary through a similar magnitude range though the intraoctahedral transition states are often higher in energy. Interoctahedral transfer transition states tend to have the highest energies. The diversity of minima and transition states suggests that proton conduction in 6.25% yttrium doped barium zirconate is more complex than conduction in 12.5% doped barium zirconate. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Acceptor doped perovskite oxides, especially yttrium doped barium zirconate, show increasing promise as fuel cell high to intermediate temperature proton conductors [1–3]. Central to improving performance of these materials is understanding how the energy landscape for proton conduction changes with dopant concentration. High proton concentration and fast conduction paths for individual protons lead to high conductivity. When barium zirconate is doped with yttrium, Y3+ occupies Zr4+ sites through a maximum load of 16% [3,4], forming charge-compensating oxygen vacancies. In humid conditions, the water can fill oxygen vacancies. Electron density shifts allow an oxide ion to remain in the former vacancy while two protons conduct through the perovskite moving from one electron rich oxide ion to another. Thus, the greater the yttrium content, the greater the expected proton concentration. However, since the yttrium dopant distorts barium zirconate effectively changing the
* Corresponding author. E-mail address: [email protected] (M.A. Gomez).
http://dx.doi.org/10.1016/j.ssi.2017.03.027 0167-2738/© 2016 Elsevier B.V. All rights reserved.
landscape in which protons conduct, increasing yttrium content may not lead to consistent increases in proton conduction. Many experiments have elucidated specific details of the proton conduction landscape. Braun et al. used X-ray and neutron diffraction studies to show that yttrium is arranged in an ordered superstructure in the barium zirconate lattice and noted that proton conduction occurs on a landscape set by the dopant [5]. Giannici et al. [4] proposed a landscape including octahedral tilting in light of the appearance of additional Raman peaks in the yttrium doped barium zirconate system and noted that these led to a permanent lowering from the cubic symmetry seen in previous studies [6]. Yamazaki et al. showed that to conduct, protons need to overcome a barrier to escape from the dopant and also a barrier to conduct once free from the dopant trap giving an overall migration barrier of 0.47 eV below 470 K [7]. Several other experiments have found migration barriers in this range for temperatures in the range 350–900 K [3,8,9]. Both Braun et al. and Yamazaki et al. have found that migration barriers decrease significantly at high temperatures suggesting a change in mechanism which lowers the activation energy significantly though they differ on the specific temperature critical to the change [7,10]. Characterizing structure flexibility and energetics may shed light on potential mechanism change at higher temperatures.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
Simulations see similar proton conduction landscape patterns and can shed further light on the specific proton conduction mechanisms in barium zirconate. Octahedral tilting has been seen to increase upon yttrium doping using Density Functional Theory (DFT) [11–13]. Very different proton conduction pathways have been found for different acceptor dopants that varied the extent of octahedral tilting [14,15]. Björketun et al. have found that larger dopants result in a larger trapping region in the conduction landscape decreasing energies of proton binding sites with high hydrogen bonding in second coordination shells. Using kinetic modelling, they also found that larger dopants had decreased activation barriers to conduction. A migration barrier of 0.30 eV was found for 12.5% yttrium doped barium zirconate in the 300–600 K temperature range [16,17]. In 12.5% Y doped BaZrO3 , Merinov and Goddard [18] used DFT to find activation energies of 0.48 to 0.49 eV and 0.41 eV for a variety of intraoctahedral and interoctahedral transfers, respectively. van Duin et al. [19] found an activation energy of 0.45 eV with a yttrium dopant concentration of 12.5% using molecular dynamics with a ReaxFF reactive force field. Using a combination of DFT, kinetic Monte Carlo, and graph theory, we have found that long range periodic pathways in the same system have limiting barriers of 0.3 eV [13] while kinetic Monte Carlo pathways which allow for the indirect meandering of diffusion have barriers of 0.39 eV [14] for a single proton system and 0.45 eV in an eight proton system at 1000 K [20]. These findings are in general agreement with experimental proton landscape details and low temperature migration barriers. Several experiments have highlighted changes in the conduction landscape with dopant concentration. A study by Kreuer et al. shows that proton mobility increases with increasing yttrium dopant concentration through a range of 2% to 20% though the mobility increase saturates by 10% doping. At the same time, the activation barrier remains 0.43–0.44 eV for 2–10 % and goes up to 0.48 eV by 20% in the temperature range of 313–625 K [3]. Hydration shows a similar pattern to mobility but starts to decrease from 20% to 25% dopant concentration [3]. Another study by Fabbri et al. [9] probed the yttrium concentration range of 20–60% in barium zirconate which is beyond the 16% concentration where yttrium occupies exclusively zirconium sites [4]. Proton concentration increased with dopant concentration and decreased with temperature [9] as expected and as seen in the earlier Kreuer et al. study [3]. Migration barriers [9] at the lower concentration range were also consistent with the earlier study [3]. This paper considers the effect of decreased yttrium concentration from 12.5% [13] to 6.25% in cases where there is a single acceptor dopant ordered through periodic boundary conditions. There are other groups [21] currently enumerating many possible multiple dopant and oxygen vacancy configurations. Such studies are critical as at least one experiment suggests that yttrium is arranged in an ordered super structure in barium zirconate [5]. In our prior work with 12.5% Y doped BaZrO3 [13], one in eight Zr sites contained yttrium in a periodic simulation box. As a result, every other octahedron in each direction contained a yttrium ion, a highly ordered acceptor dopant situation. In this study, the system is twice as large in one direction and again one Zr site is occupied by a yttrium ion which is then periodically replicated. Increasing size in all directions would lead to needing to consider multiple relative yttrium configurations and an overall significantly higher computational cost. Section 2 describes all the methods used in this contribution. Section 3.1 describes how lowering the dopant concentration from 12.5% to 6.25% affects the overall perovskite structure, tilting angle distribution, and low energy normal modes. Sections 3.2 and 3.3 describe the coupled effect of yttrium induced octahedral tilting and proton induced local distortions of the structure and energy of the binding sites and transition states, respectively. Finally, Sections 4 and 5 contain the discussion and conclusion.
127
2. Methods Bilic´ and Gale [12] show that even the structure of the undoped barium zirconate global energy minimum has alternating layers with opposite octahedral tilting in each direction corresponding to (−, −, −) in Glazer notation [22]. Bennet, Grinberg and Rappe also found octahedral tilting in the undoped barium zirconate structure [11]. In an earlier study [13], we found that the structure of 12.5% yttrium doped BaZrO3 showed greater octahedral distortions than the undoped system. In this study, we lower the yttrium concentration to 6.25%. To obtain the initial structure, we double the size of our earlier studied 2 × 2 × 2 unit cell system with 12.5% yttrium at the zirconium site and replaced one dopant ion with a zirconium ion thereby setting the doping level to 6.25%. The initial structure was in the lowest energy Glazer distortions [22] namely (−, −, −). Density functional theory (DFT) calculations with the PBE functional and a generalized gradient approximation (GGA) have been performed to find energies using the ab-initio total-energy and molecular-dynamics program VASP (Vienna ab-initio simulation program) developed at the Fakultät für Physik of the Universität Wien [23–28]. The projector augment wave (PAW) method [29] with the valence states of Ba(5s, 5p, 6s, 5d), Zr(4s, 4p, 5s, 4d), Y(4s, 4p, 5s, 4d), and O(2s, 2p) [30] is used. Partial occupancies are calculated with Gaussian smearing starting at 0.2 eV and automatically extrapolated to zero using VASP. Bands are optimized iteratively using the preconditioned residual minimization scheme. The first 10 iterations are non-self consistent. A 2×2×2 Monkhorst-Pack k-point mesh with zero shift is used. Accurate level precision and real space projections are done. The initial dimensions of the system were set to 2a×2a×4a where the a is the lattice constant of the undoped barium zirconate system of 4.26 Å found in our earlier study with the PBE functional [13]. This lattice constant is larger than the experimental lattice size of 4.19 Å [4] but similar to other DFT GGA VASP studies as seen in Table 1 of reference [12]. The LDA functional does match the experimental value better and gives similar structure patterns to the GGA functions [31] though GGA functionals are expected to do better with transition states. The conjugate gradient minimization method was used to find the optimum configuration of the nuclei in the 2a×2a×4a system with one yttrium dopant ion at a zirconium site. A uniform positively charged cloud is used to neutralize the charge. This uniform cloud naturally misses the non-uniform effect of localized charge compensating oxygen vacancies but is a common approximation in periodic systems when one does not want to enumerate all the possible configurations of two to one yttrium ion to oxygen vacancy ratios at each yttrium concentration. This structure was further allowed to relax through changing of both ionic positions and cell shape and volume. While the symmetry group was not explicitly set, VASP by default identifies the symmetry of the original structure and maintains it through the optimization. In this case, the minimal C1 symmetry is identified likely because earlier work [13] used random angle octahedral tilting in different directions to explore the different possible Glazer configurations [22]. Glazer does identify the (−, −, −) tilt system as having F 1¯ symmetry but we did not explicitly fit to or constrain for those symmetry elements. Further, as seen in Section 3.1, the length and angles of the simulation box suggest nearly cubic average unit cells. Once the lowest energy structure is found, normal modes are calculated. Our earlier 0% and 12.5% yttrium concentration system had fixed the lattice parameter at 4.26 Å [13] and box angles at 90◦ while the eight octahedra were allowed to tilt. To see the effect of dopant concentration on lattice size, the simulation box for these two concentrations is optimized with the same conditions as our earlier work [13] except that we allow lattice shape and size to relax in addition. Note, however, that all other comparisons with earlier data are comparing the earlier work with the 4.26 Å lattice constant [13].
128
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
Proton binding sites at 6.25% yttrium concentration are found using the conjugate gradient minimization method to relax ion positions but not cell shape or volume. Four initial proton locations are considered about each oxygen with the OH bond pointing into each cubic plane as in our earlier work [13]. Optimization is halted when forces are less than 0.01 eV/Å. Energies are expected to be accurate to about 0.02 eV. Relaxed transition state energies are found using a climbing nudged elastic band method [32]. Initially, a standard nudged elastic band method is used followed by a climbing nudged elastic band method where the highest image is set to climb. In an effort to find as many transition states as possible, the looser convergence criteria of converging forces at the transition state to 0.1 eV/Å was set. For some transition states, convergence in force was set to 0.01 eV/Å. The energy changes to the barriers are expected to have errors around 0.04 eV. 3. Results 3.1. Octahedral tilting in the optimized doped perovskite without protonation The final optimized simulation box dimensions for the 6.25% yttrium concentration case were 8.46 Å by 8.47 Å by 16.97 Å with angles 90.0◦ , 90.0◦ , and 90.1◦ for this 2×2×4 unit cell system. Table 1 shows the average lattice constants and angles for 0%, 6.25%, and 12.5% yttrium simulations where the box size and shape is allowed to change. The experimental lattice constant column shows individual lattice constants for cubic systems and two lattice constants for tetragonal systems. The calculated lattice constant values overestimate experimental values [3,4] though they do show the same increasing trend with dopant concentration. Notice that the 12.5% yttrium system shows a larger lattice constant when the unit cell is allowed to vary in shape and size than in our earlier work which held it at 4.26 Å [13]. The remainder of the 12.5% comparisons in this paper use the data from reference [13] which kept a box size that is closer to the experimental value at the higher concentration. Including oxygen vacancies and starting with different Glazer octahedral tilting can impact calculated lattice constants. Visual inspection of the structure shows a Glazer octahedral tilting pattern of (−, −, −). Fig. 1 compares one side of the optimized 6.25% doped structure with the earlier 12.5% doped structure [13] suggesting that octahedral tilting roughly decreases with decreasing dopant concentration. Octahedral tilting brings two oxide ions closer together while bringing two other oxide ions further apart. This results in a variation of O-O distances across cubes as shown in Fig. 1. The greater range in O-O distances in the higher doped system is a result of the larger octahedral tilting. For a more quantitative comparison, the distribution of ZrOZr and ZrOY angles is explored. In a perfectly cubic system without octahedral tilting these angles would be 180◦ . Fig. 2 shows the distribution of ZrOX angles with X=Zr and X=Y for both 12.5% and 6.25% yttrium doped systems. All angles are less than 180◦ highlighting that both Table 1 The average unit cell lengths are shown for optimized structures allowing the cell shape and size to relax. Experimental lattice constants are shown in the last column for comparison. The experiments with a single lattice constant found cubic symmetry while those with two lattice constants found tetragonal symmetry. %Y
Average lattice constants (Å)
Angles (◦ )
Experiment lattice constants (Å)
0 5 6.25 12.5 15
4.23, 4.22, 4.22
90.0, 90.1, 90.0
4.19 [4] 4.21 [3]
4.23, 4.24, 4.24 4.29, 4.29, 4.29
90.0, 90.0, 90.1 90.1, 90.1, 90.1 a=4.23, c=4.21 [3]
systems exhibit some octahedral tilting. The average ZrOZr angles are 170◦ and 172◦ with range widths of 6◦ and 10◦ for 12.5% and 6.25% yttrium, respectively. The average ZrOY angles are 164◦ and 165◦ with range widths of 1◦ and 8◦ for 12.5% and 6.25% yttrium, respectively. Both averages show a decrease in octahedral tilting or an increase in angle as the yttrium concentration decreases. The details of the distributions reveal more significant subtleties than the averages. The ZrOZr angle distribution in the lower yttrium concentration system has more larger angles i.e. angles showing less tilt or distortion from the cubic case. However, there is a small number of angles that show increased tilting. The lower yttrium concentration ZrOY angle distribution shows two smaller angles and two larger angles with an overall higher average angle showing overall decreased tilting. It is important to note that this 2 × 2 × 4 system is only lower in concentration in one direction. The contrasting higher concentration in the perpendicular plane may be leading to the increased tilting peaks. An earlier experimental study [4] noted that EXAFS data could be described by a ZrOY angle of about 160◦ . Both concentrations considered here have most ZrOY angles in the range of 160–165◦ . O-Zr distances are 2.13(3) Å and 2.12(2) Å at 12.5% and 6.25% doping, respectively. O-Y distances are 2.23(1) Å and 2.22(1) Å at 12.5% and 6.25% doping, respectively. The value in parenthesis is one standard deviation of the distribution at the last significant figure. XRD and EXAFS data by Giannici [4] in Fig. 10 shows distances within a standard deviation of these values. It is worthwhile noting that the tilt system used to fit the data in Giannici [4] is (0, 0, +) while our earlier work [13] predicted a (−, −, −) tilt system for 12.5% doping and had the (0, 0, +) case as 0.17 eV higher in energy. The (−, −, −) tilt system was also found in the undoped case in another study [12]. The zirconium and yttrium ions seem to stay very close to the center of the octahedra in the minimum energy structures. Normal mode analysis of the minimum energy structures shows many low energy vibrations involving zirconium and barium oscillations and higher energy octahedral tilting vibrations. One imaginary mode with a frequency corresponding to an energy less that 0.01 eV was found for the low and high dopant concentrations. Since the energy of this mode is less than the error expected in these calculations and moving in the mode direction did not yield a lower energy, the structure was considered to be a true minimum. The lowest energy normal mode at 6.25% and 12.5% yttrium concentrations has frequencies of 38 cm −1 and 40 cm −1 , respectively. These frequencies are again lower than the accuracy of the energy calculations and correspond to 55 K and 58 K. Low frequency modes involve vibrations of the barium ions and some zirconium ions within their octahedra. Modes at 200 cm −1 and 199 cm −1 for 6.25% and 12.5% yttrium concentrations, respectively, involve ZrOX angle changes and are indicative of octahedral tilting. These higher frequencies correspond to temperatures of 288 K and 286 K. Other higher frequency normal modes also involve tilting. This suggests that at typical proton conduction temperatures that low frequency vibrational modes will involve moving of cations in their sites and higher frequency though still thermally accessible modes will show octahedral tilting oscillations. Higher frequency modes (up to 787 cm −1 and 763 cm −1 in low and high dopant systems) show even more octahedral tilting vibrations but are less thermally accessible. Overall, even at moderate temperatures, there is enough thermal energy that the structure is flexible and small normal mode oscillations around the lowest energy structure are possible. 3.2. Structure and trends for proton binding sites As shown in Fig. 3, the relative energy for proton binding sites (black squares) roughly increases as a function of distance away from the dopant up to about 5.5 Å. The energy shown is relative to the lowest energy proton binding site. The trend is not monotonic due to proton induced local distortions which can change the proximity of
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
129
Fig. 1. (a) shows a section of the 6.25% yttrium doped barium zirconate system while (b) shows the 12.5% doped system. Oxide ions are shown in red. Barium ions are shown in pink. Octahedra around zirconium ions are in green and octahedra around yttrium ions are in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the proton and non-bonded oxide ions beyond the dopant induced tilting. We will refer to the decreasing of the distance between the proton and non-bonded oxide ions as increasing hydrogen bonding. Fig. 4 shows the corresponding data from our earlier work [13] on 12.5% doping to be clustered. The trend is for increasing relative energy as a function of proton distance from the dopant. The one exception to the increase in energy occurs at slightly under 4.0 Å where local distortions bring some binding sites on oxide ions that are second nearest neighbors to the dopant to a lower energy. The low energy for protons on some first and some second nearest neighbor oxide ions are more spread out in the 6.25% yttrium case than the 12.5% case. The increased octahedral tilting in the 12.5% yttrium system results in a reduction of binding site diversity. Fig. 5 shows three binding sites around the same oxide ion for the 6.25% yttrium case. In each case, the proton induces significant local distortions. In binding site (a), the proton brings two oxide ions which were 4.35 Å apart in the unprotonated system to a distance of 3.57 Å apart. This allows both oxide ions to be closer to the proton.
Since all octahedra in the system are connected by vertices and the long side of this simulation box is just four octahedra long, this proton induced distortion reverses the dopant induced distortion in this configuration. This is an artificial effect of the boundary conditions. Releasing the periodic constraints by freezing the plane of octahedra in the center of this simulation box would have prohibited long range tilting changes. However, the artificial system result highlights that when the dopant concentration is low, local proton induced tilting can be significant. In binding site (b), the proton increases the original dopant induced tilting bringing two oxide ions from 4.07 Å apart to 3.51 Å apart. Again, this brings both oxide ions closer to the proton. Binding site (c) shows that the proton can also bring the oxide ion in the direction of the OH bond closer. The distance between the oxygen in the hydroxide and the opposite oxide ion started at 4.05 Å and contracted to 3.55 Å. In each of these cases, the proton brings oxide ions closer enhancing hydrogen bonding. Overall, protons can induce oxide ions to approach each other or move apart through local tilting of octahedra around zirconium
12 12.5% X = Zr 12.5 % X = Y 6.25% X = Zr 6.25% X = Y
Number of Observations
10
8
6
4
2
0 160
165
170
175
180
Fig. 2. The distribution of ZrOX angles is shown for both 12.5% and 6.25% yttrium concentrations. Traditional histogram rectangles are used for the higher concentration while points are used for the lower concentration. ZrOZr angle distributions are shown in black while ZrOY angle distributions are shown in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
130
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
0.90 Inter TS Intra TS Rot TS Minima
0.80
0.70
Relative Energy (eV)
0.60
0.50
0.40
0.30
0.20
0.10
0.00 2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Distance from dopant (Å) Fig. 3. The relative energy of minima and transition states is shown as a function of proton-dopant distance for 6.25% dopant concentration. All energies are relative to the global minimum for this concentration. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
and yttrium ions. The reversal in the tilting induced by the dopant in some cases is likely just an artifact of small system size and periodic boundary conditions. Figs. 6 and 7 highlight that the ZrOX distribution with X=Zr and Y can widen significantly upon addition of a proton in both the 6.25% and 12.5% yttrium doped systems. The three binding sites considered in the higher dopant concentration case have a hydroxide in a plane without dopant. While both distributions broaden relative to the no proton cases, the lower concentration ZrOX distribution has lower octahedral tilting or angles closer to 180◦ for the binding sites far from the dopant.
3.3. Transition states connecting proton binding sites Fig. 3 shows the energy of all the transition states between proton binding sites relative to the binding site global minimum. In general, all transition state energies roughly increase as the proton-dopant distance increases until about 5.5 Å. Transition states energies for rotations (green plus symbols) and intraoctahedral transfers (red triangles) are similar with the latter more often being higher in energy. The highest transition state energies are for interoctahedral transfers (blue circles). As we saw in Section 3.2, the backbone changes due to different local distortions in the different minima.
0.90 Inter TS Intra TS Rot TS Minima
0.80
0.70
Relative Energy (eV)
0.60
0.50
0.40
0.30
0.20
0.10
0.00 2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
Distance from dopant (Å) Fig. 4. The relative energy of minima and transition states is shown as a function of proton-dopant distance for 12.5% dopant concentration. All energies are relative to the global minimum for this concentration.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
131
Fig. 5. Three different binding site structures on the same oxide ion are shown in (a), (b), and (c).
Not all possible backbone distortions were considered to find the transition states so our data represents only a subset of possible minimum energy configurations and transition states. Fig. 4 shows the corresponding transition states for the 12.5% dopant concentration system. The overall pattern is the same but the energies for binding and transition states are more localized in the higher concentration system. A variety of rotational, intraoctahedral transfer, and interoctahedral transfer transitions are seen in the 6.25% yttrium concentration system. A small sampling is presented in Figs. 8, 9, and 10. Fig. 8 highlights a rotation including a backbone distortion. The distance between two oxide ions is shown to highlight the face in which the backbone change is most prominently seen. Notice that the two oxide ions tracked go from being the closer two oxide ions on that face to being the further two oxide ions on that cube face. At the same time, the dopant-proton distance increases from 2.32 Å to 3.06 Å to 3.38 Å while the proton-barium distance goes from 2.54 Å to 2.42 Å to 2.50 Å. Rotations tend to bring the proton in close proximity to barium at the transition state. To try to overcome this repulsion, the system distorts allowing both hydroxide and barium to move away as in other studies [33]. These two oxide ions near the proton in (a), (b), and (c) with average distances to the proton of 2.17 Å, 2.58 Å, and 2.15 Å, respectively. The array of competing factors makes it difficult to predict whether the forward or backward barrier is higher for this transition. Both forward direction increases in dopant-proton distance as well as sharper decreases in barium-proton distance
suggest a higher forward barrier. The average hydrogen bond length, however, has a larger increase in the backward direction suggesting a higher backward barrier. No doubt, there are other factors not yet considered as well. The actual forward barrier of 0.21 eV is slightly lower than the reverse barrier of 0.26 eV. Fig. 9 shows an intraoctahedral transfer that retains the same basic backbone with smaller proton induced changes. This transfer occurs on a face with higher dopant concentration and has a larger barrier to leave the dopant than to reach it namely 0.32 and 0.25 eV, respectively. Fig. 10 shows an interoctahedral transfer in the higher concentration plane and near the dopant. Forward and backward barriers are 0.29 and 0.23 eV, respectively, with the barrier for leaving the dopant (forward) being slightly larger. The broad change in relative energy of both minima and transition states as the proton/dopant distance increases seen in Fig. 3 for the 6.25% yttrium concentration case leads to a very broad pattern for barriers to moving closer and further from the dopant. Individual cases highlighted show that while sometimes the energy decreases as the dopant is approached, there are local proton induced changes that may counter this effect. This can be a result of increased hydrogen bonding or decreased proton-barium repulsion. At 12.5% yttrium concentration, Fig. 4 shows significantly more clustering of binding sites and transition state energies leading to clearer patterns in activation barriers as the proton-dopant distance increases. The patterns are still prone to oscillations due to local distortions caused by hydrogen bond enhancement and proton-barium repulsion.
12 no protons X = Zr no protons X = Y (a) X = Zr (b) X = Zr (c) X = Zr (a) X=Y (b) X = Y (c) X = Y
Number of Observations
10
8
6
4
2
0 155
160
165
170
175
180
Fig. 6. The distribution of ZrOX with X=Zr and Y for Fig. 5 (a), (b), and (c) proton binding site structures are shown.
132
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134 12 no protons X = Zr no protons X = Y (a) X = Zr (b) X = Zr (c) X = Zr (a) X = Y (b) X = Y (c) X = Y
Number of Observations
10
8
6
4
2
0 155
160
165
170
175
180
Fig. 7. The distribution of ZrOX with X=Zr and Y for three binding sites (a), (b) and (c) on an oxide ion on a plane between the dopant planes in the 12.5% yttrium doped system.
4. Discussion The optimized 6.25% yttrium doped barium zirconate structure shows slight octahedral tilting while still retaining roughly cubic average unit cells. Visual comparison with the 12.5% yttrium system shows that decreasing the yttrium concentration decreases octahedral tilting. Comparison of distributions of ZrOX angle with X=Zr and Y between 6.25% and 12.5% yttrium doped barium zirconate shows roughly increased tilting with increasing yttrium concentration. Low energy normal modes with zirconium and barium ions oscillating in their lattice positions are seen in addition to higher energy normal modes with some octahedral tilting. Some of these octahedral tilting normal modes have thermally accessible energies at proton conduction temperatures suggesting that small octahedral tilting oscillations maintaining overall tilt direction are possible during proton conduction. Octahedral tilting vibrations at higher frequencies have been calculated before for undoped barium zirconate by the Rappe group [11]. Since there is more octahedral tilting in the doped systems, a lowering in the frequency of these vibrations is not unexpected. Transient octahedral distortions were suggested for the undoped system in addition to non-transient distortions for doped systems by Giannici et al. [4]. ZrOY angles, and O-Zr and O-Y distances are in also in good agreement with experiment [4].
The energy of proton minima roughly increases as the proton is moved away from the dopant. The increased diversity of backbone shifts leads to an increase in the spread of minimum energies in the 6.25% dopant case in comparison to the more clustered 12.5% dopant energies. Both show low energies for protons bound to oxide ions in both the first and second neighbor shells of the yttrium as seen in other studies [16,34,35]. Overall, energy trends are not only affected by the dopant induced octahedral tilting but also by proton induced distortions. Proton induced distortions include decreasing proton and non-bonded oxide ion distances enhancing hydrogen bonding as seen in other studies [16,34]. Occasionally, the periodicity in this small system artificially allows octahedral tilting to be reversed across a plane. While this is not expected in the real system, this highlights some of the deficiencies of this small periodic model and at the same time shows that proton induced distortions are more significant when the dopant induced octahedral tilting is less. Considering larger systems and different dopant distributions would increase physical reality and allow further testing of current conclusions with an increase in computational effort. Rotational, intraoctahedral transfer and interoctahedral transfer transition states also roughly increase in energy with increasing proton-dopant distance for the 6.25% yttrium concentration system.
Fig. 8. The initial minimum (a), rotational transition state (b), and the final minimum (c) for a rotation that involves a change in backbone is shown.
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
133
Fig. 9. The initial minimum (a), intraoctahedral transfer transition state (b), and the final minimum (c) for an intraoctahedral transfer with a smaller backbone change is shown.
However, as with the binding sites, the trend is not monotonic. Competition between dopant induced octahedral tilting and proton induced local distortions leads to energy diversity in these different states. Proton induced distortions can enhance hydrogen bonding for transition states as they did for binding sites. However, in addition, because the rotational transition state places the proton closer to the central barium, distortions moving the barium and hydroxide are also seen as in other perovskites [33]. The 12.5% concentration transition state energies are more clustered and show a clearer increasing energy trend as a function of increasing proton-dopant distance though oscillations still occur. A full map of all possible binding sites and transition states was not found at the lower concentration as the many proton induced distortions makes the size of such a map beyond our computational resources. Without the full map, it is not possible to find the best pathways or overall migration barriers. However, it is clear that the smaller dopant concentration leads to decreased octahedral tilting and can make proton induced distortions more significant. This leads to a greater diversity in binding and transition states though the range of barriers is only slightly expanded. As a result, it is likely that the limiting barrier is still of similar magnitude for the 6.25% yttrium doped system as for the 12.5% yttrium doped system though the distribution of paths in the lower concentration system should be more diverse. Our earlier work on the higher concentration system where we could make a full map gave an overall proton migration barrier of 0.39 eV [14] for a single proton system and 0.45 eV in an eight proton system [20]. Based on the similar range in activation barriers for the lower concentration system, we would expect similar overall migration barriers. In fact, experiments show that the activation barrier remains 0.43 − 0.44 eV for 2–10% and goes up to 0.48 eV by 20% yttrium doping [3]. Since the kinetic Monte Carlo calculations we have used in the past use zero temperature energies and harmonic vibrational frequencies in the pre-exponential factors,
we expect that the limiting migration barriers found are only valid for the lower temperature regime. However, at higher temperatures, more of the normal modes including significant octahedral tilting should become accessible and these may increase proton conduction through enhancing trap escape and maintaining protons in a trap free region lowering activation barriers at high temperatures as shown in experiments [7,10]. More detailed calculations of phonon modes including the proton and avoiding low temperature approximations could help understand these effects. 5. Conclusion Comparison of barium zirconate structures with 6.25% and 12.5% ordered yttrium doping at the zirconium site show increased octahedral tilting with increasing dopant concentration. The energy of proton binding sites and transition states roughly increases with proton-dopant distance. The greater diversity of possible proton binding sites in the lower concentration system leads to a greater spread in energies. The energy trends are affected by both octahedral tilting and proton induced distortions. Overall, the greater diversity of minima and transition states in the lower doped system suggests a more complex conduction mechanism in this system than in our earlier studies of the 12.5% yttrium doped system. Acknowledgments We would like to thank Frederick G. Haibach as well as anonymous reviewers for their insightful comments. This research was supported by an NSF RUI Grant No. CHE-1111474. Computational resources were provided in part by the MERCURY supercomputer consortium http://mars.hamilton.edu under NSF MRI No. CHE1229354.
Fig. 10. The initial minimum (a), interoctahedral transfer transition state (b), and the final minimum (c).
134
M. Gomez et al. / Solid State Ionics 304 (2016) 126–134
References [1] C. Duan, J. Tong, M. Shang, S. Nikodemski, M. Sanders, S. Ricote, A. Almansoori, R. O’Hayre, Readily processed protonic ceramic fuel cells with high performance at low temperatures, Science 349 (6354) (2017) 1321–1326. [2] K.D. Kreuer, Proton-conducting oxides, Annu. Rev. Mater. Res. 33 (2003) 333–359. [3] K.D. Kreuer, S. Adams, W. Münch, A. Fuchs, U. Klock, J. Maier, Proton conducting alkaline earth zirconates and titanates for high drain electrochemical applications, Solid Sate Ionics 145 (2001) 295. [4] F. Giannici, M. Shirpour, A. Longo, A. Martorana, R. Merkle, J. Maier, Long-range and short-range structure of proton-conducting Y:BaZrO3 , Chem. Mater. 23 (2011) 2994–3002. [5] A. Braun, A. Ovalle, S. Erat, V. Pomjakushin, A. Cervellino, W. Stolte, T. Graule, Yttrium and hydrogen superstructure and correlation of lattice expansion and proton conductivity in the BaZr0.9 Y0.1 O2.95 proton conductor, Appl. Phys. Lett. 95 (2009) 224103. [6] T. Schober, H.G. Bohn, Water vapor solubility and electrochemical characterization of the high temperature proton conductor BaZr0.9 Y0.1 O2.95 , Solid State Ionics 127 (3-4) (2000) 351–360. [7] Y. Yamazaki, F. Blanc, Y. Okuyama, L. Buannic, J.C. Lucio-Vega, C.P. Grey, S.M. Haile, Proton trapping in yttrium-doped barium zirconate, Nature Materials 12 (2013) 647–651. [8] H.G. Bohn, T. Schober, Electrical conductivity of high-temperature proton conductor BaZr0.9 Y0.1 O2.95 , J. Am. Ceram. Soc. 83 (4) (2000) 768–772. [9] E. Fabbri, D. Pergolesi, S. Licoccia, E. Traversa, Does the increase in Y-dopant concentration improve the proton conductivity of BaZr1 −x Yx O3 −d fuel cell electrolytes? Solid State Ionics 181 (2010) 1043–1051. [10] A. Braun, S. Duval, J. Embs, F. Juranyi, P. Ried, P. Holtappels, R. Hempelmann, U. Stimming, T. Graule, Proton diffusivity in the BaZr0.9 Y0.1 O3 −d proton conductor, J. Appl. Electr. Chem. 39 (4) (2009) 471–475. [11] J.W. Bennett, I. Grinberg, A.M. Rappe, Effect of symmetry lowering on the dielectric response of BaZrO3 , Phys. Rev. B 73 (2006) 180102. ´ J.D. Gale, Ground state structure of BaZrO3 : a comparative first[12] A. Bilic, principles study, Phys. Rev. B 79 (2009) 174107. [13] M.A. Gomez, M. Chunduru, L. Chigweshe, L. Foster, S.J. Fensin, K.M. Fletcher, L.E. Fernandez, The effect of yttrium dopant on the proton conduction pathways of BaZrO3 , a cubic perovskite, J. Chem. Phys. 132 (2010) 214709. [14] R. Krueger, F. Haibach, D. Fry, M. Gomez, Centrality measures highlight proton traps and access points to proton highways in kinetic Monte Carlo trajectories, J. Chem. Phys. 142 (2015) 154110. [15] M.A. Gomez, F.-J. Liu, Protons in Al doped BaZrO3 escape dopant traps to access long range proton conduction highways, Solid State Ionics 252 (2013) 40–47. [16] M.E. Björketun, P.G. Sundell, G. Wahnström, Effect of acceptor dopants on the proton mobility in BaZrO3 : a density functional investigation, Phys. Rev. B 76 (2007) 054307. [17] M.E. Björketun, P.G. Sundell, G. Wahnström, D. Engberg, A kinetic Monte Carlo study of proton diffusion in disordered perovskite structured lattices based on first-principles calculations, Solid State Ionics 176 (2005) 3035. [18] B. Merinov, W. Goddard, III, Proton diffusion pathways and rates in Y-doped BaZrO3 solid oxide electrolyte from quantum mechanics, J. Chem. Phys. 130 (2009) 194707.
[19] A.C.T. van Duin, B.V. Merinov, S.S. Han, C.O. Dorso, W.A. Goddard, III, ReaxFF reactive force field for Y-doped BaZrO3 proton conductor with applications to diffusion rates for multigranular systems, J. Chem. Phys. 112 (2008) 11414–11422. [20] M.A. Gomez, D. Fry, M. Sweet, Effects on the proton conduction limiting barriers and trajectories in BaZr0.875 Y0.125 O3 due to the presence of other protons, J. Korean Ceram. Soc. 53 (2016) 521–528. [21] A. Kuwabara, M. Shiga, C.A.J. Fisher, H. Moriwake, First principles calculations of solid solution states in Y-doped BaZrO3 , Solid State Proton Conductors Meeting, 2016. [22] A.M. Glazer, The classification of tilted octahedra in perovskites, Acta Crystallogr. Sect. B: Struct. Sci. 28 (1972) 3384. [23] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186. http://link.aps.org/doi/10.1103/PhysRevB.54.11169. http://dx.doi.org/10.1103/ PhysRevB.54.11169. [24] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (1999) 1758–1775. http://link.aps.org/doi/10.1103/PhysRevB.59.1758. http://dx.doi.org/10.1103/ PhysRevB.59.1758. [25] G. Kresse, Ab initio molekular dynamik für flüssige metalle, Ph.D. thesis. Technische Universität at Wien. 1993. [26] G. Kresse, J. Hafner, Ab initio molecular-dynamics for liquid-metals, Phys. Rev. B 47 (1993) RC558. [27] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1996) 15–50. [28] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169. [29] G. Kresse, J. Joubert, From ultrasoft pseudopotentials to the projector augmented wave method, Phys. Rev. B 59 (1999) 1758. [30] The labels on these PAW method files within the VASP library are Ba_sv, Zr_sv, Y_sv, Al, O, and H. [31] M.A. Gomez, M.A. Griffin, S. Jindal, K.D. Rule, V.R. Cooper, The effect of octahedral tilting on proton binding sites and in pseudo-cubic perovskite oxides, J. Chem. Phys. 123 (2005) 094703. [32] G. Henkelman, B.P. Uberuaga, H. Jonsson, A climbing image nudge elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (2000) 9901. [33] J.M. Polfus, T. Norby, R. Bredesen, Protons in oxysulfides, oxysulfates, and sulfides: a first-principles study of La2 O2 S, La2 O2 SO4 , SrZrS3 , and BaZrS3 , J. Phys. Chem. C 119 (42) (2015) 23875–23882. http://dx.doi.org/10.1021/acs. jpcc.5b08278. [34] F. Blanc, L. Sperrin, D. Lee, R. Dervisoglu, Y. Yamazaki, S.M. Haile, G. De Paëpe, C.P. Grey, Dynamic nuclear polarization NMR of low-gamma nuclei: structural insights into hydrated yttrium-doped BaZrO3 , J. Phys. Chem. Lett. 5 (2014) 2431–2436. [35] Y. Yamazaki, R. Hernandez-Sanches, S.M. Haile, High total proton conductivity in large-grained yttrium-doped barium zirconate, Chem. Mater. 21 (2009) 2755–2762.