Modeling Oxygen Ion Migration in the CeO2-ZrO2-Y2O3 Solid Solution

Subscriber access provided by Kaohsiung Medical University

C: Energy Conversion and Storage; Energy and Charge Transport 2

2

2

3

Modeling Oxygen Ion Migration in the CeO-ZrO-YO Solid Solution Steffen Grieshammer, Sebastian Eisele, and Julius Koettgen J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b04361 • Publication Date (Web): 27 Jul 2018 Downloaded from http://pubs.acs.org on July 30, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Modeling Oxygen Ion Migration in the CeO2-ZrO2-Y2O3 Solid Solution Steffen Grieshammer, a,b,c* Sebastian Eisele, a,b Julius Koettgenb,c a

Helmholtz-Institut Münster (IEK-12), Forschungszentrum Jülich GmbH, Corrensstraße 46, 48149 Münster, Germany b

c

Institute of Physical Chemistry, RWTH Aachen University, Landoltweg 2, 52056 Aachen, Germany

JARA-HPC, RWTH Aachen University and Forschungszentrum Jülich GmbH, Germany

*Corresponding author: [email protected]

Abstract Zirconium doping in cerium oxide, as a result of intentional material engineering or unintentional impurities, impacts material properties like the reduction behavior and defect migration. In this study, we investigate the influence of zirconium doping on the conductivity of yttrium doped ceria using DFT+U calculations. We calculate the migration energies of oxygen ions for different jump environments containing yttrium and zirconium ions and compare the results to a simplified migration energy model. The small zirconium ions lead to strong distortions of the lattice, which results in deviation between calculated and modeled energies. Both the calculated and the modeled migration energies are used in Kinetic Monte Carlo simulations to obtain the ionic conductivity for various dopant fractions. We identify three major influences on the ionic conductivity: the trapping of oxygen vacancies by dopant ions, the blocking effect, which alters the migration barriers around defects and the lattice contraction due to zirconium doping.

I Introduction Doped ceria is of scientific interest due to its application in automotive three-way catalysts for the purification of exhaust gas, in industrial catalysis and its promising properties for the production of solar fuels.1 In addition, it exhibits high oxygen ion conductivity at intermediate temperatures making it a suitable electrolyte for solid oxide fuel cells, electrolyzer cells and high temperature batteries.2 The conductivity is enabled by doping with rare-earth oxides like Y2O3 leading to the formation of mobile oxygen vacancies (VO∙∙ in Kröger-Vink-notation):

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

CeO2

Y2 O3 →

′ 2YCe +VO∙∙ + 3OxO

Page 2 of 25

(1)

In contrast, doping ceria with zirconia has no effect on the number of ionic charge carriers (cf. eq. 2). Nevertheless, the influence of zirconia on the properties of doped ceria is of interest for several reasons: 1) Ceria doped with zirconia exhibits excellent performance in the production of solar fuels as well as in catalytic processes such as in automotive three-way catalysts due to the increased oxygen storage capacity (OSC).3-6 2) In solid electrochemical cells, doped ceria is used as a buffer layer between yttria stabilized zirconia and the cathode material and a zirconia/ceria solid solution might form at the interface during operation. However, this solid solution shows a lower ionic conductivity than the initial materials and thus lowers the efficiency of the cell.7 3) Zirconia is a common impurity in ceria and could alter properties in nominally undoped ceria due to the interaction of Zr-ions with other defects.8 CeO2

ZrO2 →

x ZrCe +2OxO

(2)

In previous computational studies, a strong interaction between zirconium ions and oxygen vacancies was found and consequently even small amounts of zirconium can affect properties like the conductivity.9 For example, Eufinger et al.7 showed that the conductivity in Y-doped ceria is significantly decreased by the addition of zirconia. The ionic conductivity in doped ceria can be predicted by Kinetic Monte Carlo (KMC) simulations as demonstrated in our earlier publications.10-11 Since the migration energy of each individual jump depends on the local environment, a suitable energy model with an appropriate interaction range has to be applied. Considering only the nearest ions around the initial, final and transition state of the migrating oxygen ion leads to the ‘6-cation environment’ shown in Fig. 1. Each cation position can be occupied by a cerium ion, zirconium ion or yttrium ion leading to 36=729 configurations, where 270 unique jump configurations are possible due to symmetry. In the process of one individual jump, the oxygen ion moves from the initial to the final state with the configurational energy difference ∆𝐸conf (depending on the occupation of position 3-6 in Fig. 1) and passes through the edge formed by two cation sites (positions 1 and 2 in Fig. 1), which could be occupied by either cerium or dopant ions with an associated migration edge energy 𝐸edge . Therefore, the edge energy is defined as the migration energy for jump environments 2 ACS Paragon Plus Environment

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where only the edge positions are doped and all other positions are undoped. A first approximation to the migration energy for any configuration is given by:10-14

𝐸mig = 𝐸edge +

∆𝐸conf 2

(3)

