- Email: [email protected]
Effect of strain on the oxygen diffusion in yttria and gadolinia co-doped ceria
Solid State Ionics 230 (2013) 37–42
Contents lists available at SciVerse ScienceDirect
Solid State Ionics journal homepage: www.elsevier.com/locate/ssi
Effect of strain on the oxygen diffusion in yttria and gadolinia co-doped ceria M.J.D. Rushton b, A. Chroneos a,⁎, S.J. Skinner b, J.A. Kilner b, R.W. Grimes b a b
Materials Engineering, The Open University, Milton Keynes, MK7 6AA, United Kingdom Department of Materials, Imperial College London, London SW7 2AZ, United Kingdom
a r t i c l e
i n f o
Article history: Received 23 May 2012 Received in revised form 25 August 2012 Accepted 12 September 2012 Available online 6 October 2012 Keywords: Doped ceria Molecular dynamics Strain
a b s t r a c t Atomic scale simulations have been used to investigate the impact of co-doping (yttrium and gadolinium) and strain on oxygen diffusion and binding of dopant–vacancy clusters in ceria. Doped ceria in its relaxed or strained form is used as an electrolyte for solid oxide fuel cell applications. For unstrained co-doped ceria we calculate an activation energy for migration of 0.70–0.75 eV in the temperature range of 973–1873 K. Co-doping with yttrium and gadolinium only affected oxide ion diffusion to a small degree when compared to single doping. The diffusion coefficient was substantially increased by tensile strain while compressive strain caused a decrease. To gain further insight why tensile strain leads to higher diffusivity static simulations were employed. It is calculated that tensile strain reduces the binding energies of clusters between oxygen vacancies and trivalent dopant atoms while compressive strain leads to higher binding energies. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Solid oxide fuel cells (SOFCs) face two major obstacles to a large scale application. These are the high cost of current SOFC designs and their degradation over the long timescales needed for commercial application. One route for SOFCs to overcome these barriers is through operation within the intermediate temperature range (500–700 °C) [1]. This presents difficulties for traditional ceramic electrodes and electrolytes as it requires materials with enhanced activity, which leads to a requirement for higher oxygen diffusivities [2]. One way of achieving this is by identifying new materials and a number of candidate oxides are being investigated for cathode and electrolyte applications in the next generation intermediate temperature SOFCs (IT-SOFCs) [3–8]. An alternative way of enhancing ionic transport may be through the presence of an interface between dissimilar oxides, since the effects of interfacial strain are thought to affect oxide ion transport [9–15]. More conventional methods for enhancing conductivity involve co-doping, the idea being to minimise the effect of any defect interactions through the introduction of heterogeneous dopant species [16]. Trivalent-doped ceria is a well established electrolyte material due to its high oxygen ion diffusivity and relatively lower reduction temperatures [17]. Doping with oxides such as R2O3 (here R = Gd and/or Y) improves oxygen diffusion due to the formation of oxygen vacancies in the resultant Ce1 −xRxO2−x/2 solid solution [18]. Doped ceria has received considerable attention due to its application as the electrolyte material in SOFCs as it exhibits up to 2–3 orders of magnitude higher oxygen ion conductivity than the common high temperature fuel cell
electrolyte, and yttria stabilised zirconia, at intermediate temperatures [19–21]. Co-doping may offer an attractive method to increase the conductivity in CeO2 even further [22–25]. Atomic scale simulations have been carried out in order to gain insight into the effects of strain and/or doping regimes on oxygen transport in SOFC materials [26–30]. Recently De Souza et al. considered the effect of strain on oxygen diffusivity within ceria using a static lattice method and predicted a considerable reduction in the activation energy for oxygen vacancy migration on dilation of the lattice [29]. Here, we additionally consider the combined effects of single and co-doping on strained and relaxed lattices. In a similar manner to De Souza's study, static Mott–Littleton calculations are used to consider the interaction of point defects within strained ceria, furthermore we use molecular dynamics (MD) calculations to observe how strain affects oxygen ion diffusivity as a function of temperature. 2. Methodology Both the MD and the static Mott–Littleton simulation methodologies employed in this study rely on a simple description of interatomic forces. Here the Born description of ionic crystals was used [31]: forces acting between atoms are resolved into long-range coulombic and shortrange pair-wise interactions described using parameterized Buckingham pair potentials [32]. The lattice energy, EL, is given by,
EL ¼
N −1 X N X i¼1 j>i
⁎ Corresponding author. E-mail addresses: [email protected] (M.J.D. Rushton), [email protected] (A. Chroneos). 0167-2738/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ssi.2012.09.015
! qi qj r ij C ij þ Aij exp − − 6 4πε0 r ij ρij r ij
ð1Þ
where N is the number of ions in the system, ε0 is the permittivity of free space, rij is the separation between ions i and j of charges qi and qj
38
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
Table 1 Atomic composition of the simulated systems with the chemical formula Ce0.835X0.165O1.96. X
Number of atoms
Gd:Y
Gd
Y
Ce
O
1:0 2:1 1:1 1:2 0:1
338 225 169 113 0
0 113 169 225 338
1710 1710 1710 1710 1710
3927 3927 3927 3927 3927
dielectric continuum approximation, with the ionic relaxations taking place in a single step. The interaction energies between ions in Regions I and IIa were calculated explicitly. Finally, the outer Region IIb consists of a point charge array providing the Madelung field of the remaining lattice, and the response energy in Region IIb was calculated using the Mott–Littleton approximation [36]. The GULP code [37] was used for all the static atomistic simulations. The interaction of point defects within the lattice can be quantified by the calculation of their binding energies, Eb, defined by Eb ¼ ∑Eisolateddefects −Edefectcluster :
respectively, and Aij, ρij and Cij are the short-range Buckingham potential parameters specific to particular pairs of interacting species. The potential parameters have already been employed successfully in static simulations of doped ceria [27]. Here they are used without the shell model of the previous work to facilitate computationally tractable MD simulations. A short-range potential cut-off of 6.5 Å was used. As full periodic boundaries [33] were employed, long range Coulombic forces were calculated using the Ewald sum [34,35].
ð2Þ
According to this definition a positive binding energy implies that the defect cluster is more energetically favourable as compared to the isolated defects. 2.2. Molecular dynamics
Static calculations were used to determine the binding energies of point defects within the ceria lattice. The response of the crystal to charged defects was determined using the Mott–Littleton approach. Within this, a perfect crystal lattice is initially constructed: periodic boundaries were employed whereby a unit-cell is tessellated in space to represent a continuous crystal. The energy of the system is minimised by relaxing atomic positions and cell parameters via a conjugate gradient approach. The relaxed structure is then used to create concentric spherical regions. Point defects are introduced at the centre of the innermost region (Region I). Here the ions are treated explicitly and are relaxed subject to forces calculated from the interionic potentials. In Region IIa, which surrounds Region I, the ionic forces are calculated via a
MD simulations were used to investigate the effects on oxygen ion diffusion of co-doping ceria with gadolinia and yttria, as a function of temperature (between 973 K and 1873 K). In addition strain in the range of −2.5 to 2.5% was applied to doped compositions in order to examine its influence on diffusion. During MD simulations, Newton's second law of motion is solved iteratively. This means that the positions and velocities of a set of atoms can evolve as a function of time from a description of interatomic forces [33,38]. Here the velocity Verlet integration method was used [35,38] and an integration time step of 1 fs was used throughout. Molecular dynamics calculations were performed using the DL_POLY code (version 3.09) [35,39]. Each simulated system contained CeO2 doped with 9 mol% of X2O3 (where X could be Gd3+ and/or Y 3+). For the constant pressure (i.e. zero strain) simulations, five Gd:Y ratios were considered ranging from doping only with Gd (1:0), through co-doped compositions with ratios of 2:1, 1:1, 1:2, to 0:1 (the last corresponding to doping only
Fig. 1. Arrhenius plot of the calculated oxygen diffusion coefficient (D) against temperature (T) for 9 mol% X2O3 in CeO2 for X = Gd3+ and/or Y3+ with ratios as indicated in the figure. Error bars indicate the highest and lowest D values for the five random configurations considered while points show the mean. For clarity, at each temperature points have been spread along the x-axis with grey drop lines representing the temperature each group actually represents.
Fig. 2. Oxygen diffusion coefficient (D) against temperature (T) for 9 mol% X2O3 in CeO2 for X = Gd3+ and/or Y3+ with ratios as indicated in the figure. Compared with experimentally determined diffusion coefficients for singly doped CeO2 compositions. For the purposes of comparison, the 9 mol% composition studied here, has the chemical formula: Ce0.835X0.165O1.918.
2.1. Static atomistic simulations
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
T (K) 10−5
1273
1473
1873
The defect equation (in Kröger Vink notation) [40] for this process is as follows:
1073
a) Gd:Y = 1:0
CeO2
X 2 O3 → 2X
D (cm2s−1)
10−6
10−7
10−8
10−9 6
7
8
9
10
10000/T (K−1) T (K) 10−5
1273
1473
1873
1073
b) Gd:Y = 1:1
D (cm2s−1)
10−6
10−7
10−8
7
8
9
10
′
••
Ce
X
þ VO þ 3OO :
ð3Þ
Every two trivalent ions introduced into the cation sublattice lead to the formation of a charge compensating oxygen vacancy. With this in mind, each simulation cell was obtained from an 8 × 8 × 8 CeO2 super-cell (originally containing 6144 atoms). To this 338 dopant ions were randomly assigned to the Ce sites in the desired Gd:Y ratio while 169 oxygen ions were removed, at random, from the anion sublattice to give the compositions detailed in Table 1. Five different random configurations were created for each composition. For both the unstrained and strained calculations averages were obtained from all five random structures. Before calculating the diffusion data, each structure was first equilibrated for 16 ps: during the first 8 ps molecular dynamics was performed in the NVT (constant volume and constant temperature) thermodynamic ensemble [33]. In order to ensure rapid convergence of temperature with respect to time, the Berendsen thermostat was used to control the temperature of the system during this period [41]. A further 158 ps of dynamics was performed in the NPT (constant pressure, constant temperature) ensemble in which both lattice volume and atom positions are allowed to relax with time (temperature and pressure were controlled using the Nosé–Hoover thermostat and barostat respectively) [42–44]. Ionic diffusion can be examined by considering the evolution of the mean squared displacement (MSD) of atoms in a molecular dynamics simulation as a function of time. In order to obtain statistically significant values squared displacements are averaged over a large number of ions to give mean squared displacements (in the present case averages were made over all the oxygen ions in the simulated systems):
10−9 6
39
r
E D 2 MSDðt Þ ¼ r i ðt Þ−r i ð0Þ
ð4Þ
10000/T (K−1) where ri(t) is the position of ion i at time t. The diffusion coefficient, D, can be related to the MSD through:
T (K)
10
1273
1473
1873 −5
E D 2 r i ðt Þ−r i ð0Þ ∝6Dt:
1073
c) Gd:Y = 0:1
D (cm2s−1)
10−6
10−7
10−8
10−9 6
7
8
9
10
10000/T (K−1) Strain −0.025
−0.010
0.000
0.010
0.025
Fig. 3. Arrhenius plots showing the effect of temperature and strain on diffusion coefficient for CeO2 doped with 9 mol% X2O3 when a) X = Gd3+ (1:0), b) X = Gd3+,Y3+ with the co-doping ratio of 1:1 and c) X = Y3+ (0:1). For clarity, the ±0.020 and ±0.015 strain results are not shown here.
with Y). For the strained systems the 2:1, 1:1 and 0:1 systems were considered. Doping CeO2 with trivalent cations leads to an oxygen deficient structure where the trivalent cations substitute at Ce lattice sites.
ð5Þ
Hence, for each combination of composition and temperature, oxygen ion diffusion coefficients were obtained as the gradient of a linear least squares fit to a plot of MSD against time. A further 150 ps of NVT dynamics was performed during which MSDs were extracted. Due to the slower lattice dynamics experienced at the lowest temperature (973 K), it was found necessary to equilibrate each configuration for a further 150 ps (giving a total of 300 ps NVT MD), in order to fully equilibrate the structures before extracting the diffusion data, as for the higher temperatures, MSDs were extracted from the final 150 ps of this run. MSDs were extracted from the final 150 ps of a constant pressure (NPT) dynamics with atomic coordinates being sampled at 0.05 ps intervals. In addition to the unstrained lattice diffusion calculations described above, a series of diffusion calculations were performed in which strain was induced in the material. Both compressive (negative) and tensile (positive) homogeneous triaxial strains (which are not achievable in practice) were considered in the range of −0.025 to 0.025 for each 0.005 strain increment. Strained lattices were generated for 1:0, 1:1 and 0:1 Gd:Y ratios using the following procedure. The relaxed system volume, obtained from the constant pressure equilibration (NPT) procedure (see above), was used to establish the simulation cell's equilibrium dimensions, at zero strain, for each combination of temperature and composition. Using this as a starting point, the system was then strained iteratively, by increasing (or decreasing) the x, y and z dimensions of the cell in 0.005 increments (at each step atom positions were also scaled appropriately). In order to equilibrate the structure at each strain
40
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
a) 1073 K
b) 1273 K D (1 × 107cm2s−1)
D (1 × 107cm2s−1)
Gd:Y 15
1:0 1:1 0:1
10
5
0
40 30 20 10 0
−2.5
−1.0
0.0
1.0
−2.5
2.0
−1.0
Strain (%)
0.0
1.0
2.0
1.0
2.0
Strain (%)
c) 1473 K
d) 1673 K 150
D (1 × 107cm2s−1)
D (1 × 107cm2s−1)
80 60 40 20 0
100
50
0 −2.5
0.0
−1.0
1.0
2.0
−2.5
−1.0
Strain (%)
0.0
Strain (%)
Fig. 4. The effect of strain on the diffusion coefficient of CeO2 doped with 9 mol% X2O3 when X = Gd3+ (1:0), X = Y3+ (0:1) and when X = Gd3+,Y3+ with the co-doping ratio of 1:1 for temperatures of (a) 1073 K, (b) 1273 K, (c) 1473 K, and (d) 1673 K.
level, 5 ps of constant volume MD was performed (using a Nosé– Hoover thermostat to control temperature). Subsequently, oxygen MSDs were obtained from each equilibrated, strained lattice by performing 150 ps of constant volume dynamics; again atomic coordinates were sampled at 0.05 ps intervals.
for all the Gd:Y ratios considered is almost isotropic (MSD for all directions are equal). As expected for an oxygen deficient fluorite structure oxygen transport is mediated by the oxygen vacancies, which charge balance the trivalent dopants in the CeO2 lattice. 3.2. Impact of strain
3. Results and discussion 3.1. Yttria and gadolinia co-doped ceria Fig. 1 shows the calculated oxygen diffusion coefficient, D, versus the reciprocal temperature (T) for 9 mol% X2O3 in CeO2 where X = Gd 3+ and/or Y 3+ for the range of Gd:Y ratios previously specified (i.e. 1:0, 2:1, 1:1, 1:2 and 0:1), across the 973–1873 K temperature range. The diffusivities are similar irrespective of the Gd:Y ratios. In particular, there is only a very small range in the activation energies of diffusion, Ea, from 0.70 eV (9 mol% Gd2O3 in CeO2) to 0.75 eV (9 mol% Y2O3 in CeO2). The diffusivity results are plotted with the available experimental data for yttria and gadolinia single doped compositions in Fig. 2 and are seen to be in good general agreement [45–49]. The oxygen diffusion
Table 2 Activation energies (in eV) for oxygen migration obtained from the gradient of plots given in Fig. 3. Gd:Y ratio
1:0 1:1 0:1
Strain −0.025
0.00
0.025
0.731 0.743 0.711
0.665 0.710 0.706
0.524 0.545 0.527
Fig. 3 shows an Arrhenius plot of the calculated oxygen diffusion coefficient for 9 mol% (a) Gd2O3 doped CeO2, (b) (Gd,Y)2O3 co-doped CeO2 and (c) Y2O3 doped CeO2 for up to −2.5% compressive and 2.5% tensile homogeneous triaxial strains. It is found that tensile strain increases the oxygen diffusion coefficient considerably, whereas compressive strain has the opposite effect. For example, at 1273 K a 2.5% tensile strain leads to an almost five fold increase in diffusivity whereas a 2.5% compressive strain leads to a similar decrease in diffusivity (see Fig. 4). At temperatures lower than 1273 K the changes in diffusivity due to strain are slightly smaller, around a factor of three, while at higher temperatures the changes remain around a factor of five or perhaps slightly less. Activation energies were calculated by taking the gradient of Fig. 3 and are given in Table 2. To gain further insight into the mechanism responsible for the observed change in diffusivity with strain and co-doping, we have investigated the formation of defect clusters in doped CeO2. Specifically, we have used static atomistic simulations to calculate the binding energies (refer to Eq. 2) of defect clusters for a range of strain conditions. We con/ / sidered the positively charged pairs, {YCe :VO••}• and {GdCe : VO••}•, and also / •• / × / •• / × / the charge neutral clusters: {YCe :VO :YCe} , {YCe : VO :GdCe} and {GdCe : •• / × VO : GdCe} . We calculated the binding energies for different possible configurations for all the different sites. Consequently, there are numerous cluster configurations for every strain condition and composition but here we only report the minimum energy cluster in each case (see
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
41
Gd
Gd Gd NN 2NN 3NN
Binding Energy (eV)
1.2
1.0
0.8
Binding Energy (eV)
1.4 2.2
NN−NN NN−2NN 2NN−2NN
2.0 1.8 1.6 1.4 1.2 1.0
0.6 −2
−1
0
1
2
1
2
1
2
Strain (%) 0.4 −2
−1
0
1
Gd Y
2
Binding Energy (eV)
Strain (%)
Y 1.4
Binding Energy (eV)
1.2
1.0
2.2 2.0 1.8 1.6 1.4 1.2 1.0
0.8
−2
−1
0
Strain (%) 0.6
YY
−2
−1
0
1
2
Strain (%) / •• • / •• • Fig. 5. Binding energies of the {GdCe : VO } (top) and {YCe : VO } (bottom) at NN, 2NN and third nearest neighbour configurations (3NN).
/ / Figs. 5 and 6). Fig. 5 reports the binding energies of {YCe :VO••}• and {GdCe : VO••}• at the nearest neighbour (NN), second nearest neighbour (2NN) / and third nearest neighbour configurations (3NN). For the {YCe :VO••}• cluster (or dimer) it was found that the NN pairs are the strongest / bound. For the {GdCe : VO••}• cluster, the NN configuration is again most / stable but only up to 1.5% tensile strain above which the 2NN, {GdCe : VO••}•, becomes favourable. As expected the binding energies are higher under compression as there is a reduction in the defect volume of a cluster compared to its components (see discussion in Ref. [27]). / / × / / × Fig. 6 reports the binding energies of {YCe :VO•• :YCe } , {YCe :VO•• :GdCe } / •• / × and {GdCe :VO :GdCe} clusters. These energies are significantly greater than the pairs and therefore will provide a significant driving force for the oxygen vacancies to bind with the trivalent dopants. Interestingly, / / × for {GdCe :VO•• :GdCe } the NN–NN configuration is most energetically favoured (where both Gd ions are first neighbour with respect to the oxygen vacancy) under negative strains, whereas the 2NN–2NN configuration is preferred (where the Gd atoms are both at 2NN positions with / respect to the oxygen vacancy) under a positive strain. For the {YCe : / × / / × VO•• :YCe } and {YCe :VO•• :GdCe } clusters there is an analogous crossover from the NN–NN to the 2NN–2NN but at higher strains (1% and 2% respectively). Despite the greater stability of the neutral trimer clusters over equivalent dimer clusters, the pairs may be more populous at higher temperatures due to their lower configurational entropy. Irrespective of the relative concentrations of dimers and trimers, both Figs. 5 and 6 indicate that compressive strain significantly increases the association of oxygen vacancies for all the defect clusters
Binding Energy (eV)
0.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 −2
−1
0
Strain (%) / •• / × / •• / × / •• / × Fig. 6. Binding energies of (a) {YCe :VO :YCe } , (b) {YCe :VO :GdCe } and (c) {GdCe :VO :GdCe } clusters.
considered: meaning that compressing the ceria lattice increases the stability of defect clusters containing oxygen vacancies. This association represents an impediment to oxygen vacancy diffusion and consequently, is expected to be an important effect in the decrease in diffusion coefficient predicted by MD for negative strains and shown in Fig. 3. 4. Conclusions The present study has considered the complexities of defect processes in co-doped relaxed and strained ceria using static and dynamic atomic scale computer simulation techniques. It has been shown here that in CeO2, over the whole co-doping range there is very little effect of co-doping on the oxygen diffusion. Conversely, strain has a significant effect upon the diffusivity with compressive strain leading to lower diffusivity while tensile strain leads to higher diffusivity. The static simulations indicated that tensile strain reduces the binding energies of clusters between oxygen vacancies and trivalent dopant
42
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
atoms while compressive strain leads to higher binding energies. This provides one explanation as to why tensile strain leads to the higher diffusivity seen in the molecular dynamics simulations. Acknowledgements Computing resources were provided by the HPC Facility of Imperial College London. References [1] N.P. Brandon, S. Skinner, B.C.H. Steele, Annu. Rev. Mater. Res. 33 (2003) 183. [2] J. Fleig, Annu. Rev. Mater. Res. 33 (2003) 36. [3] D. Rupasov, A. Chroneos, D. Parfitt, J.A. Kilner, R.W. Grimes, S.Y. Istomin, E.V. Antipov, Phys. Rev. B 79 (2009) 172102. [4] Z. Shao, S.M. Haile, J. Ahn, P.D. Ronney, Z. Zhan, S.A. Barnett, Nature 435 (2005) 795. [5] S.B. Adler, Solid State Ionics 111 (1998) 125. [6] A. Kushima, D. Parfitt, A. Chroneos, B. Yildiz, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 13 (2011) 2242. [7] A. Chroneos, D. Parfitt, J.A. Kilner, R.W. Grimes, J. Mater. Chem. 20 (2010) 266. [8] D. Parfitt, A. Chroneos, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 12 (2010) 6834. [9] N. Sata, K. Eberman, K. Eberl, J. Maier, Nature 408 (2000) 946. [10] T. Suzuki, I. Kosacki, H.U. Anderson, Solid State Ionics 151 (2002) 111. [11] I. Kosacki, C.M. Rouleau, P.F. Becher, J. Bentley, D.H. Lowndes, Solid State Ionics 176 (2005) 1319. [12] X.X. Guo, I. Matei, J.S. Lee, J. Maier, Appl. Phys. Lett. 92 (2007) 103102. [13] J.G. Barriocanal, A.R. Calzada, M. Varela, Z. Sefrioui, E. Iborra, C. Leon, S.J. Pennycook, J. Santamaria, Science 321 (2008) 676. [14] J.A. Kilner, Nat. Mater. 7 (2008) 838. [15] A. Kushima, B. Yildiz, J. Mater. Chem. 20 (2010) 4809. [16] X. Xie, R.V. Kumar, J. Sun, L.J. Henson, J. Power Sources 195 (2010) 5660. [17] B.C.H. Steele, A. Heinzel, Nature 411 (2001) 345. [18] M. Yashima, T. Takizawa, J. Phys. Chem. C 114 (2010) 2385. [19] H.L. Tuller, J. Phys. Chem. Solids 55 (1994) 1343. [20] A. Lashtabeg, S.J. Skinner, J. Mater. Chem. 16 (2006) 3161.
[21] P.P. Dholabhai, J.B. Adams, P. Crozier, R. Sharma, Phys. Chem. Chem. Phys. 12 (2010) 7904. [22] J.M. Ralph, J. Przydatek, J.A. Kilner, T. Seguelong, Ber. Bunsen Ges. Phys. Chem. 101 (1997) 1403. [23] S. Omar, E.D. Wachsman, J.C. Nino, Solid State Ionics 177 (2006) 3199. [24] F. Wang, A. Chroneos, C. Jiang, U. Schwingenschlögl, Phys. Chem. Chem. Phys. 14 (2012) 11737. [25] J.A. Kilner, Chem. Lett. 37 (2008) 1012. [26] I. Seymour, A. Chroneos, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 13 (2011) 15305. [27] L. Minervini, M.O. Zacate, R.W. Grimes, Solid State Ionics 116 (1999) 339. [28] G. Dezanneau, J. Hermet, B. Dupé, Int. J. Hydrogen Energ. 37 (2012) 8081. [29] R.A. De Souza, A. Ramadan, S. Hörner, Energy Environ. Sci. 5 (2012) 5445. [30] A. Chroneos, B. Yildiz, A. Tarancon, D. Parfitt, J.A. Kilner, Energy Environ. Sci. 4 (2011) 2774. [31] M. Born, J.E. Mayer, Z. Phys. 75 (1932) 1. [32] R.A. Buckingham, Proc. R. Soc. Lond. A: Math. Phys. Sci. 168 (1938) 264. [33] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1989. [34] P.P. Ewald, Ann. Phys. 64 (1921) 253. [35] W. Smith, I.T. Todorov, Mol. Simul. 32 (2006) 935. [36] N.F. Mott, M.J. Littleton, Trans. Faraday Soc. 34 (1938) 485. [37] J.D. Gale, J. Chem. Soc., Faraday Trans. 93 (1997) 629. [38] D.C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2004. [39] W. Smith, T.R. Forester, J. Mol. Graphics 14 (1996) 136. [40] F.A. Kröger, H.J. Vink, Solid State Phys. 3 (1956) 307. [41] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [42] S. Nosé, J. Chem. Phys. 81 (1984) 511. [43] W.G. Hoover, Phys. Rev. A 31 (1985) 1695. [44] S. Nosé, Mol. Phys. 52 (1984) 255. [45] H. Inaba, H. Tagawa, Solid State Ionics 83 (1996) 1. [46] J.M. Floyd, Indian J. Technol. 11 (1973) 589. [47] P.S. Manning, J.D. Sirman, J.A. Kilner, Solid State Ionics 93 (1997) 125. [48] E. Ruiz-Trejo, J.D. Sirman, Yu.M. Baikov, J.A. Kilner, Solid State Ionics 113 (1998) 565. [49] I.E.L. Stephens, J.A. Kilner, Solid State Ionics 177 (2006) 669.
Contents lists available at SciVerse ScienceDirect
Solid State Ionics journal homepage: www.elsevier.com/locate/ssi
Effect of strain on the oxygen diffusion in yttria and gadolinia co-doped ceria M.J.D. Rushton b, A. Chroneos a,⁎, S.J. Skinner b, J.A. Kilner b, R.W. Grimes b a b
Materials Engineering, The Open University, Milton Keynes, MK7 6AA, United Kingdom Department of Materials, Imperial College London, London SW7 2AZ, United Kingdom
a r t i c l e
i n f o
Article history: Received 23 May 2012 Received in revised form 25 August 2012 Accepted 12 September 2012 Available online 6 October 2012 Keywords: Doped ceria Molecular dynamics Strain
a b s t r a c t Atomic scale simulations have been used to investigate the impact of co-doping (yttrium and gadolinium) and strain on oxygen diffusion and binding of dopant–vacancy clusters in ceria. Doped ceria in its relaxed or strained form is used as an electrolyte for solid oxide fuel cell applications. For unstrained co-doped ceria we calculate an activation energy for migration of 0.70–0.75 eV in the temperature range of 973–1873 K. Co-doping with yttrium and gadolinium only affected oxide ion diffusion to a small degree when compared to single doping. The diffusion coefficient was substantially increased by tensile strain while compressive strain caused a decrease. To gain further insight why tensile strain leads to higher diffusivity static simulations were employed. It is calculated that tensile strain reduces the binding energies of clusters between oxygen vacancies and trivalent dopant atoms while compressive strain leads to higher binding energies. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Solid oxide fuel cells (SOFCs) face two major obstacles to a large scale application. These are the high cost of current SOFC designs and their degradation over the long timescales needed for commercial application. One route for SOFCs to overcome these barriers is through operation within the intermediate temperature range (500–700 °C) [1]. This presents difficulties for traditional ceramic electrodes and electrolytes as it requires materials with enhanced activity, which leads to a requirement for higher oxygen diffusivities [2]. One way of achieving this is by identifying new materials and a number of candidate oxides are being investigated for cathode and electrolyte applications in the next generation intermediate temperature SOFCs (IT-SOFCs) [3–8]. An alternative way of enhancing ionic transport may be through the presence of an interface between dissimilar oxides, since the effects of interfacial strain are thought to affect oxide ion transport [9–15]. More conventional methods for enhancing conductivity involve co-doping, the idea being to minimise the effect of any defect interactions through the introduction of heterogeneous dopant species [16]. Trivalent-doped ceria is a well established electrolyte material due to its high oxygen ion diffusivity and relatively lower reduction temperatures [17]. Doping with oxides such as R2O3 (here R = Gd and/or Y) improves oxygen diffusion due to the formation of oxygen vacancies in the resultant Ce1 −xRxO2−x/2 solid solution [18]. Doped ceria has received considerable attention due to its application as the electrolyte material in SOFCs as it exhibits up to 2–3 orders of magnitude higher oxygen ion conductivity than the common high temperature fuel cell
electrolyte, and yttria stabilised zirconia, at intermediate temperatures [19–21]. Co-doping may offer an attractive method to increase the conductivity in CeO2 even further [22–25]. Atomic scale simulations have been carried out in order to gain insight into the effects of strain and/or doping regimes on oxygen transport in SOFC materials [26–30]. Recently De Souza et al. considered the effect of strain on oxygen diffusivity within ceria using a static lattice method and predicted a considerable reduction in the activation energy for oxygen vacancy migration on dilation of the lattice [29]. Here, we additionally consider the combined effects of single and co-doping on strained and relaxed lattices. In a similar manner to De Souza's study, static Mott–Littleton calculations are used to consider the interaction of point defects within strained ceria, furthermore we use molecular dynamics (MD) calculations to observe how strain affects oxygen ion diffusivity as a function of temperature. 2. Methodology Both the MD and the static Mott–Littleton simulation methodologies employed in this study rely on a simple description of interatomic forces. Here the Born description of ionic crystals was used [31]: forces acting between atoms are resolved into long-range coulombic and shortrange pair-wise interactions described using parameterized Buckingham pair potentials [32]. The lattice energy, EL, is given by,
EL ¼
N −1 X N X i¼1 j>i
⁎ Corresponding author. E-mail addresses: [email protected] (M.J.D. Rushton), [email protected] (A. Chroneos). 0167-2738/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ssi.2012.09.015
! qi qj r ij C ij þ Aij exp − − 6 4πε0 r ij ρij r ij
ð1Þ
where N is the number of ions in the system, ε0 is the permittivity of free space, rij is the separation between ions i and j of charges qi and qj
38
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
Table 1 Atomic composition of the simulated systems with the chemical formula Ce0.835X0.165O1.96. X
Number of atoms
Gd:Y
Gd
Y
Ce
O
1:0 2:1 1:1 1:2 0:1
338 225 169 113 0
0 113 169 225 338
1710 1710 1710 1710 1710
3927 3927 3927 3927 3927
dielectric continuum approximation, with the ionic relaxations taking place in a single step. The interaction energies between ions in Regions I and IIa were calculated explicitly. Finally, the outer Region IIb consists of a point charge array providing the Madelung field of the remaining lattice, and the response energy in Region IIb was calculated using the Mott–Littleton approximation [36]. The GULP code [37] was used for all the static atomistic simulations. The interaction of point defects within the lattice can be quantified by the calculation of their binding energies, Eb, defined by Eb ¼ ∑Eisolateddefects −Edefectcluster :
respectively, and Aij, ρij and Cij are the short-range Buckingham potential parameters specific to particular pairs of interacting species. The potential parameters have already been employed successfully in static simulations of doped ceria [27]. Here they are used without the shell model of the previous work to facilitate computationally tractable MD simulations. A short-range potential cut-off of 6.5 Å was used. As full periodic boundaries [33] were employed, long range Coulombic forces were calculated using the Ewald sum [34,35].
ð2Þ
According to this definition a positive binding energy implies that the defect cluster is more energetically favourable as compared to the isolated defects. 2.2. Molecular dynamics
Static calculations were used to determine the binding energies of point defects within the ceria lattice. The response of the crystal to charged defects was determined using the Mott–Littleton approach. Within this, a perfect crystal lattice is initially constructed: periodic boundaries were employed whereby a unit-cell is tessellated in space to represent a continuous crystal. The energy of the system is minimised by relaxing atomic positions and cell parameters via a conjugate gradient approach. The relaxed structure is then used to create concentric spherical regions. Point defects are introduced at the centre of the innermost region (Region I). Here the ions are treated explicitly and are relaxed subject to forces calculated from the interionic potentials. In Region IIa, which surrounds Region I, the ionic forces are calculated via a
MD simulations were used to investigate the effects on oxygen ion diffusion of co-doping ceria with gadolinia and yttria, as a function of temperature (between 973 K and 1873 K). In addition strain in the range of −2.5 to 2.5% was applied to doped compositions in order to examine its influence on diffusion. During MD simulations, Newton's second law of motion is solved iteratively. This means that the positions and velocities of a set of atoms can evolve as a function of time from a description of interatomic forces [33,38]. Here the velocity Verlet integration method was used [35,38] and an integration time step of 1 fs was used throughout. Molecular dynamics calculations were performed using the DL_POLY code (version 3.09) [35,39]. Each simulated system contained CeO2 doped with 9 mol% of X2O3 (where X could be Gd3+ and/or Y 3+). For the constant pressure (i.e. zero strain) simulations, five Gd:Y ratios were considered ranging from doping only with Gd (1:0), through co-doped compositions with ratios of 2:1, 1:1, 1:2, to 0:1 (the last corresponding to doping only
Fig. 1. Arrhenius plot of the calculated oxygen diffusion coefficient (D) against temperature (T) for 9 mol% X2O3 in CeO2 for X = Gd3+ and/or Y3+ with ratios as indicated in the figure. Error bars indicate the highest and lowest D values for the five random configurations considered while points show the mean. For clarity, at each temperature points have been spread along the x-axis with grey drop lines representing the temperature each group actually represents.
Fig. 2. Oxygen diffusion coefficient (D) against temperature (T) for 9 mol% X2O3 in CeO2 for X = Gd3+ and/or Y3+ with ratios as indicated in the figure. Compared with experimentally determined diffusion coefficients for singly doped CeO2 compositions. For the purposes of comparison, the 9 mol% composition studied here, has the chemical formula: Ce0.835X0.165O1.918.
2.1. Static atomistic simulations
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
T (K) 10−5
1273
1473
1873
The defect equation (in Kröger Vink notation) [40] for this process is as follows:
1073
a) Gd:Y = 1:0
CeO2
X 2 O3 → 2X
D (cm2s−1)
10−6
10−7
10−8
10−9 6
7
8
9
10
10000/T (K−1) T (K) 10−5
1273
1473
1873
1073
b) Gd:Y = 1:1
D (cm2s−1)
10−6
10−7
10−8
7
8
9
10
′
••
Ce
X
þ VO þ 3OO :
ð3Þ
Every two trivalent ions introduced into the cation sublattice lead to the formation of a charge compensating oxygen vacancy. With this in mind, each simulation cell was obtained from an 8 × 8 × 8 CeO2 super-cell (originally containing 6144 atoms). To this 338 dopant ions were randomly assigned to the Ce sites in the desired Gd:Y ratio while 169 oxygen ions were removed, at random, from the anion sublattice to give the compositions detailed in Table 1. Five different random configurations were created for each composition. For both the unstrained and strained calculations averages were obtained from all five random structures. Before calculating the diffusion data, each structure was first equilibrated for 16 ps: during the first 8 ps molecular dynamics was performed in the NVT (constant volume and constant temperature) thermodynamic ensemble [33]. In order to ensure rapid convergence of temperature with respect to time, the Berendsen thermostat was used to control the temperature of the system during this period [41]. A further 158 ps of dynamics was performed in the NPT (constant pressure, constant temperature) ensemble in which both lattice volume and atom positions are allowed to relax with time (temperature and pressure were controlled using the Nosé–Hoover thermostat and barostat respectively) [42–44]. Ionic diffusion can be examined by considering the evolution of the mean squared displacement (MSD) of atoms in a molecular dynamics simulation as a function of time. In order to obtain statistically significant values squared displacements are averaged over a large number of ions to give mean squared displacements (in the present case averages were made over all the oxygen ions in the simulated systems):
10−9 6
39
r
E D 2 MSDðt Þ ¼ r i ðt Þ−r i ð0Þ
ð4Þ
10000/T (K−1) where ri(t) is the position of ion i at time t. The diffusion coefficient, D, can be related to the MSD through:
T (K)
10
1273
1473
1873 −5
E D 2 r i ðt Þ−r i ð0Þ ∝6Dt:
1073
c) Gd:Y = 0:1
D (cm2s−1)
10−6
10−7
10−8
10−9 6
7
8
9
10
10000/T (K−1) Strain −0.025
−0.010
0.000
0.010
0.025
Fig. 3. Arrhenius plots showing the effect of temperature and strain on diffusion coefficient for CeO2 doped with 9 mol% X2O3 when a) X = Gd3+ (1:0), b) X = Gd3+,Y3+ with the co-doping ratio of 1:1 and c) X = Y3+ (0:1). For clarity, the ±0.020 and ±0.015 strain results are not shown here.
with Y). For the strained systems the 2:1, 1:1 and 0:1 systems were considered. Doping CeO2 with trivalent cations leads to an oxygen deficient structure where the trivalent cations substitute at Ce lattice sites.
ð5Þ
Hence, for each combination of composition and temperature, oxygen ion diffusion coefficients were obtained as the gradient of a linear least squares fit to a plot of MSD against time. A further 150 ps of NVT dynamics was performed during which MSDs were extracted. Due to the slower lattice dynamics experienced at the lowest temperature (973 K), it was found necessary to equilibrate each configuration for a further 150 ps (giving a total of 300 ps NVT MD), in order to fully equilibrate the structures before extracting the diffusion data, as for the higher temperatures, MSDs were extracted from the final 150 ps of this run. MSDs were extracted from the final 150 ps of a constant pressure (NPT) dynamics with atomic coordinates being sampled at 0.05 ps intervals. In addition to the unstrained lattice diffusion calculations described above, a series of diffusion calculations were performed in which strain was induced in the material. Both compressive (negative) and tensile (positive) homogeneous triaxial strains (which are not achievable in practice) were considered in the range of −0.025 to 0.025 for each 0.005 strain increment. Strained lattices were generated for 1:0, 1:1 and 0:1 Gd:Y ratios using the following procedure. The relaxed system volume, obtained from the constant pressure equilibration (NPT) procedure (see above), was used to establish the simulation cell's equilibrium dimensions, at zero strain, for each combination of temperature and composition. Using this as a starting point, the system was then strained iteratively, by increasing (or decreasing) the x, y and z dimensions of the cell in 0.005 increments (at each step atom positions were also scaled appropriately). In order to equilibrate the structure at each strain
40
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
a) 1073 K
b) 1273 K D (1 × 107cm2s−1)
D (1 × 107cm2s−1)
Gd:Y 15
1:0 1:1 0:1
10
5
0
40 30 20 10 0
−2.5
−1.0
0.0
1.0
−2.5
2.0
−1.0
Strain (%)
0.0
1.0
2.0
1.0
2.0
Strain (%)
c) 1473 K
d) 1673 K 150
D (1 × 107cm2s−1)
D (1 × 107cm2s−1)
80 60 40 20 0
100
50
0 −2.5
0.0
−1.0
1.0
2.0
−2.5
−1.0
Strain (%)
0.0
Strain (%)
Fig. 4. The effect of strain on the diffusion coefficient of CeO2 doped with 9 mol% X2O3 when X = Gd3+ (1:0), X = Y3+ (0:1) and when X = Gd3+,Y3+ with the co-doping ratio of 1:1 for temperatures of (a) 1073 K, (b) 1273 K, (c) 1473 K, and (d) 1673 K.
level, 5 ps of constant volume MD was performed (using a Nosé– Hoover thermostat to control temperature). Subsequently, oxygen MSDs were obtained from each equilibrated, strained lattice by performing 150 ps of constant volume dynamics; again atomic coordinates were sampled at 0.05 ps intervals.
for all the Gd:Y ratios considered is almost isotropic (MSD for all directions are equal). As expected for an oxygen deficient fluorite structure oxygen transport is mediated by the oxygen vacancies, which charge balance the trivalent dopants in the CeO2 lattice. 3.2. Impact of strain
3. Results and discussion 3.1. Yttria and gadolinia co-doped ceria Fig. 1 shows the calculated oxygen diffusion coefficient, D, versus the reciprocal temperature (T) for 9 mol% X2O3 in CeO2 where X = Gd 3+ and/or Y 3+ for the range of Gd:Y ratios previously specified (i.e. 1:0, 2:1, 1:1, 1:2 and 0:1), across the 973–1873 K temperature range. The diffusivities are similar irrespective of the Gd:Y ratios. In particular, there is only a very small range in the activation energies of diffusion, Ea, from 0.70 eV (9 mol% Gd2O3 in CeO2) to 0.75 eV (9 mol% Y2O3 in CeO2). The diffusivity results are plotted with the available experimental data for yttria and gadolinia single doped compositions in Fig. 2 and are seen to be in good general agreement [45–49]. The oxygen diffusion
Table 2 Activation energies (in eV) for oxygen migration obtained from the gradient of plots given in Fig. 3. Gd:Y ratio
1:0 1:1 0:1
Strain −0.025
0.00
0.025
0.731 0.743 0.711
0.665 0.710 0.706
0.524 0.545 0.527
Fig. 3 shows an Arrhenius plot of the calculated oxygen diffusion coefficient for 9 mol% (a) Gd2O3 doped CeO2, (b) (Gd,Y)2O3 co-doped CeO2 and (c) Y2O3 doped CeO2 for up to −2.5% compressive and 2.5% tensile homogeneous triaxial strains. It is found that tensile strain increases the oxygen diffusion coefficient considerably, whereas compressive strain has the opposite effect. For example, at 1273 K a 2.5% tensile strain leads to an almost five fold increase in diffusivity whereas a 2.5% compressive strain leads to a similar decrease in diffusivity (see Fig. 4). At temperatures lower than 1273 K the changes in diffusivity due to strain are slightly smaller, around a factor of three, while at higher temperatures the changes remain around a factor of five or perhaps slightly less. Activation energies were calculated by taking the gradient of Fig. 3 and are given in Table 2. To gain further insight into the mechanism responsible for the observed change in diffusivity with strain and co-doping, we have investigated the formation of defect clusters in doped CeO2. Specifically, we have used static atomistic simulations to calculate the binding energies (refer to Eq. 2) of defect clusters for a range of strain conditions. We con/ / sidered the positively charged pairs, {YCe :VO••}• and {GdCe : VO••}•, and also / •• / × / •• / × / the charge neutral clusters: {YCe :VO :YCe} , {YCe : VO :GdCe} and {GdCe : •• / × VO : GdCe} . We calculated the binding energies for different possible configurations for all the different sites. Consequently, there are numerous cluster configurations for every strain condition and composition but here we only report the minimum energy cluster in each case (see
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
41
Gd
Gd Gd NN 2NN 3NN
Binding Energy (eV)
1.2
1.0
0.8
Binding Energy (eV)
1.4 2.2
NN−NN NN−2NN 2NN−2NN
2.0 1.8 1.6 1.4 1.2 1.0
0.6 −2
−1
0
1
2
1
2
1
2
Strain (%) 0.4 −2
−1
0
1
Gd Y
2
Binding Energy (eV)
Strain (%)
Y 1.4
Binding Energy (eV)
1.2
1.0
2.2 2.0 1.8 1.6 1.4 1.2 1.0
0.8
−2
−1
0
Strain (%) 0.6
YY
−2
−1
0
1
2
Strain (%) / •• • / •• • Fig. 5. Binding energies of the {GdCe : VO } (top) and {YCe : VO } (bottom) at NN, 2NN and third nearest neighbour configurations (3NN).
/ / Figs. 5 and 6). Fig. 5 reports the binding energies of {YCe :VO••}• and {GdCe : VO••}• at the nearest neighbour (NN), second nearest neighbour (2NN) / and third nearest neighbour configurations (3NN). For the {YCe :VO••}• cluster (or dimer) it was found that the NN pairs are the strongest / bound. For the {GdCe : VO••}• cluster, the NN configuration is again most / stable but only up to 1.5% tensile strain above which the 2NN, {GdCe : VO••}•, becomes favourable. As expected the binding energies are higher under compression as there is a reduction in the defect volume of a cluster compared to its components (see discussion in Ref. [27]). / / × / / × Fig. 6 reports the binding energies of {YCe :VO•• :YCe } , {YCe :VO•• :GdCe } / •• / × and {GdCe :VO :GdCe} clusters. These energies are significantly greater than the pairs and therefore will provide a significant driving force for the oxygen vacancies to bind with the trivalent dopants. Interestingly, / / × for {GdCe :VO•• :GdCe } the NN–NN configuration is most energetically favoured (where both Gd ions are first neighbour with respect to the oxygen vacancy) under negative strains, whereas the 2NN–2NN configuration is preferred (where the Gd atoms are both at 2NN positions with / respect to the oxygen vacancy) under a positive strain. For the {YCe : / × / / × VO•• :YCe } and {YCe :VO•• :GdCe } clusters there is an analogous crossover from the NN–NN to the 2NN–2NN but at higher strains (1% and 2% respectively). Despite the greater stability of the neutral trimer clusters over equivalent dimer clusters, the pairs may be more populous at higher temperatures due to their lower configurational entropy. Irrespective of the relative concentrations of dimers and trimers, both Figs. 5 and 6 indicate that compressive strain significantly increases the association of oxygen vacancies for all the defect clusters
Binding Energy (eV)
0.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 −2
−1
0
Strain (%) / •• / × / •• / × / •• / × Fig. 6. Binding energies of (a) {YCe :VO :YCe } , (b) {YCe :VO :GdCe } and (c) {GdCe :VO :GdCe } clusters.
considered: meaning that compressing the ceria lattice increases the stability of defect clusters containing oxygen vacancies. This association represents an impediment to oxygen vacancy diffusion and consequently, is expected to be an important effect in the decrease in diffusion coefficient predicted by MD for negative strains and shown in Fig. 3. 4. Conclusions The present study has considered the complexities of defect processes in co-doped relaxed and strained ceria using static and dynamic atomic scale computer simulation techniques. It has been shown here that in CeO2, over the whole co-doping range there is very little effect of co-doping on the oxygen diffusion. Conversely, strain has a significant effect upon the diffusivity with compressive strain leading to lower diffusivity while tensile strain leads to higher diffusivity. The static simulations indicated that tensile strain reduces the binding energies of clusters between oxygen vacancies and trivalent dopant
42
M.J.D. Rushton et al. / Solid State Ionics 230 (2013) 37–42
atoms while compressive strain leads to higher binding energies. This provides one explanation as to why tensile strain leads to the higher diffusivity seen in the molecular dynamics simulations. Acknowledgements Computing resources were provided by the HPC Facility of Imperial College London. References [1] N.P. Brandon, S. Skinner, B.C.H. Steele, Annu. Rev. Mater. Res. 33 (2003) 183. [2] J. Fleig, Annu. Rev. Mater. Res. 33 (2003) 36. [3] D. Rupasov, A. Chroneos, D. Parfitt, J.A. Kilner, R.W. Grimes, S.Y. Istomin, E.V. Antipov, Phys. Rev. B 79 (2009) 172102. [4] Z. Shao, S.M. Haile, J. Ahn, P.D. Ronney, Z. Zhan, S.A. Barnett, Nature 435 (2005) 795. [5] S.B. Adler, Solid State Ionics 111 (1998) 125. [6] A. Kushima, D. Parfitt, A. Chroneos, B. Yildiz, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 13 (2011) 2242. [7] A. Chroneos, D. Parfitt, J.A. Kilner, R.W. Grimes, J. Mater. Chem. 20 (2010) 266. [8] D. Parfitt, A. Chroneos, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 12 (2010) 6834. [9] N. Sata, K. Eberman, K. Eberl, J. Maier, Nature 408 (2000) 946. [10] T. Suzuki, I. Kosacki, H.U. Anderson, Solid State Ionics 151 (2002) 111. [11] I. Kosacki, C.M. Rouleau, P.F. Becher, J. Bentley, D.H. Lowndes, Solid State Ionics 176 (2005) 1319. [12] X.X. Guo, I. Matei, J.S. Lee, J. Maier, Appl. Phys. Lett. 92 (2007) 103102. [13] J.G. Barriocanal, A.R. Calzada, M. Varela, Z. Sefrioui, E. Iborra, C. Leon, S.J. Pennycook, J. Santamaria, Science 321 (2008) 676. [14] J.A. Kilner, Nat. Mater. 7 (2008) 838. [15] A. Kushima, B. Yildiz, J. Mater. Chem. 20 (2010) 4809. [16] X. Xie, R.V. Kumar, J. Sun, L.J. Henson, J. Power Sources 195 (2010) 5660. [17] B.C.H. Steele, A. Heinzel, Nature 411 (2001) 345. [18] M. Yashima, T. Takizawa, J. Phys. Chem. C 114 (2010) 2385. [19] H.L. Tuller, J. Phys. Chem. Solids 55 (1994) 1343. [20] A. Lashtabeg, S.J. Skinner, J. Mater. Chem. 16 (2006) 3161.
[21] P.P. Dholabhai, J.B. Adams, P. Crozier, R. Sharma, Phys. Chem. Chem. Phys. 12 (2010) 7904. [22] J.M. Ralph, J. Przydatek, J.A. Kilner, T. Seguelong, Ber. Bunsen Ges. Phys. Chem. 101 (1997) 1403. [23] S. Omar, E.D. Wachsman, J.C. Nino, Solid State Ionics 177 (2006) 3199. [24] F. Wang, A. Chroneos, C. Jiang, U. Schwingenschlögl, Phys. Chem. Chem. Phys. 14 (2012) 11737. [25] J.A. Kilner, Chem. Lett. 37 (2008) 1012. [26] I. Seymour, A. Chroneos, J.A. Kilner, R.W. Grimes, Phys. Chem. Chem. Phys. 13 (2011) 15305. [27] L. Minervini, M.O. Zacate, R.W. Grimes, Solid State Ionics 116 (1999) 339. [28] G. Dezanneau, J. Hermet, B. Dupé, Int. J. Hydrogen Energ. 37 (2012) 8081. [29] R.A. De Souza, A. Ramadan, S. Hörner, Energy Environ. Sci. 5 (2012) 5445. [30] A. Chroneos, B. Yildiz, A. Tarancon, D. Parfitt, J.A. Kilner, Energy Environ. Sci. 4 (2011) 2774. [31] M. Born, J.E. Mayer, Z. Phys. 75 (1932) 1. [32] R.A. Buckingham, Proc. R. Soc. Lond. A: Math. Phys. Sci. 168 (1938) 264. [33] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1989. [34] P.P. Ewald, Ann. Phys. 64 (1921) 253. [35] W. Smith, I.T. Todorov, Mol. Simul. 32 (2006) 935. [36] N.F. Mott, M.J. Littleton, Trans. Faraday Soc. 34 (1938) 485. [37] J.D. Gale, J. Chem. Soc., Faraday Trans. 93 (1997) 629. [38] D.C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2004. [39] W. Smith, T.R. Forester, J. Mol. Graphics 14 (1996) 136. [40] F.A. Kröger, H.J. Vink, Solid State Phys. 3 (1956) 307. [41] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [42] S. Nosé, J. Chem. Phys. 81 (1984) 511. [43] W.G. Hoover, Phys. Rev. A 31 (1985) 1695. [44] S. Nosé, Mol. Phys. 52 (1984) 255. [45] H. Inaba, H. Tagawa, Solid State Ionics 83 (1996) 1. [46] J.M. Floyd, Indian J. Technol. 11 (1973) 589. [47] P.S. Manning, J.D. Sirman, J.A. Kilner, Solid State Ionics 93 (1997) 125. [48] E. Ruiz-Trejo, J.D. Sirman, Yu.M. Baikov, J.A. Kilner, Solid State Ionics 113 (1998) 565. [49] I.E.L. Stephens, J.A. Kilner, Solid State Ionics 177 (2006) 669.