This model is based on the consideration that the energy difference between the initial and final state can be linearly interpolated and the migration barrier, which depends on the edge occupation, is added at the geometric and energetic center of the jump (cf. Fig. S1). Additionally, the migration energy is limited to positive values (𝐸mig ≥ 0), although eq. 3 could predict negative values. In a further approximation, the contribution ∆𝐸conf can be estimated from the pair interactions of the defects, ∆𝐸conf =

∑ shell 𝑖 species 𝑗

final 𝑁𝑖,𝑗 𝜀𝑖,𝑗 −

∑ shell 𝑖 species 𝑗

initial 𝑁𝑖,𝑗 𝜀𝑖,𝑗

(4)

where 𝑁𝑖,𝑗 and 𝜀𝑖,𝑗 are the number of interactions and the interaction energies, respectively and the sums includes all possible interaction types and could include further distances beyond the nearest neighbor positions. This model allows a fast evaluation of the migration energy 𝐸mig with a limited number of parameters and automatically retains the energy difference between the initial and final state. Conversely, the model has some limitations. Firstly, for a symmetric configuration, i.e. the same number of doping ions on both sides of the migration edge, the migration energy is always equal to the edge energy, which cannot be expected to be the general case. Secondly, the transition state is always assumed to be in the geometrical and energetic center between initial and finial state, which is unlikely for asymmetric configurations.

Fig. 1: 6-cation environment (yellow spheres) around the jump of an oxygen ion (red sphere) to an oxygen vacancy (red cube). Numbers correspond to the nomenclature in the text. Positions 1 and 2 form the migration edge that has to be passed during the ion jump. Positions 3-6 are the nearest neighbor positions (NN) to initial and final position of the vacancy.



Page 4 of 25

In this paper, we check the validity of the model for Zr/Y-co-doped ceria by explicitly calculating the energy barriers for different occupations of the six cations around the jump using density functional theory (DFT). 178 of the 270 possible jump configurations are calculated in this study. The individual configurations will be termed according to the numbering in Fig. 1, e.g. the configuration CeZr-YCeCeY, where one Zr-ion is at the migration edge and a Y-ion is on both the left and the right side of the 6-cation environment. The calculated energies are applied in Kinetic Monte Carlo simulations to assess the impact of Zr-doping on the ionic conductivity of Y-doped ceria.

II Computational methods 1. Density functional theory All DFT+U calculations were performed with the Vienna Ab initio Simulation Package15-16 (VASP) in the generalized gradient approximation (GGA) according to Perdew, Burke, and Ernzerhof17 (PBE) with the projected augmented wave (PAW) method.18 The energy cutoff of the plane waves was set to 500 eV and the convergence parameters for electronic and ionic relaxation were set to 10-4 eV and 0.01 eV/Å, respectively. For the 2 × 2 × 2 supercell of the ceria fluorite structure a 3 × 3 × 3 Monkhorst-Pack k-point mesh was applied.19 The 2s22p4, 5s25p66s25d14f1, 4s24p65s24d1, and 4s24p65s24d2 electrons of oxygen, cerium, yttrium and zirconium atoms were treated as valence electrons, respectively. A Hubbard U-parameter was introduced to account for the localization of strongly correlated electrons by the simplified rotational invariant approach with U = 5 eV for the 4f-orbitals of cerium.20 The lattice parameter of the unit cell was fixed at the value of bulk ceria (5.49 Å) according to previous calculations.10 The total number of electrons in the cell was adapted for all defective cells to reproduce the actual charge state of the defects, e.g. (Ce31YO63)+ for a 2 × 2 × 2 supercell containing one Y ion and one oxygen vacancy. Though charge-neutral cells containing defects according to eqs. 1 and 2 without adjustment of the number of electrons would be preferable; in this work interactions between defects shall be limited. Therefore, charge-neutral cells with large distances between defects are virtually divided into oppositely charged cells. Charged cells are calculated by VASP assuming a neutralizing background charge, which is a valid approach as shown in literature.21-23 The migration energies 4 ACS Paragon Plus Environment

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


were determined by the nudged elastic band – climbing image (NEB-CI) method24 after initial and final state were fully relaxed.

2. Monte Carlo The distribution of the defects was simulated by the Metropolis Monte Carlo (MMC) algorithm19 as described in an earlier publication.9 The lattices were constructed from 14 × 14 × 14 ceria unit cells and the defect distributions were equilibrated until the averaged energy was constant. After the equilibration, the configurations were sampled over 106 subsequent steps. The ionic conductivities were simulated by Kinetic Monte Carlo (KMC) using the in-house program iCon. The lattices were constructed from 14 × 14 × 14 ceria unit cells. The attempt frequency was set to a commonly accepted value of 1013 Hz.8, 25 After a prerun of 200 Monte Carlo steps per oxygen ion (MCSP), the conductivities were subsequently recorded over 1000 MCSP. Results were averaged over ten independent simulations. Pair interactions applied in MMC and KMC simulations were cut off at a distance of 5.49 Å and are given in Table 1 and Table 2. The energies are shifted in a way that the energy of the first shell beyond the interaction radius is set to zero.11 For the contribution of the 6-cation environment, three different models were tested while contributions of all further distances were included by pair interactions according to eq. 4. The simulated system is described as a solution of dopant ions and oxygen vacancies in a matrix of cerium oxide. Therefore, it is expected that the model loses accuracy at high dopant concentration and thus the total fraction of dopants is limited to be below 0.5. Oxidizing conditions were considered for all simulations, with the result that no reduction of the lattice takes place and the concentration of oxygen vacancies is purely due to doping. Table 1: Pair interaction between dopant ions and an oxygen vacancy. x ′ interaction shell distance (Å) ZrCe -VO∙∙ interaction (eV) YCe -VO∙∙ interaction (eV)

1NN

2.38

-0.60

-0.27

2NN

4.55

-0.05

-0.06



Page 6 of 25

Table 2: Pair interaction between two oxygen vacancies. Interactions 3NNa and 3NNb are without or with a cation on the interaction axis, respectively.

interaction shell distance (Å) VO∙∙ -VO∙∙ interaction (eV) 1NN

2.75

0.84

2NN

3.88

0.27

3NNa

4.75

0.26

3NNb

4.75

0.40

4NN

5.49

0.28

III Results and Discussion 1. Defect configurations For lower doping levels, one may argue that only configurations with a very limited number of dopants occur. However, the interaction of oxygen vacancies with the dopants can lead to a clear deviation from the random distribution. We checked this hypothesis by performing Metropolis Monte Carlo (MMC) simulations for a doping level of 𝑥Zr = 0.1 and 𝑦Y = 0.2. Cations were distributed randomly and the distribution of oxygen ions and vacancies was simulated at 1000 K applying the MMC algorithm. The fractions of the five most frequent configurations are given in Fig. 2 in comparison to the results for the random distribution. The results clearly show that the configurations with at least one zirconium ion are the most frequent ones. The fractions of these jump configurations clearly deviate from the random distribution due to the strong interaction between dopant ions and oxygen vacancies.


Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 2: Fraction of most frequent jump configurations in thermodynamic equilibrium simulated at 1000 K and in the random distribution for 𝒙𝐙𝐫 = 𝟎. 𝟏 and 𝒚𝐘 = 𝟎. 𝟐.

2. Edge energies Fig. 3 shows the edge energies for different configurations of the edge. For the calculation of the edge energy, the occupation of the edge is varied while all other cations are cerium ions. The energies were calculated for the volume of the bulk cell (constant volume) as well as the relaxed volume of the initial and final state (constant pressure). The energies show a clear dependence on the sum of the ionic radii of the edge ions with larger ionic radii leading to higher edge energies. This finding is in agreement with the early postulation by Kilner and Brook26 and also the results by Plata et al.,27 who calculated migration energies at the ceria-zirconia interface. However, the results are in contrast to the assumption of Schriever et al.,28 who stated that the energy barrier depends on the charge of the edge ions rather than the radius. The relaxation of the volume has some effect on the migration barrier especially for lower radii, but the general trend is not changed.



Page 8 of 25

29

Fig. 3: Edge energies for different configurations depending on the combined ionic radius of the edge ions. The energies are calculated at the volume of the bulk cell or the relaxed volume of the initial state.

3. Migration configurations Calculation of ∆𝑬𝐜𝐨𝐧𝐟 In the model, the energy difference ∆𝐸conf is calculated from pair interaction energies. To check the validity of this assumption, Fig. 4 shows the energy difference ∆𝐸conf for the configurations with Ce-Ce edge between the initial and final state as calculated by DFT and by the pair interaction model using energies from Table 1. The agreement of the model with a maximum deviation of 0.06 eV for one of the Zr-rich configurations is remarkable. Consequently, the energy difference ∆𝐸conf in eq. 3 can be well estimated by the corresponding pair interaction energies in agreement with our previous findings.11


Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 4: ∆𝑬𝐜𝐨𝐧𝐟 for configurations with a Ce-Ce edge as calculated by DFT and the pair interaction model

Migration energies for Ce-Ce edge The migration energies for the different configurations calculated by DFT and the ‘standard’ model according to eq. 3 are shown in Fig. 5 for the Ce-Ce edge. The configuration CeCeZrZrCeCe was not calculated as the initial state is not stable in this case. This is due to the strong association of the oxygen vacancy with the two zirconium ions. Fig. 5 reveals that the model reproduces the trend of the DFT migration energies correctly, but predicts higher migration energies than calculated by DFT. In addition, the configurations without Zr-ions are better predicted than the configurations with Zr-ions and the largest deviation is found for the configuration CeCe-YZrZrY. The overall standard deviation for the model is 0.15 eV. The larger deviation for Zr-containing configurations is probably due to the strong association between zirconium ions and oxygen vacancies, which leads to a strong distortion of the local structure.9 A closer look at Fig. 5 and the associated configurations reveals further limitations of the model. It is apparent that a large number of configurations have a model energy of 0.52 eV but the corresponding DFT energies vary between 0.22 eV and 0.52 eV. These configurations are all symmetric with the same number of dopant ions in initial and final state and, therefore, the migration energies are equal to the Ce-Ce edge energy according to the model. However, the 9 ACS Paragon Plus Environment


Page 10 of 25

DFT energies are very divers and consequently the migration barriers are not only influenced by the energy difference between the states. Within this group the configurations CeCe-ZrCeZrCe and CeCe-ZrCeCeZr exhibit a special behavior. These configurations only differ by the fact that the two Zr ions are either in ‘cis’ (ZrCeCeZr) or ‘trans’ (ZrCeZrCe) position as shown in Fig. 6. Nevertheless, the migration energy for the cis configuration is 0.25 eV lower than for the trans configuration. This difference can be explained by the interaction of the zirconium ion with the oxygen ion at the transition state. The Zr ion shows an attractive interaction with the oxygen vacancy or, in different words, repels the oxygen ions. In the cis configuration, the migrating oxygen ion follows a curved path and the transition state is stabilized by a larger distance to both zirconium ions (Fig. 6 left). This stabilization is not possible in the trans case where the migration path is almost linear and the migration energy is consequently higher (Fig. 6 right). This effect can be found for other cis/trans configurations as well but is not as pronounced as in the discussed case.

Fig. 5: Comparison of the migration energies for the Ce-Ce edge as calculated by DFT and the model.


Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 6: Configuration ZrCeCeZr (cis) and ZrCeZrCe (trans) and the path of the oxygen ion between initial and final state. Cerium, zirconium and oxygen ions are represented in yellow, green and red respectively.

The impact of finite size effects on the migration energy has been discussed before and it was shown that the migration energy depends on the size of the calculated supercell.11 To investigate this effect, the same calculations were performed in 3 × 3 × 3 supercells. Some differences in the energies were found between 2 × 2 × 2 and 3 × 3 × 3 supercell but the general trends remain the same (cf. Fig. S2). Therefore, further edges were calculated within 2 × 2 × 2 supercells.

Further edges Calculations for the other edges were performed for a total of 134 jump configurations. The comparison of the DFT results and model predictions is given in Fig. 7. The combined results are similar to the results for the Ce-Ce edge and the model generally overestimates the migration energies with a standard deviation of 0.149 eV. Exceptions are mainly found for very small migration energies. In these cases, ∆𝐸conf has a negative sign and has a large contribution compared to the migration edge. Since the transition state is always assumed to be located in the center between initial and final state, the resulting migration energy is zero (due to the constrain 𝐸mig ≥ 0). For comparison the migration path and energy profile of YZr-YZrCeCe are shown in Fig. 8. The transition state is clearly shifted to the higher energy state as expected from Hammond’s postulate. This results in a small albeit positive migration energy in one direction, whereas it is predicted to be zero by the model.



Page 12 of 25

Fig. 7: Comparison of migration energies for doped edges as calculated by DFT and the model.

Fig. 8: Migration path (left) and relative energy (right) for configuration YZr-YZrCeCe. Cerium, zirconium, yttrium and oxygen ions are represented in yellow, green, gray and red respectively. The energy is given relative to the final state and the oxygen ion position is projected to the vector between the initial and final position of the ion. The line is a guide to the eye only.

As shown above, the migration energy depends on the migration edge and the cation occupation of the NN positions for the initial and final state. In the standard model, the second term is approximated by pair interactions. Alternatively, we can extract this contribution for every α,CeCe configuration α from the explicitly calculated migration energies 𝐸mig of the Ce-Ce edge by

subtracting the edge energy 𝐸edge,CeCe. Within this approach, we can then calculate the migration energy for every configuration α with the edge 𝑋𝑌 from this contribution and the corresponding edge energy:

α,𝑋𝑌 α,CeCe 𝐸mig = 𝐸edge,𝑋𝑌 + (𝐸mig − 𝐸edge,CeCe )

(5)


Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The results of this model approach are shown in Fig. 9. The modelled energies are in good agreement with the calculated DFT values with an overall standard deviation of 0.058 eV. We would like to point out that a larger number of configurations lead to a migration energy of 0 eV (due to the constrain 𝐸mig ≥ 0). The reason for this finding is the low migration edge energies of Zr-containing edged leading to a large energy difference compared to the Ce-Ce edge.

Fig. 9: Comparison of migration energies for different configurations as calculated directly by DFT and based on the Ce-Ce migration energies for all migration edges as described in the text.

This result shows that it is indeed possible to separate the influence of the migration edge and the NN positions, i.e. the two contributions do not depend on each other. This might be due to the fact that the geometry of these positions is different with regard to the migration path. While the edge cations are next to the center of the jump, the NN positions are next to the initial and final state.

4. Structural analysis The standard model can reproduce the correct trends and the energy difference ∆𝐸conf can be approximated by the pair interaction energy. However, the model overestimates the migration barriers. To investigate this deviation in more detail, the relation between composition, structure, and energy is analyzed. For the analysis, a variety of compositional and structural descriptors is extracted from the calculated structures given in Table S1 in the supporting information. 13 ACS Paragon Plus Environment


Page 14 of 25

The descriptors are related to the migration energies by a decision tree algorithm30 with 140 randomly selected migration energies serving as training set and 36 energies as a test set. The results of the fit are given in Fig. 10. The comparison between the DFT energies and the model shows a good representation of the energies by the model. In terms of predictive power, the results are limited as, on the one hand, the predictors include parameters accessible only from NEB calculations and, on the other hand, the formal relation between descriptors and target value, i.e. migration energy, is not quite clear. Nevertheless, the fit also gives the relative importance of the predictors for the model as shown on the right side of Fig. 10. The most important descriptors are the migration length between initial and transition state (a9), i.e. the distance that the oxygen ions move between these states, and the displacement of the oxygen ion in the initial state (a7). Of course, these two quantities are connected as a displacement of the oxygen ion also leads to a change of the migration length.

Fig. 10: Left: Comparison between calculated values and fit values. Right: Relative importance of descriptors.

For a detailed analysis, the relation between migration length and migration energy for the Ce-Ce edge is shown in Fig. 11. Except for some zirconium containing configurations, a clear relation between both quantities is found with increasing migration length leading to increasing migration energy. Therefore, configurations with a small distance between initial and final state lead to small migration energies. The change of the migration length is due to two effects. On the one hand, the transition state shifts to the state with the higher energy as seen in Fig. 8. On the other 14 ACS Paragon Plus Environment

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hand, dopants on the positions 3 to 6 lead to a compression of the octahedron and a shortage of the migration path. Similar results as in Fig. 11 are found for the other migration edges.

Fig. 11: Calculated migration energy depending on the migration distance between initial state and transition state.

5. Extended model The analysis of the configurations showed that the standard model (eq. 3) overestimates the migration energies and that additional doping leads to a shortage of the migration length and thereby a decrease of the migration energy. For example, the standard model predicts a migration value of 0.52 eV for the configurations CeCe-CeCeCeCe, CeCe-YCeYCe and CeCe-YYYY while the DFT calculations lead to values of 0.52 eV, 0.43 eV and 0.35 eV, respectively. Therefore, an increase of the number of dopants in NN position leads to a decrease of the migration barrier for forward and backward jump. An extension of the model is possible by including not only the migration edge energy and the energy difference between initial and final state in the model but also the total number of dopants in the NN positions. The migration energy for each configuration is then given by:

𝐸mig = 𝐸edge +

∆𝐸conf + 𝑁Y 𝜀Y + 𝑁Zr 𝜀Zr 2

(6)



where 𝑁Y /𝑁Zr are the number of respective dopant ions on NN positions and

Page 16 of 25

𝜀Y /𝜀Zr the

corresponding energy contributions. We extracted the values 𝜀Y and 𝜀Zr from the configurations CeCe-YYYY and CeCe-ZrZrZrZr leading to values of -0.042 eV and -0.059 eV, respectively. The resulting model energies are shown in Fig. 12. The extension of the model clearly improves the energies compared to the standard model giving an overall standard deviation of 0.061 eV. We emphasize here that the values for 𝜀Y and 𝜀Zr were not obtained by fitting all the configurations but extracted only from two additional configurations. In total the model includes three contributions: The migration barrier depending on the occupation of the migration edge, the contribution due to the energy difference between initial and final state and the contribution to the migration barrier due to the total number of dopant ions in NN positions, which lead to a distortion of the local structure. It should be noted that the contribution ∆𝐸conf determines the energy difference of the initial and final states while the other contributions only influence the height of the transition state relative to these states. For the 6-cation environment a total of ten energy parameters are necessary; including six edge energies, two pair interactions and two energy contributions for the total number of dopants.

Fig. 12: Comparison of migration energies for doped edges as calculated by DFT and the extended model.


Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Simulated conductivity The ionic conductivity in Y/Zr-doped ceria was simulated with the Kinetic Monte Carlo method applying three different models: -

‘Standard’: The migration energy is calculated according to eq. 3 with the individual edge energies and further contributions from pair interactions.

-

‘Extended’: The migration energy for the 6-cation environment is calculated according to eq. 6 including the consideration of the total number of dopant ions around initial and final position of the migrating oxygen ion. Further contributions are added by pair interactions.

-

‘Explicit’: The migration energy for the 6-cation environment is included explicitly from the calculated configurations. Configurations not calculated are approximated by eq. 5. All other contributions are calculated from pair interactions.

Fig. 13 shows the results for the three models at 773 K depending on the yttrium fraction without Zr-doping. Each curve shows the maximum in ionic conductivity as observed in experiments.31 This maximum originates from competing effects: On the one hand, the conductivity increases due to the formation of additional oxygen vacancies with doping. On the other hand, the additional dopant ions lead to ‘trapping’, i.e. the vacancies associate with the dopant ions, and ‘blocking’, i.e. the edge energy increases if dopant ions are located at the migration edge.11 Compared to the standard model, the extended and explicit models show a shift of the maximum to higher yttrium fractions and an increase in the conductivity. The extended model leads to similar results as the more sophisticated explicit model. In comparison to the experimental values by Zhang et al.31 the simulations overestimate the conductivity. This could be due to the selection of the attempt frequency, which is taken here as a commonly accepted value of 1013 Hz.8 Moreover, further experimental literature data for the bulk conductivity in Y-doped ceria shows scattering of the ionic conductivity.11 Interestingly, the standard model results in the best agreement with the experimental data. One reason for this result might be the chosen supercell size and its resulting energy parameters for the KMC simulations as well as the use of the commonly excepted attempts frequency of 1013 Hz.11



Page 18 of 25

Fig. 13: Simulated conductivity depending on the yttrium fraction at 773 K using the standard, extended and explicit 31 migration energy model in comparison with experimental values by Zhang et al.

In Fig. 14, the results for the conductivity of ceria doped with a constant dopant fraction of 0.2 yttrium and a varying amount of zirconium are shown. For all models the conductivity first decreases with increasing zirconium fraction. This is due to the strong trapping of the oxygen vacancies by the zirconium ions. However, for the explicit model, the conductivity strongly increases for zirconium fractions above 0.05 in contrast to experiment.7 The reason is the low energy barrier for oxygen migration of configurations containing zirconium, i.e. a strong ‘negative blocking’ effect. A similar but much weaker effect is observed for the standard model with a slight increase of the conductivity for zirconium fractions above 0.15. Despite including the strong negative blocking effect, the extended model does not agree well with the explicit model for higher zirconium fractions in contrast to the simulations without zirconium doping (Fig. 13). Instead, for zirconium fractions above 0.05, the conductivity decreases linearly. This can be explained by the large number of invalid jump attempts. Due to the energy model, some jumps can lead to configurations where the probability of the return jump is one, i.e. where no migration barrier exists. These configurations are regarded as unstable and the jump attempts leading to this type of configuration are automatically rejected. However, invalid jump attempts still increase the simulation time. The fraction of invalid jump attempts is given in Figure S3 in the supporting information, showing that the extended model exhibits a 18 ACS Paragon Plus Environment

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


large fraction of invalid jump attempts compared to the other models. This behavior originates in the fact that the additional energy correction in the extended model often leads to vanishing migration energies for many jump configurations especially at high zirconium fractions. In case a model predicts a large number of migration barriers incorrectly as zero, a large number of invalid jump attempts occurs and the motion of the oxygen ions is not described correctly anymore. Consequently, this results in the behavior observed in Fig. 14 for the extended model.

Fig. 14: Simulated conductivity depending on the level of zirconium doping for ceria doped with 0.2 yttrium at 773 K using the standard, extended and explicit migration energy model.

None of the applied models is able to describe the experimental behavior correctly. This finding suggests that another effect besides trapping and blocking influences the migration energies. Indeed, the small size of the zirconium ion leads not only to a strong trapping and a negative blocking effect, but it also has an effect on the lattice parameter as the calculated lattice parameter of ZrO2 (5.15 Å) is considerably smaller than the calculated lattice parameter of CeO2 (5.49 Å). Therefore, zirconia doping decreases the lattice parameter. However, our DFT calculations were all performed using the lattice parameter in defect-free ceria. Since the lattice parameter is known to have a clear impact on the migration energy,32 this effect is considered in the following for all migration energy models. For this purpose, we calculated the migration energies with a lattice parameter of 5.39 Å corresponding to a zirconium content of 0.3 according to Vegard’s law. Further test calculations showed that the migration energy depends linearly on the lattice parameter for 0 < 𝑥 < 0.3. Therefore, all edge energies were calculated for the smaller 19 ACS Paragon Plus Environment


Page 20 of 25

lattice parameter and, in the following, the edge energies were described as a linear function of the zirconium content. For the migration energies of the explicit configurations and all the pair interactions, a set of configurations with the smaller lattice parameter was calculated and a linear model was fitted to predict the migration energies for all configurations depending on the zirconium content. More details can be found in the supporting information. In the new simulations the ionic conductivities were then simulated with migration energies explicitly depending on the zirconium fraction. Additional influence of the lattice parameter on the calculation of the defect concentration and the jump length in the KMC simulations were neglected as the effect on the conductivity is below 2%. The results of the simulations are shown in Fig. 15. With the new models, the conductivities decrease with increasing zirconium content in the considered range in agreement with the experimental findings by Eufinger et al.7 Once again, the standard model shows the best agreement with the experimental values. However, all simulated curves show a similar form and vary mainly by a shift in the absolute values. Since the commonly excepted value of 1013 Hz was used for the attempt frequency, a shift along the y-axis is possible if the actual attempt frequency is different from this value. Furthermore, it should be noted that additional effects might be important, such as the influence of doping on the attempt frequency or thermal effects on the migration energies.8, 33 We also neglected the influence of yttrium doping on the lattice parameter since the change in the considered doping range is much smaller than for zirconium doping. In this context, only small differences between the models are found and all models show the same trend in agreement with experiment. We emphasize that the experimental decrease in ionic conductivity with increasing zirconia doping can only be reproduced by including the change of the lattice parameter with dopant fraction


Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 15: Conductivities simulated with the volume depending models for different levels of zirconium doping for ceria doped 7 with 0.2 yttrium at 1000 K. Experimental values by Eufinger et al. are provided for comparison.

IV Conclusion We calculated the migration energy of oxygen ions for 178 out of 270 unique occupations of the 6-cation environment around the jump center in Y/Zr-doped ceria. The results are compared to different models, which estimate the migration energies by a limited set of parameters. The ‘standard’ model combines the explicitly calculated ‘edge’ energies with the energy difference between initial and final state obtained from pair interactions. Eight energy parameters are required for this model. However, the strong distortion of the lattice, due to the small radius of the zirconium ion, leads to a pronounced deviation between explicitly calculated DFT results and the standard model. The structural analysis shows a clear relation between the distance that the oxygen ion has to travel between initial and final state and the corresponding migration energy, which is due to the local distortion introduced by the dopant ions. From the analysis of the calculated data it is apparent that the migration energies do not only depend on the migration edge and the energy difference between initial and final state but also on the total number of dopants at initial and final position. Therefore, an ‘extended model’ is suggested, which additionally includes the total number of dopants at the NN position of initial and final state This model leads to a better agreement of modeled and explicitly calculated energy 21 ACS Paragon Plus Environment


Page 22 of 25

barriers. In the extended model, ten energy parameters are required to describe the 270 possible configurations for the 6-cation environment with good accuracy. In addition, we could show that it is possible to group jump configurations according to their migration edge. As a result, all jump configurations can be calculated by combining the jump configurations with a Ce-Ce edge and the edge energies (the ‘explicit’ model). Ionic conductivities were simulated using the KMC method with the standard, extended and explicit model to describe the influence of the 6-cation environment and additional pair interactions to include further distances. The simulation of the ionic conductivity of Y-doped ceria produces a maximum in the conductivity as a function of the dopant fraction corresponding to experiments, which can be ascribed to ‘trapping’ and ‘blocking’ effects. For Y/Zr-co-doped ceria, simulations with the explicit model and a constant lattice parameter show an increase of the conductivity with increasing zirconium fraction for 𝑥 > 0.05 in Ce0.8-xZrxY0.2O1.9 in contrast to experimental observations. The origin of this behavior is the ‘negative blocking’ effect, i.e. the decrease of the migration energy for zirconium containing environments. The experimental finding is correctly reproduced only if the decrease of the lattice parameter due to doping with the smaller zirconium ions is taken into account in the migration energy model. Our study shows that besides trapping and blocking the change of the lattice parameter due to doping has a substantial influence on the ionic conductivity and has to be included if doping considerably changes the lattice parameter.

Supporting Information Additional data, as mentioned in the text, is given in “supporting-information.pdf”.

Acknowledgment The authors like to thank Randy Jalem from the National Institute for Materials Science, Tsukuba, Japan for the creation of the decision tree model. The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARAHPC Partition part of the supercomputer CLAIX at RWTH Aachen University. Structural models in Fig. 1, Fig. 6, and Fig. 8 were created with VESTA.34 22 ACS Paragon Plus Environment

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References 1.

2. 3.

4.

5.

6. 7.

8.

9. 10. 11.

12. 13. 14. 15. 16. 17. 18. 19. 20.

21.

Chueh, W. C.; Falter, C.; Abbott, M.; Scipio, D.; Furler, P.; Haile, S. M.; Steinfeld, A., High-Flux Solar-Driven Thermochemical Dissociation of CO2 and H2O Using Nonstoichiometric Ceria. Science 2010, 330, 1797–1801. Steele, B. C. H.; Heinzel, A., Materials for Fuel-Cell Technologies. Nature 2001, 414, 345–352. Hao, Y.; Yang, C.-K.; Haile, S. M., Ceria–Zirconia Solid Solutions (Ce1– XZrxO2−δ , x ≤ 0.2) for Solar Thermochemical Water Splitting: A Thermodynamic Study. Chem. Mater. 2014, 26, 6073– 6082. Call, F.; Roeb, M.; Schmücker, M.; Bru, H.; Curulla-Ferre, D.; Sattler, C.; Pitz-Paal, R., Thermogravimetric Analysis of Zirconia-Doped Ceria for Thermochemical Production of Solar Fuel. Am. J. Analyt. Chem. 2013, 04, 37–45. Call, F.; Roeb, M.; Schmücker, M.; Sattler, C.; Pitz-Paal, R., Ceria Doped with Zirconium and Lanthanide Oxides to Enhance Solar Thermochemical Production of Fuels. J. Phys. Chem. C 2015, 119, 6929–6938. Trovarelli, A., Catalytic Properties of Ceria and CeO2 -Containing Materials. Cat. Rev. 1996, 38, 439–520. Eufinger, J.-P.; Daniels, M.; Schmale, K.; Berendts, S.; Ulbrich, G.; Lerch, M.; Wiemhöfer, H.-D.; Janek, J., The Model Case of an Oxygen Storage Catalyst - Non-Stoichiometry, Point Defects and Electrical Conductivity of Single Crystalline CeO2-ZrO2-Y2O3 Solid Solutions. Phys. Chem. Chem. Phys. 2014, 16, 25583–25600. Koettgen, J.; Zacherle, T.; Grieshammer, S.; Martin, M., Ab Initio Calculation of the Attempt Frequency of Oxygen Diffusion in Pure and Samarium Doped Ceria. Phys. Chem. Chem. Phys. 2017, 19, 9957-9973. Grieshammer, S., Defect Interactions in the CeO2-ZrO2-Y2O3 Solid Solution. J. Phys. Chem. C 2017, 121, 15078-15084. Grieshammer, S.; Grope, B. O. H.; Koettgen, J.; Martin, M., A Combined DFT U and Monte Carlo Study on Rare Earth Doped Ceria. Phys. Chem. Chem. Phys. 2014, 16, 9974. Koettgen, J.; Grieshammer, S.; Hein, P.; Grope, B.O.H.; Nakayama, M.; Martin, M., Understanding the Ionic Conductivity Maximum in Doped Ceria: Trapping and Blocking. Phys. Chem. Chem. Phys. 2018, 20, 14291-14321. Lee, E.; Prinz, F.; Cai, W., Enhancing Ionic Conductivity of Bulk Single-Crystal Yttria-Stabilized Zirconia by Tailoring Dopant Distribution. Phys. Rev. B 2011, 83. Ven, A. V. d.; Ceder, G.; Asta, M.; Tepesch, P., First-Principles Theory of Ionic Diffusion with Nondilute Carriers. Phys. Rev. B 2001, 64, 184307. Ven, A. V. d.; Ceder, G., First Principles Calculation of the Interdiffusion Coefficient in Binary Alloys. Phys. Rev. Lett. 2005, 94, 45901. Kresse, G.; Furthmüller, J., Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. Kresse, G., From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758–1775. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1997, 78, 3865–3868. Blöchl, P. E., Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. Monkhorst, H. J.; Pack, J. D., Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188–5192. Dudarev, S. L.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P., Electron-Energy-Loss Spectra and the Structural Stability of Nickel Oxide:An Lsda+U Study. Phys. Rev. B 1998, 57, 1505– 1509. Zacherle, T.; Schmidt, P. C.; Martin, M., Ab Initio Calculations on the Defect Structure of βGa2O3. Phys. Rev. B 2013, 87. 23 ACS Paragon Plus Environment


22.

23.

24. 25. 26. 27.

28.

29. 30. 31. 32. 33. 34.

Page 24 of 25

Zheng, J. X.; Ceder, G.; Maxisch, T.; Chim, W. K.; Choi, W. K., Native Point Defects in Yttria and Relevance to Its Use as a High-Dielectric-Constant Gate Oxide Material: First-Principles Study. Phys. Rev. B 2006, 73. Li, X.; Finnis, M. W.; He, J.; Behera, R. K.; Phillpot, S. R.; Sinnott, S. B.; Dickey, E. C., Energetics of Charged Point Defects in Rutile TiO2 by Density Functional Theory. Acta Mater. 2009, 57, 5882-5891. Henkelman, G.; Uberuaga, B. P.; Jónsson, H., A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901–9904. Grope, B. O. H.; Zacherle, T.; Nakayama, M.; Martin, M., Oxygen Ion Conductivity of Doped Ceria: A Kinetic Monte Carlo Study. Solid State Ionics 2012, 225, 476–483. Kilner, J.; Brook, R., A Study of Oxygen Ion Conductivity in Doped Non-Stoichiometric Oxides. Solid State Ionics 1982, 6, 237–252. Plata, J. J.; Márquez, A. M.; Sanz, J. F., Understanding the Interplay of Dopants, Interfaces, and Anionic Conductivity in Doped Ceria/Zirconia Heteroepitaxial Structures. Chem. Mater. 2014, 26, 3385–3390. Genreith-Schriever, A. R.; Hebbeker, P.; Hinterberg, J.; Zacherle, T.; Souza, R. A., Understanding Oxygen-Vacancy Migration in the Fluorite Oxide CeO2: An Ab Initio Study of Impurity-Anion Migration. J. Phys. Chem. C 2015, 119, 28269–28275. Shannon, R. D., Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides. Acta Crystall. A 1976, 32, 751–767. Rokach, L.; Maimon, O., Data Mining with Decision Trees: Theory and Applications, Second edition. ed.; World Scientific: Hackensack, New Jersey, 2015, p xxi, 305 pages. Zhang, T.; Ma, J.; Kong, L.; Chan, S.; Kilner, J., Aging Behavior and Ionic Conductivity of CeriaBased Ceramics: A Comparative Study. Solid State Ionics 2004, 170, 209–217. Hinterberg, J.; Zacherle, T.; Souza, R. A., Activation Volume Tensor for Oxygen-Vacancy Migration in Strained Ceo2 Electrolytes. Phys. Rev. Lett. 2013, 110, 205901. Koettgen, J.; Schmidt, P. C.; Bucko, T.; Martin, M., Ab Initio Calculation of the Migration Free Energy of Oxygen Diffusion in Pure and Samarium-Doped Ceria. Phys. Rev. B 2018, 97, 024305. Momma, K.; Izumi, F., Vesta. J. Appl. Crystallogr. 2008, 41, 653–658.


Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


Modeling Oxygen Ion Migration in the CeO2-ZrO2-Y2O3 Solid Solution

Recommend Documents