Understanding the Enhancement of Ionic Transport in

Publication Date (Web): September 10, 2018. Copyright © 2018 American ... Abstract. Zirconia-based ceramics have been the most promising oxide electr...
0 downloads 0 Views 681KB Size
Subscriber access provided by University of South Dakota

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Understanding the Enhancement of Ionic Transport in Heterogeneously-Doped Zirconia by Heterointerface Engineering Mehmet Emin Kilic, and Aloysius Soon J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b05111 • Publication Date (Web): 10 Sep 2018 Downloaded from http://pubs.acs.org on September 15, 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 41 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

Understanding the Enhancement of Ionic Transport in Heterogeneously-Doped Zirconia by Heterointerface Engineering Mehmet Emin Kilic and Aloysius Soon∗ Department of Materials Science and Engineering, Yonsei University, Seoul 03722, Korea E-mail: [email protected] Abstract Zirconia-based ceramics have been the most promising oxide electrolyte material with high ionic conductivity for solid oxide fuel cell (SOFC) applications. Even though yttria-stabilized zirconia (YSZ) and scandia-stabilized zirconia (ScSZ) are typically used for the SOFC at high temperatures, their performance is not optimal at operating temperatures with respect to its ionic conductivity and stability. Previous literature have focused largely on ionic diffusion dynamics in bulk YSZ and ScSZ, while their heterogeneously-doped alloy and heterolayered superlattices are less investigated. In this work, using molecular dynamics simulations and diffusion dynamics analysis, we examine and consider five main mechanisms that may contribute to the enhancement of the overall ionic conductivity of these doped zirconia, namely the influence of cation size, concentration, distribution, the crystal orientation and direction, and lastly, the degree of atomic roughness at the interface in the heterolayered structures. Our results support that heterointerface engineering at the atomic scale greatly reduces local lattice distortions (commonly seen in the bulk phases) while inducing an in-plane strain, and thus leading to an overall enhancement of the ionic conductivity and stability for SOFC applications.

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

1

Introduction

Solid oxide fuel cells (SOFC), as electrochemical devices that convert chemical energy into electrical energy, consist of three main elements – the cathode, the electrolyte, and the anode. The electrolyte must be stable under reducing and oxidizing environments, in addition to having high ionic conductivity with low electronic conductivity. Zirconia-based ceramics have been the most favored electrolyte with high ion conductivity for SOFC. 1 Pure zirconia (ZrO2 ) is monoclinic (m) at low temperatures. Upon increasing the temperature, the material exhibits phase transitions from the m phase to the tetragonal (t) phase, and then the cubic (c). Eventually, melting occurs at 2870 K. 2 These phase transitions, under SOFC operating conditions, give rise to drastic changes in both volume and conductivity, which makes the pure material unsuitable for applications. For most practical applications, it is often doped with other oxides, e.g. yttria (Y2 O3 ) and scandia (Sc2 O3 ), to aid stabilization, resulting in the formation of yttria-stabilized zirconia (YSZ) and scandia-stabilized zirconia (ScSZ). By considering their formal oxidation states, they will include oxygen anion vacancies to be charge balanced. Indeed, these oxygen vacancies, and their interactions with other vacancies or cations, dramatically affect the structural, thermal, mechanical, and electrical properties of doped zirconia. Hence, these inclusions provide not only a way to gain stabilization in zirconia but also a means to direct ionic conduction where the ion transportion is directly attributed to the movement of these atomic defects (oxygen vacancies). 3–7 YSZ and ScSZ are extensively used as the electrolyte because of the high oxide ion and low electronic conductivity, as well as their stability under reducing and oxidizing atmospheres. 8–10 It is well known that the ionic conductivity of single crystal c-phase YSZ and ScSZ depends on the content of doped cations, exhibiting a maximum at around 8 mol % for YSZ 11–14 and 8-9 mol % for ScSZ. 15 The ionic conductivity of ScSZ (reported as 0.30 S/cm) is higher than that of YSZ (reported 0.13 S/cm) at high operating temperatures (typically 1000 ◦ C). This is most probably 2

ACS Paragon Plus Environment

Page 2 of 41

Page 3 of 41 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

due to the close similarity of the ionic radius of Sc3+ and Zr4+ (i.e. rSc3+ /rZr4+ = 1.03). 16 On the other hand, to improve long-term chemical and mechanical stability, and reducing operational cost, it is beneficial to lower the operating temperature of SOFCs from this technical point-of-view. In this aspect, the utilization of YSZ is not optimal at intermediate temperature. At intermediate temperatures, although ScSZ also exhibit high ionic conductivity (0.13 S/cm at 800 ◦ C), 17 there are some challenges such as high cost (due to its rarity and difficulties in separation), poor long-term stability, and polyphase features when compared to YSZ. Thus, ScSZ has also not been extensively explored in SOFC applications in spite of its potential. To improve the overall ionic conductivity and stability of electrolyte materials, researchers are seeking more sustainable alternatives such as heterogeneous doping and heterolayered structures with the following requirements satisfied: Long term stability, short term start-up response, and low cost in addition to high ionic conductivity. There is indeed considerable scientific and practical interest to study different compositions of YSZ-ScSZ for the improvement and enhancement of ionic conductivity and overall stability. For instance, to avoid the unwanted phase transitions and to further stabilize the c phase, the partial substitution of Sc3+ by Y3+ in ScSZ has been explored and studied. 18–23 The characteristics of the local environment of the oxygen ion is crucial for stability and ion transport, which may also affect the total number of mobile vacancy. 24,25 Of late, Xue et al. have suggested that interfaces formed in ScSZ electrolytes contribute to good cubic phase stability, low grain boundary resistance, and high ionic conductivity at intermediate temperature. 26 In this study, using an atomistic model, we examine the ion diffusion in Y/Sc-doped zirconia-based systems via performing molecular dynamics (MD) simulations with well-tested (i.e. determined based on the lattice parameters, lattice elastic properties, and high- and low-dielectric constants) interatomic potentials for YSZ and ScSZ. 27–29 We will focus on the oxygen ion conductivity of inhomogeneously- and heterogeneously-doped zirconia, and

3

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

Page 4 of 41

heterolayered doped zirconia, by considering five main mechanisms contributing to the enhancement of ionic diffusion at the interface: (1) the influence of type of cations with respect to the ionic size (2) the influence of cation concentration (3) the influence of cation distribution, (4) the influence of crystal structure, and (5) the nature of the interfaces.

2

Methodology

2.1

Interatomic potentials

The classical MD and lattice dynamics (LD) simulations are performed using the LAMMPS 30 and GULP code, 31 respectively. The effective potential of interatomic interaction between i and j ions is considered as a Buckingham-type with a long-range Coulombic term: ( ) qi qj rij Cij V (rij ) = + Aij exp − − 6 rij ρij rij

,

(1)

where qi holds for the charge of ion i (and likewise, qj for ion j), rij is the interionic distance, and other Aij , ρij and Cij parameters are used to represent the short-range terms for each element. The potential parameters for Zr−O, Y−O, Sc−O and O−O pairs are obtained from Reference 27 and Reference 29, which are given in Table 1. In the case of cation-cation interactions, forces are assumed to be purely Coulombic. The Coulomb interactions are summed up using Ewald’s method 32 with the formal charges of the ions. The analysis of the electronic structure of zirconia shows that zirconia is fairly ionic. 33,34 For this consideration, it is assumed that the considered zirconia-based systems are fully ionic, where the ionic charges are chosen to be 4+, 3+, 3+, and 2− for Zr, Y, Sc, and O, respectively. The time step used for the integration of Newton’s equations of motion is set to 0.5 femtoseconds (fs). The temperature and the pressure are controlled via the Nose-Hoover thermostat and barostat, respectively. 35–37 The Buckingham-type potential model has already

4

ACS Paragon Plus Environment

Page 5 of 41 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

Table 1: Interionic potential parameters used for the doped zirconia systems. i−j

Aij (eV)

ρij (˚ A)

O–O Zr–O Y–O Sc–O

9547.96 1502.11 1366.35 1575.85

0.224 0.345 0.348 0.3211

Cij (eV. ˚ A6 ) Reference 32.0 5.1 19.6 0.0

27 27 27 29

been successfully employed to describe the dynamic properties of zirconia-based materials and in previous literature. 11,38–40

2.2

Statistical averaging approach

In the construction of models, the statistical approach is assumed to better reconcile with experimental findings since the distribution of ions and defects in crystalline electrolytes are strongly inhomogeneous (or random). For this in mind, the (inhomogeneous) positions of dopant ions in zirconia are considered via a random number generation method. The distribution of doped ions in homogeneous (order) and inhomogeneous (random) systems is illustrated in Figure S1 of the Supporting Information. Hence, a series of randomly distributed doped zirconia-based structures are constructed instead of ordered doped zirconia systems. To average and minimize statistical fluctuations, five different random configurations (i.e. each one is generated as randomly distributions of doped ions) are set up for all models, and the calculations are then repeated and averaged. From this averaging procedure, we suggest that five randomized models for each considered system may provide a good estimate of its intrinsic diffusivity. Thus, to determine the statistical average, error bars are set for the five different random structures for each model. In the present work, each reported data is thus taken as the average over the populations of five different configurations.

5

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

2.3

Atomic geometries of bulk models

For the bulk systems, inhomogeneously-doped YSZ and ScSZ, and heterogeneously-doped YScSZ, the MD simulation cell is composed of a (9 × 9 × 9) supercell (taken with respect to the conventional fluorite crystal) with cell length of 46.04 ˚ A, and contains 8748 atoms with periodic boundary conditions. These bulk systems are studied within different dopant concentrations, i.e. 6, 8, 10, and 12 mol % for YSZ, 8 and 9 mol % for ScSZ, and 8 mol % for YScSZ. For the case of 8 mol % YScSZ, three different relative compositions of Y:Sc are modelled: (1) 2 mol % Y2 O3 and 6 mol % Sc2 O3 , (2) 4 mol % Y2 O3 and 4 mol % Sc2 O3 , and (3) 6 mol % Y2 O3 and 2 mol % Sc2 O3 . More details about the explicit number of ions for the considered systems are given in Table 2. For instance, for 8 mol % YSZ structure (denoted as YSZ), 432 of Zr4+ ions are randomly substituted by Y3+ doped ions, and 216 of O2− ions are randomly substituted by oxygen vacancies (VO ). Every two 3+ dopant ions in the system will require one VO to maintain the electrical neutrality of the system.

2.4

Atomic geometries of interface models

For the interface models, three-dimensional fluorite-based superlattice structures are considered. They compose of (9 × 9 × N ); where N represents the number of unit cells in the z-axis of the corresponding superlattice. Two types of interface models (namely, sharp and rough) are deliberated, and are based on the nature of interfaces constructed. For Y3+ ions that are abruptly separated from that Sc3+ ions at the interface region, we denote this interface model as the sharp interface model. For the case where the Y3+ and Sc3+ ions are inhomogeneously distributed within the same region near the interface, such an interface model is denoted as the rough interface model. The number of dopant atoms and the concentration ratio for both superlattice models are taken same. Both sharp and rough interface models are constructed from [001], [110] 6

ACS Paragon Plus Environment

Page 6 of 41

Page 7 of 41 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

and [111] of the bulk cubic fluorite crystal, and hence denoted as (100)s , (110)s and (111)s , and (100)r , (110)r and (111)r respectively (see Table 2). For instance, the (100)s superlattice, which is constructed from the [001] crystal direction when considering a (9 × 9 × 18) supercell, contains 4968 Zr4+ , 432 Y3+ , 432 Sc3+ , 11232 O2− ions, and 432 VO . The regions which are YSZ-dominant, interface and ScSZ-dominant within the (100)s are then further denoted as (100)s YSZ, (100)s I and (100)s ScSZ, respectively. To simulate superlattice models, the elementary domain is divided into five zones for both sharp and rough interface models as illustrated in Figure S2 of the Supporting Information. The concentration of each zone is taken the same as 8 mol %, i.e. for (100)s sharp model, two half layers of (100)s YSZ surround one (100)s ScSZ layer, in which the doped concentration of (100)s YSZ, (100)s I, and (100)s ScSZ is 8 mol %. For each superlattice, periodic boundary conditions are also applied in all three directions.

2.5

Radial distribution function and coordination number

The radial distribution functions (RDF) for these zirconia-based systems are derived to characterize their atomic geometries. The averaged number of target atoms (α) is counted when they are positioned within a distance range of r to r + δr, centering for each given α atom. In this approach, δr is taken as a small perturbation to the distance, r. In other words, the RDF (gαβ (r)) is a mathematically rigorous approach to describe the density probability for an atom of the α species to have a neighboring β species atom at a given distance r. This can be elegantly expressed as,

gαβ (r) =

dnαβ (r) Nα 4πr2 dr V

.

(2)

Here, V represents as the volume of the periodic simulation cell and Nα is equal to cα N , where cα is the concentration of species α and N the total number of atoms in the cell.

7

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

Page 8 of 41

Table 2: Number of ions and vacancies in different bulk and superlattice(with interface regions) models for doped zirconia-based systems. Here, M2 O3 (in mol %) is the percentage mole fraction of dopant (where M denotes either Y3+ or/and Sc3+ ) in the host ZrO2 matrix, and dz (in ˚ A) is the length of considered simulation cell in z-axis, defined as the thickness for superlattice illustrated in Figure S2 of the Supporting Information.

Model

M2 O3 (mol %)

Zr4+

Y3+

Sc3+

O2−

VO

dz (˚ A)

6.04 8.00 8.00 9.01 10.04 12.07 8.00 8.00 8.00

2584 2484 2484 2434 2384 2288 2484 2484 2484

332 432 0 0 532 628 108 216 324

0 0 432 482 0 0 324 216 108

5666 5616 5616 5591 5566 5518 5616 5616 5616

166 216 216 241 266 314 216 216 216

46.04 46.04 46.04 46.04 46.04 46.04 46.04 46.04 46.04

8.00 8.00 8.00 8.00 8.00 8.00

1104 1656 2760 4948 4416 3680

96 144 240 432 384 320

96 144 240 432 384 320

2496 3744 6240 11232 9984 8320

96 144 240 432 384 320

20.46 30.70 51.16 92.09 86.82 88.61

8.00 8.00 8.07

552 552 490

48 48 43

48 48 43

1248 1248 1109

48 48 43

10.27 10.85 11.82

Bulk 6YSZ YSZ ScSZ 9ScSZ 10YSZ 12YSZ YScSZ 4Y4ScSZ 6Y2ScSZ Superlattice (100)s d20 (100)s d30 (100)s d50 (100)s (110)s (111)s Interface region (100)s I (110)s I (111)s I

The coordination number (CN) for these zirconia-based systems is defined as the number of neighbor atoms with specified atom type which are within the specified cutoff distance from the central atom. Here, the cutoff distance is determined by the RDF analysis.

8

ACS Paragon Plus Environment

Page 9 of 41 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

2.6

Quenching process

Ionic diffusion analysis is performed using the thermally-quenched structures. More specifically, all structures are first exposed to total energy minimization and then the quenching process. Firstly, each initial configurations has been optimized by the conjugate gradient energy minimization, which relaxes the position of ions to a lower energy configuration. After the energy minimization process, a quenching process is set to accelerate the equilibration process for each system under a temperature gradient. This approach directs the positions of vacancies and ions to the global minimum in avoidance of artificial defect formation. The quenching process is performed as follows: The system is firstly heated to 2000 K using a constant pressure and temperature ensemble (NPT) to ensure thermal equilibrium is reached for the desired temperature. After that, the constant volume and energy ensemble (NVE) is enforced and the NVE MD simulation runs for 500 ps. Subsequently, the system is then rapidly cooled down (hence, quenched) from 2000 to 300 K with the NPT ensemble for 50 ps, and then further relaxed for another 200 ps with the NVE ensemble. Finally, the system is equilibrated with NPT and NVE ensembles at 300 K temperature. Figures S3, S4, and S5 indicate that the considered systems preserve their crystalline structure and cubic phase stability during each step of the quenching process. With this quenching process, the quenched atomic structures for a given target temperature are then used for further analysis.

2.7

Mean square displacement and diffusion coefficient

The calculated ion diffusion coefficient (D) for all zirconia-based models is computed for the quenched structures via the consideration of the mean square displacement (MSD), < R2 >. The < R2 > is calculated and collected throughout each MD run as a function of the simulation time (t). The D for each system is subsequently determined according to the Einstein relation: 41 9

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

< R2 >=

N 1 ∑ < [Ri (t0 + t) − Ri (t0 )]2 > = 6D(t)t + B N i=1

where B is a constant, the brackets and



Page 10 of 41

,

(3)

represent taking an average over all simulation

time, and the summation over all ions, respectively. The diffusion coefficient, D can be related to the conductivity, σ and it is then calculated according to the Nernst-Einstein relation given in Equation 4:

σ=−

Z 2ρ D kB T HR

,

(4)

where Z and ρ are the ionic charge and the number density of the conducting ions, respectively, and HR the Haven ratio. For cubic fluorite structures with a low vacancy/defect concentration, 42 HR may be taken as 0.65 and this value is commonly reported in literature for stabilized zirconia. 11,43,44 Next, the activation energy barrier, Ea , of oxygen ion diffusion are determined using the Arrhenius equation, 45 D = D0 exp(−

Ea ) , kB T

(5)

where T corresponds to the target temperature, kB the Boltzmann constant, and D0 for the temperature-independent pre-factor. This can be linearlized to the following:

log D = −

Ea 1 + log D0 kB T

10

ACS Paragon Plus Environment

.

(6)

Page 11 of 41 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

3

Results and discussion

3.1

Inhomogeneously- and heterogeneously-doped zirconia

3.1.1

Structural properties and analysis

Before discussing the doped zirconia systems, it is necessary to mention about the phase property of the pristine ZrO2 and M2 O3 (where M = Y or Sc). It is known that the pristine zirconia has three different crystal polymorphic expressions, depending on its temperature and pressure. From experiments, at room temperature, the most stable phase of zirconia is the monoclinic structure (m-phase), and stays stable over a very large temperature range of up to 1170 K; under an ambient pressure and high temperatures, t- and c-phases (at 1478 and above 2650 K, respectively) become thermodynamically stable. 46 However, at room temperature, the effective stabilization of the c and t phases can be achieved with aliovalent doping. Likewise, M2 O3 , with M = Y or Sc, crystallizes in three different lattice structures – namely, the m-, h- (hexagonal), and c- phases. As can be seen from Figures 1(a) and 1(b), the energetically most stable phases (at T = 0 K) for ZrO2 and M2 O3 are the m-phase and c-phase, respectively. Given that our discussion will be centered around the technical operating conditions of SOFC (i.e. ambient pressures and high temperatures), the cubic phases of both ZrO2 and M2 O3 (as illustrated in Figures 1 (c) and 1 (d)) will be used to build the atomic models in this study. The optimization of the lattice parameters of the zirconia-based structures is first performed at both T = 0 K and T ̸= 0 K using LD and MD simulations, respectively. Firstly, the lattice constants (a, b, c , α, β, γ), volume (V ), density (d), and energy (E) of the systems are studied at T = 0 K using LD simulations. It is found that the lattice constant a of the ˚ 48 pure c-ZrO2 is 5.08 ˚ A, which agrees well with the experimental values of 5.08 ˚ A , 47 5.09 A (at T = 300 − 600 K) and theoretical value of 5.11 ˚ A at T = 0 K. 49 As listed in the Table 3, the corresponding structural parameters for the 5 randomized 11

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

Figure 1: (a) The total energy (E) versus volume (V ) relations for (a) ZrO2 with three different crystal structures, namely the monoclinic (m), tetragonal (t), and cubic (c) phases, and (b) M2 O3 , with M = Y or Sc, with three different crystal structures, namely, the monoclinic (m), hexagonal (h), and cubic (c) phases. The m, t (h), and c phases are represented by orange, blue, and red colors, respectively. Each marker on the graphs represents a calculated data point, while the lines show the best fit. (c) The unit cell of c-ZrO2 , with the three different crystal directions (i.e. [100], [110], and [111]) represented by the yellow, orange, and blue planes, respectively. (d) The unit cell of c-M2 O3 , with M = Y or Sc. Here, the Zr, O, dopant (Y or Sc) atoms are shown in blue, red, and orange, respectively. models of YSZ, ScSZ and YScSZ are provided. The lattice parameter a for YSZ, ScSZ and YScSZ is found to be in the range of 5.11 to 5.12, 5.08 to 5.09, and 5.09 to 5.10 ˚ A, respectively, which indicates the influence of the ionic size of the doped cations on the crystal lattice. Kim et al, have reported that there is indeed a correlation between the lattice parameter and the ionic size of the dopant cations for fluorite-structured oxides. 50 Given that the ionic radius of Sc3+ is very similar to that of Zr4+ (rSc3+ / rZr4+ = 1.03), their a values are comparable. We 12

ACS Paragon Plus Environment

Page 12 of 41

Page 13 of 41 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

Table 3: Structural parameters of inhomogeneously-doped YSZ, ScSZ, and heterogeneouslydoped YScSZ bulk systems with 8 mol % dopant concentration. Here, the superscript label (from 1 to 5) represents a series of 5 randomized models for each system. a (˚ A)

b (˚ A)

c (˚ A)

α (◦ )

β (◦ )

γ (◦ )

V (˚ A3 )

d E 3 (g/cm ) (eV/atom)

c–ZrO2 YSZ1 YSZ2 YSZ3 YSZ4 YSZ5

5.082 5.114 5.113 5.114 5.117 5.114

5.082 5.114 5.113 5.121 5.112 5.116

5.082 5.110 5.114 5.117 5.116 5.116

90.0 89.9 90.1 90.2 90.2 90.0

90.0 90.0 90.0 90.0 89.9 90.0

90.0 89.9 90.0 89.9 90.0 89.9

131.289 133.637 133.710 133.978 133.834 133.851

6.234 6.049 6.045 6.033 6.040 6.039

-37.334 -36.061 -36.062 -36.062 -36.061 -36.065

ScSZ1 ScSZ2 ScSZ3 ScSZ4 ScSZ5

5.078 5.082 5.081 5.081 5.080

5.080 5.082 5.083 5.078 5.081

5.080 5.085 5.086 5.083 5.081

90.0 90.0 89.9 90.2 89.8

90.0 90.0 90.0 89.9 90.0

89.9 90.1 90.0 90.1 89.9

131.062 131.314 131.358 131.123 131.134

5.837 5.826 5.824 5.835 5.834

-36.256 -36.250 -36.251 -36.250 -36.252

YScSZ1 YScSZ2 YScSZ3 YScSZ4 YScSZ5

5.094 5.099 5.104 5.099 5.099

5.096 5.097 5.099 5.101 5.095

5.096 5.097 5.102 5.104 5.098

90.0 89.9 90.0 90.1 89.9

89.8 90.1 90.1 90.1 89.9

90.1 90.0 90.0 89.8 89.9

132.297 132.477 132.792 132.761 132.447

5.946 5.938 5.924 5.926 5.940

-36.160 -36.157 -36.152 -36.151 -36.153

also notice that the a and E values for the randomized models do not change dramatically for each system. Using these structural parameters (as determined from the LD simulations at T = 0 K), we next proceed to perform MD simulations at different elevated temperatures. The temperature for each MD run is raised from room temperature (300 K) to 1800 K under constant pressure within a NPT ensemble (which accommodates volumetric changes). The a lattice parameter increases linearly with temperature, indicating expansion of the lattice cell (as shown in Figure 2). We note that higher fluctuations in the a lattice parameter expansion (with increasing temperature) can be seen for the heterogeneously-doped YScSZ as the larger difference between the radii of Y3+ and Sc3+ leads to a more dramatic lattice distortion, which may reflect poorer structural integrity at the operating temperature for SOFC.

13

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

Figure 2: The a lattice parameter (in ˚ A) of YSZ, ScSZ, and YScSZ with 8 mol % dopant concentration as a function of temperature T (in K). These zirconia-based systems are heated from 300 to 1800 K with an increment of 100 K. The pressure is set to 1 atm within the NPT ensemble, which accommodates for volumetric changes in the system. The red, yellow, and blue markers represent the data points for YSZ, YScSZ, and ScSZ, respectively. The distribution of the dopant ions (dark red and dark blue spheres for Y3+ and Sc3+ , respectively) in the randomized models for each system is depicted in the atomic structures shown on the right. Next, to provide an insight into long-range (dis)order analysis in these zirconia-based systems, the pair radial distribution function (RDF), ρ(r) for each system (at elevated temperatures) is calculated and plotted in Figure 3 and Figure S6 of the Supporting Information. The RDF of c-ZrO2 at T = 0 K (in Figure S6 (a)) is provided as a reference. Referring to the RDF of c-ZrO2 at T = 1000 K (in Figure S6 (d)), we find that the ρ(r)Zr−O (shown in black lines) has high-intensity peaks at r = 2.202 and 4.218 ˚ A, while the ρ(r)O−O (in orange) exhibits the highest peak at 2.538 ˚ A and followed by much lower peaks of intensity after 3.594 ˚ A. These peaks can be explained by the nearest and the next-nearest neighbor distances in crystal lattice, with a ideal lattice constant, a. In the ideal cubic fluorite struc√ √ ture, these distances can be defined as a 3/4 and a 11/4 for the cation-oxygen distances, √ whereas a/2 and a/ 2 can be assumed for that of oxygen-oxygen distances, respectively. Likewise, in the case of the inhomogeneously- and heterogeneously-doped zirconia systems, we examine the ρ(r) (at T = 1273 K) for the constituent ions in YSZ (Figure 3 (a)), 14

ACS Paragon Plus Environment

Page 14 of 41

Page 15 of 41 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

Figure 3: The radial distribution function (RDF), ρ(r) and coordination number (CN) of doped zirconia systems: (a,d) YSZ, (b,e) ScSZ, and (c,f) YScSZ at T = 1273 K. The ρ(r) for the Zr-O, O-O, Y-O, and Sc-O bond pairs are represented by the black, orange, red, and blue, respectively. Here, the first intense peak for each pair is also shown in the inserts. The CN of oxygen to Zr4+ , Y3+ and Sc3+ ions as a function of time is represented by the black, red, and blue lines, respectively. ScSZ (Figure 3 (b)), and YScSZ (Figure 3 (c)). The ρ(r)Zr−O for YSZ, ScSZ, and YScSZ are found to have the highest peak at r = 2.094, 2.105, and 2.118 ˚ A, respectively, while the ρ(r)O−O peaks at r = 2.610, 2.575, and 2.586 ˚ A for the corresponding doped systems. The most intense peak for Y(Sc)-O bond pair RDF is located at r = 2.310 ˚ A (2.286 ˚ A), while that for the Y-O and Sc-O bond pair RDF of YScSZ are 2.185 and 2.154 ˚ A, respectively. In agreement with literature, 51 the first highest peak in ρ(r)Y−O is at a much larger distance than that in ρ(r)Zr−O and ρ(r)Sc−O . Due to the difference in ionic radii of Y3+ and Sc3+ , the relative difference in the maximum peak for the Zr-O bond pair is indeed smaller than that of the Y-O pair (0.17 ˚ A) and the Sc-O pair (0.13 ˚ A). This is clearly reflected in the RDF of ScSZ (Figure 3(b)) where the main peaks are generally left-shifted when compared to the

15

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

other systems. Thus, the RDFs of anion-cation pairs exhibit two, (almost) one, and three main peaks for YSZ, ScSZ and YScSZ systems, respectively, due to the ionic size difference between host and doped cations. In the light that only one type of cation-anion pair in cubic undoped ZrO2 , the existence of three different anion-cations pairs for the YScSZ may well account for the higher structural distortions and thus lower structural integrity in YScSZ. Furthermore, we do find that the RDF of cation-anion and anion-anion pairs generally indicates some loss of long-range order due to a mobile oxygen sub-lattice. To illustrate this, we have calculated the RDF for heterogeneously-doped zirconia system at a higher temperature of 2000 K in Figure S3 (f). Here, we observe that the vibrational amplitude of ions increases with increasing temperature (i.e. the width of the peaks), but the locations of the RDF peaks are essentially temperature-independent. The O-O bond distance is much more affected by the increase in temperature, and some re-normalized minor peaks are observed around the more intense peaks due to the existence of oxygen vacancies. As expected, all peaks in the ρ(r) are broadened as the temperature is increased. The broadening of subsequent peaks at larger bond distances suggest greater disorder and increased oxygen diffusion at higher temperatures. In addition to the RDF, the coordination number, CN has been computed for different pairs of ions within the specified cutoff distance, which is determined by the first intense peak found in the RDF. We examine the CN (at T = 1273 K) for the cations in YSZ (Figure 3 (d)), ScSZ (Figure 3 (e)), and YScSZ (Figure 3 (f)). The CN of Zr4+ ions in the ideal fluorite oxide structure is 8. The calculated CN of Zr4+ for all the doped systems is found to be lower than 8, which indicates the presence of vacancies. The distribution and preferential sites of oxygen vacancies are also investigated under finite temperature (at T = 1273 K) for 750 ps. We note that the average CN of Zr+4 ions is lower than that of Y3+ for YSZ system, which confirms that the oxygen vacancies prefer to bind closer to the Zr4+ cations, agreeing with previous literature. 14 More generally, our results indicate that the oxygen vacancies prefer to bind to the smaller

16

ACS Paragon Plus Environment

Page 16 of 41

Page 17 of 41 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

cations in doped zirconia system. For the ScSZ system, the CN of Zr4+ is very close to that of Sc3+ . An increase in the number of preferred vacancy sites is typically associated with an increase in the number of mobile vacancies – i.e. enhancing the overall structural stability and improving the oxygen ion mobility of the system. Here, the vacancies do not seem to have a clear preference, which corroborates well with the observed increased ion mobility (with minimized lattice strain). This is also very much inline with reported experimental and theoretical results 28,52,53 where the minimum binding energy for cation-oxygen vacancy is indeed for the Sc4+ cation (when compared to cations of In, Yb, Y, Gd, Sm, Ce and La). Thus, we can conclude that the total number of mobile vacancies is higher than that of YSZ system, which improves its ionic mobility. For the YScSZ system, the CN of Zr+4 and Sc3+ is found to be very similarly to that in the ScSZ system, with the CN of Y3+ increasing slightly. It is worth noticing that from the CN analysis, the total number mobile O vacancies for YScSZ is lower than YSZ and ScSZ, which suggests dominant vacancy-vacancy and vacancy-defects interactions. We argue that the immobile oxygen vacancies in heterogeneously-doped YScSZ may tend to localize due to these deleterious interactions under technical conditions which can lead to structural instability of heterogeneous doped systems. It is noted that the the fluctuation of CN for heterogeneous-doped YScSZ is higher than that of inhomogeneous-doped YSZ and ScSZ at 1273 K for 750 ps. These fluctuations are associated with larger differences in cationic charges and radii, which leads to the discussed excessive lattice strain. Thus, the RDF, CN, and binding energy (vacancy-cation) analysis clearly reveal that the radii of the dopant cations do indeed play an important role in the distribution of oxygen vacancies, crystal lattice strain, and ion mobility in doped zirconia. To account for the effect of configurational entropy accurately for these large heterogeneouslydoped YScSZ and heterolayered YSZ/ScSZ systems is non-trivial. For instance, for the case of 8 mol % YScSZ in a 3×3×3 supercell, it will contain 92 Zr ions, 16 doped-metal ions, and 208 O ions, and the total of configurations in this system will easily reach 5.55×1042 (when

17

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

doped ion is only Y or Sc) which becomes intractable very quickly. However, to estimate the possible effect of configurational entropy in our systems, we use a smaller model system of YSZ where 4 Y ions and 2 O vacancies are considered in order to reduce the number of configurations. Using the site occupancy disorder program, 54 we construct 100,000 different configurations of Y and O vacancy positions. The resulting energies of these 100,000 configurations (as a function of different Y ion distribution in the lattice) are shown in Figure S7 in the Supporting Information where each vertical column depicts almost 2,000 configurations of O vacancies (for a fixed Y ion distribution). Thus, given that the contribution of configurational entropy from this smaller model already reflects a rather small contribution, and agreeing with a recent study that claims the configurational entropy for 8 mol % YSZ at T = 1, 000 K is less than 0.03 eV, 55 we argue that its effect on the overall stability maybe negligible and will not change the conclusions drawn from this work.

3.1.2

Ionic diffusion and conductivity

The distribution of ions and defects (e.g. vacancies) in crystalline electrolytes is usually strongly inhomogeneous. To examine the influence of the type (doped ion size) and distribution of doped ions on the mobility of oxygen ions and vacancies in these zirconia-based electrolytes, inhomogeneously-doped YSZ and ScSZ, as well as the heterogeneously-doped YScSZ bulk systems are investigated using MD simulations. Firstly, the mean square displacement, < R2 >, of oxygen ions at different dopant concentrations (i.e. 6, 8, 10, and 12 mol % for YSZ, 8 and 9 mol % for ScSZ, and 8 mol % for YScSZ) are calculated at T = 1800 K for 500 ps (Figure 4). The < R2 > values increase with time for all given concentrations, reflecting the diffusion of mobile oxygen ions in our MD simulations. Most importantly, as illustrated in Figure 4 (a), the slope of the < R2 > as a function of time shows that the highest mobility is exhibited at 8 mol % dopant concentration for both YSZ and ScSZ. This is in good agreement with earlier reports in literature where a maximum ionic conductivity is observed for 8 mol % YSZ and 8 to 9 mol % ScSZ. 11–15

18

ACS Paragon Plus Environment

Page 18 of 41

Page 19 of 41 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

˚2 ) of oxygen ions Figure 4: The mean square displacement, < R2 > (cf. Equation 3 in A as a function of time (t, in ps) YSZ, ScSZ, and YScSZ at T = 2000 K. (a) 6, 8, 10, and 12 mol % YSZ (denoted as 6YSZ, 8YSZ, 10YSZ, and 12YSZ, respectively) and 8 and 9 mol % ScSZ (denoted as 8ScSZ and 9YSZ, respectively). Their corresponding data points are represented by green, red, brown, pink, blue, and sky blue, respectively. (b) 8 mol % YScSZ with different Y/Sc ratio. Here YScSZ, 4Y4ScSZ, and 6Y2ScSZ are denoted in blue, orange, and red, respectively. Indeed, naively, one may expect that increasing the number of oxygen vacancies by introducing more dopants may simply increase oxygen ion mobility. However, defect-defect interactions (namely, the vacancy-vacancy and cation-vacancy interactions) will play an important role in the mobility of ions. From this point of view, in the case of YSZ, vacancycation interaction is predicted to be more dominant since the ionic radius of Y3+ is larger than Zr4+ . As mentioned above, the oxygen vacancies in YSZ prefer to located at the nearest neighbor site of the smaller cations (Zr4+ ). Due to this preference of the position of oxygen vacancies, increasing of the dopants (beyond a certain threshold – in this case, 8 mol %) will induce a stronger attraction between the oxygen vacancies and the cations, reducing the total number of mobile vacancies, and thus reducing their overall mobility. On the other hand, in the case of ScSZ, cation-vacancy interaction is not as dominant due to the similarity in cationic sizes between Sc3+ and Zr4+ (cf. rZr = 0.79 ˚ A, rSc = 0.81 ˚ A). 16

19

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

Figure 5: Mean square displacements (< R2 >(t)) of oxygen ions in (a) 8 mol % YSZ, (b) 8 mol % ScSZ, and (c) 8 mol % YScSZ (with 2 mol % Y2 O3 and 6 mol % Sc2 O3 ). Here, the MD runs are conducted at various finite temperatures, ranging from T = 700 to 2,000 K, with an increment of 100 K. The < R2 >(t) values as a function of T for YSZ, ScSZ, and YScSZ are monochromatically scaled in red, blue, and orange, respectively. This may not dramatically impede the migration of oxygen ions. However, given that ScSZ is sensitive to the polymorphic phase transition with increasing dopant concentration, oxygen trapping at certain lattice sites may occur which is commonly found in zirconia-based ceramics. 56 The < R2 > is also calculated for 8 mol % YScSZ with a varying Y/Sc ratio (see Figure 4 (b)). It is found that the ionic diffusion increases with the lowering of Y3+ dopants due to the increased amount of free space in the crystal, as supported by previous theoretical 57 and experimental 22 studies, upholding our argument that a more severe lattice distortion for high Y3+ concentrations is expected to impede oxygen ion mobility. By simply comparing the < R2 >(t) slope computed for all the dopant concentrations of YSZ, ScSZ, and YScSZ discussed here, it is clear that expected order of the highest oxygen ion mobility should be ScSZ, followed by YSZ, and then the hetero YScSZ, where the dopant concentration of 8 mol % for each system exhibits the highest oxygen ion mobility. Now, focusing our attention only on the dopant concentration of 8 mol % for all zirconiabased systems in this work, we will simply refer to them as YSZ, ScSZ and YScSZ accordingly

20

ACS Paragon Plus Environment

Page 20 of 41

Page 21 of 41 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

Table 4: The ionic conductivity, σ (cf. Equation 4, in Scm−1 ) of YSZ, ScSZ, and YScSZ for temperatures ranging from 1000 to 1800 K. Note that for each system, five randomized structural models have been used for statistical averaging and the statistical error bar (∆σ) is also reported.

T (K) 1000 1200 1400 1600 1800

YSZ

ScSZ

YScSZ

σ±∆σ (Scm−1 ) 0.015 ± 0.001 0.041 ± 0.001 0.084 ± 0.003 0.136 ± 0.004 0.198 ± 0.005

σ±∆σ (Scm−1 ) 0.032 ± 0.002 0.076 ± 0.006 0.143 ± 0.002 0.218 ± 0.005 0.294 ± 0.004

σ±∆σ (Scm−1 ) 0.013 ± 0.001 0.036 ± 0.002 0.069 ± 0.002 0.110 ± 0.003 0.156 ± 0.003

hereafter for the ease of discussion/notation. In Figure 5, for range of temperatures (i.e. from 700 to 2000 K with an increment of 100 K), the < R2 >(t) for each system is evaluated at each temperature (for 500 ps) and is found to be varying linearly. From the linear slope of these < R2 >(t), the diffusion coefficients (D) for these systems (which are given in Table S1 of Supporting Information) are then determined via Equation 3 and the ionic conductivity, σ via Equation 4. In the perspective of ionic conductivity, it is worth mentioning that σ at low temperatures is found relatively small compared to that at high temperatures. On the average, the value of σ for YSZ and ScSZ is calculated as 0.1 and 0.2 Scm −1 , respectively, which are in good agreement with the corresponding reported experimental values of 0.13 Scm−1 for YSZ and 0.24 Scm−1 for ScSZ at technical conditions (T = 1000 ◦ C). 17 Details of the calculated σ values for other zirconia-based systems can be found in Table 4 where the statistical error bars for the 5 randomized models for each system are provided as well. In comparison, for the hetero-YScSZ system, σ is calculated to range between 0.013 to 0.156 Scm−1 , which is lower than that of YSZ and ScSZ when considering T = 1000 to 1800 K. Indeed, one would expect that the σ of YScSZ to be of an intermediate value between that of YSZ and ScSZ. The heterogeneously-doped YScSZ system includes three different types of

21

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

cations Zr4+ , Y3+ , and Sc3+ with a ionic size difference of 0.79, 0.96, and 0.81 ˚ A, respectively (as illustrated in Figure 3 (a)). Here, we argue that the existence of the different character of the anion-cations pairs (i.e. their charges, ion sizes, and bond distances) may lead to breaking of the pristine cubic symmetry and the loss of octahedral interstials for oxygen ion transportation. The excess defects and lattice distortions may then cause a larger barrier and hindrance to ionic transport of oxygen. 26,58 This can be inferred from the fluctuations in the lattice parameters (at T = 300 K and higher temperatures), measured energies, and the CN. Here, we rationalize the decrease in σ for YScSZ is primarily due to the excessive lattice distortions that hinders ionic diffusion 58 and the decrease in the total number of mobile vacancies. As already seen in Figure 2 and Figure 3, the larger a lattice fluctuations, the peak separations in RDF, and the CN fluctuations of YScSZ (as compared to that of YSZ and ScSZ) support this argument. To further examine the ionic diffusion for cations, we calculate the < R2 > for Zr4+ , Y3+ , and Sc3+ , which are found to be almost close to zero and do not vary with time, even at T = 1800 K, (as illustrated in Figures 6(a), 6(b), and 6(c)). The cation transport at operating temperatures have been reported as more than one-trillion-times slower than that for the anion in the zirconia based electrolytes. 59,60 Hence, it is reasonable to conclude that these metal cations are not expected to diffuse (even for very high temperatures) and will only exhibit vibrational features on their lattice sites. Specifically, for the YSZ, ScSZ and YScSZ systems, very small displacement and oscillation amplitude of cations are indeed observed in our MD simulations. On this note, we postulate that a higher cation oscillation in the system may then explain its lower stability. We note that the oscillation amplitude of cations in YScSZ is higher than that of YSZ and ScSZ, which then reflects its reported instability. In addition, in Figure 6 (d), by resolving the < R2 >(t) into the 3 orthogonal directions of the simulation cell – namely, the x-, y-, and z-directions, we examine the direction-dependence ionic diffusion at T = 1800 K and

22

ACS Paragon Plus Environment

Page 22 of 41

Page 23 of 41 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

Figure 6: < R2 >(t) for Zr4+ , Y3+ , and Sc3+ cations (shown in black, red, and blue, respectively ) for (a) YSZ, (b) ScSZ, and (c) YScSZ as a function of time at T = 1800 K. (d) < R2 >(t) for oxygen ions (resolved into the orthogonal x-direction (Ox ), y-direction (Oy ) and z-direction (Oz )) for YSZ (in gradated red lines), ScSZ (in gradated blue lines), and YScSZ (in gradated orange lines) at T = 1800 K. (e) The Arrhenius plots of the diffusion coefficient (log D) as a function of 1000/T for YSZ, ScSZ, and YScSZ where for each system, five randomized structural models have been used for statistical averaging. The red, blue, and orange lines are obtained from a linear regression fit for each data set of YSZ, ScSZ, and YScSZ respectively. From this, the calculated diffusion activation energy barrier, Ea (cf. Equation 6) for YSZ, ScSZ, and YScSZ is 0.62 ± 0.01, 0.55 ± 0.02, and 0.60 ± 0.01 eV, respectively. find that the oxygen ion diffusion in YSC, ScSZ, and YScSZ is rather isotropic about the orthogonal < 100 > directions considered. It has been previously reported that the ionic diffusion coefficient, D is proportional with temperature, T while the diffusion activation energy barrier, Ea may somewhat vary for different temperature regimes. 61,62 Especially at low temperatures, the calculated Ea

23

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

may fluctuate with larger statistical errors and these errors are greatly reduced at higher temperatures. This is mainly attributed to the dominant vacancy-dopant and vacancy-ion interactions at low temperatures, and with increasing higher temperatures this interaction energy barrier may be overcome to stabilize the mobility of ions. Thus, to examine this temperature-dependent diffusion phenomena, we proceed to calculate Ea (cf. Equation 6) for YSZ, ScSZ, and YScSZ. Using the Arrhenius relation, we plot the diffusion coefficient (log D) as a function of 1000/T in Figure 6(e) to determine Ea . Here, we note that for each system, five randomized structural models have been used for statistical averaging and the statistical error bars (as shown in small vertical bars) can be used to probe the dominance of vacancy-dopant and vacancy-ion interactions in each case. A larger error bar in D will indicate higher fluctuations in the calculated Ea values. The calculated Ea for YSZ, ScSZ, and YScSZ is 0.62 ± 0.01, 0.55 ± 0.02, and 0.60 ± 0.01 eV, respectively, and agreeing well with the experimentally determined values for YSZ and ScSZ. 63–65 Our results show that the fluctuations in D seen at higher temperatures are indeed significantly reduced as compared to that found for lower temperatures. Interestingly, when comparing these fluctuation trends in the statistical errors of D, the heterogeneously-doped YScSZ exhibits the smallest error bars even at the low temperature regime, suggesting that some entropic disorderedness in the metal cation distribution might assist in the overcoming of the vacancy-dopant and vacancy-ion interactions at low temperatures. To understand how the randomized heterogeneity of YSZ and ScSZ in YScSZ can affect its ionic diffusion and conductivity, we proceed to explore these transport phenomena in heterolayered YSZ and ScSZ interfaces.

3.2

Heterolayered YSZ/ScSZ interfaces

The diffusion transport properties across the YSZ and ScSZ heterolayered interface (denoted as YSZ/ScSZ) have been evaluated at various temperatures. All possible heterolayered structures (e.g. with different thicknesses of YSZ and ScSZ regions, different directions of 24

ACS Paragon Plus Environment

Page 24 of 41

Page 25 of 41 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

Figure 7: Mean square displacements (< R2 >(t)) of oxygen ions for different YSZ/ScSZ heterolayered models: (100)s d20 (shown in black), (100)s d30 (shown in red), (100)s d50 (shown in blue), and (100)s (shown in brown) where the corresponding superlattice thickness, dz = 20.46, 30.70, 51.16, and 92.09 ˚ A, respectively, for 750 ps at T = 1800 K. For each (100)s YSZ/ScSZ model, the calculated D value (in m2 s−1 ×10−12 ) is also provided in parentheses. YSZ/ScSZ superlattice structures with different thickness (as defined by dz ) are illustrated on the right. the heterointerface formed, and different interface roughness) are considered and studied systematically using a superlattice model. We first analyze the relationship between the atomic structure and ionic diffusion by considering three-dimensional superlattice models constructed from the [100], [110], and [111] directions of the cubic fluorite crystal lattice. Next, the ionic conductivity of these YSZ/ScSZ models is calculated as a function of space across the interface, and specifically the atomic roughness of the interface is described either via a sharp and rough interface profile. Four different YSZ/ScSZ models (in the [100] direction) of different superlattice thickness, dz are considered: (100)s d20 , (100)s d30 , (100)s d50 , and (100)s where dz = 20.46, 30.70, 51.16, and 92.09 ˚ A, respectively. The total dopant concentration is kept the same at 8 mol %. The < R2 >(t) has been calculated from MD runs for 750 ps at T = 1800 K and the D values as a function of dz are then determined as shown in Figure 7.

25

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

Figure 8: (a) Region-dependent < R2 >(t) of oxygen ions for YSZ/ScSZ (100)s : The dominantly YSZ region (in red and labelled as (100)s YSZ), the interface region (in orange and labelled as (100)s I), and the dominantly ScSZ region (in blue and labelled as (100)s ScSZ), respectively, for 750 ps at T = 1200 K. (b) The Arrhenius plots of the diffusion coefficient (log D) as a function of 1000/T for (100)s YSZ, (100)s I, and (100)s ScSZ where five randomized structural models have been used for statistical averaging. The red, orange, and blue lines are obtained from a linear regression fit for each data set of (100)s YSZ, (100)s I, and (100)s ScSZ, and their calculated Ea values are provided in the parentheses. We find that the D value increases with increasing values of dz . The lowest value seen in (100)s d20 may be due to an interface-interface interaction at too small a dz value, while (100)s (with dz = 92.09 ˚ A) may be taken as a good converged reference point to discuss the region-dependent transport properties across the YSZ/ScSZ interface. On the other hand, the < R2 >(t) of cations in YSZ/ScSZ is very close to zero, as indicated in Figure S8 of the Supporting Information, with a small oscillation amplitude reflecting its stability. To further understand the region-dependence of diffusion in the (100)s model, we calculate the < R2 >(t) for oxygen ions in the dominantly YSZ region (labelled as (100)s YSZ), the interface region (labelled as (100)s I), and the dominantly ScSZ region (labelled as (100)s ScSZ), as shown in Figure 8 (a) and illustrated in Figure S2 (left) of the Supporting Information. Our results show that the mobility of oxygen ions is the highest in the ScSZ-dominated region ((100)s ScSZ) while the lowest in (100)s YSZ. This is clearly reflected in the Arrhenius plots (Figure 8(b)) where the calculated Ea for (100)s YSZ, (100)s I, and (100)s ScSZ is 0.61, 26

ACS Paragon Plus Environment

Page 26 of 41

Page 27 of 41 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

Figure 9: The radial distribution function (RDF), ρ(r) of the (a) Zr-O, (b) Y(Sc)-O, and (c) O-O. Here, (100)s YSZ, (100)s I, and (100)s ScSZ are illustrated by red, yellow(black), and blue, respectively, at T = 1273 K. The first intense peaks of each pair are also shown in the inserts. (d) The model for heterolayered YSZ/ScSZ superlattice with atomically sharp interface (100)s . Here, the model is divided into five zones where the dominantly YSZ region (labelled as (100)s YSZ), the interface region (labelled as (100)s I), and the dominantly ScSZ region (labelled as (100)s ScSZ) are shown. The Y3+ , Sc3+ , and anion vacancies are indicated by red, blue and small brown. 0.59 and 0.55 eV, respectively.

3.2.1

Lattice strain and distortions

With regards to the structural parameters of these YSZ/ScSZ [100] superlattices, we find that, at elevated temperatures, the in-plane unit cell a lattice parameter of YSZ/ScSZ is calculated to be 5.15 ˚ A at T = 1000 K, while that of bulk ScSZ is 5.13 ˚ A. Thus, under this consideration, the ScSZ-dominant region ((100)s ScSZ) of the superlattice will experience an in-plane biaxial strain. It has been reported for zirconia-based systems, lattice strain can significantly improve transport properties and enhance the stability of the strained system. 66–68

27

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

Figure 10: Region-dependent oxygen coordination number for Zr4+ , Y3+ , and Sc3+ ions (shown in black, red and blue, respectively) for the heterolayered YSZ/ScSZ superlattice: (a) The dominantly YSZ region (labelled as (100)s YSZ), (b) the interface region (labelled as (100)s I), and (c) the dominantly ScSZ region (labelled as (100)s ScSZ), respectively, for 750 ps at T = 1273 K. This lends support to the idea of using mechanical strain to tune the ionic conductivity in doped zirconia via such heterolayered interfaces for ionic transport applications. Furthermore, from the RDF analysis of YSZ/ScSZ superlattice (at T = 1273 K), the peak profile for the cation-anion and the anion-anion pairs are more similar, showing the characteristic peaks near the same positions (see Figure 9 (region by region analysis for (100)s ) and Figure S9 of the Supporting Information (all interface regions-(100)s I, (100)r I, (110)s I, (110)r I, (111)s I, and (111)r I). For instance, the first peaks of ρ(r)Zr−O in (100)s YSZ, (100)s I, and (100)s ScSZ are located at 2.106, 2.094, and 2.106 ˚ A, and those of the ρ(r)O−O are at 2.586, 2.598, and 2.586 ˚ A, respectively. Moreover, from doped ion-anion RDF analysis, the first peaks for Y-O and Sc-O are now found 2.25 and 2.22, respectively, in interface region ((100)s I) of YSZ/ScSZ whereas the corresponding values are 2.286 and 2.154 ˚ A in the heterogeneously-doped bulk YScSZ. This may indicate that the excessive lattice distortions in bulk YScSZ can be greatly reduced by forming heterolayered YSZ/ScSZ interfaces. Similarly, the RDF for two other heterointerface models in the [110] and [111] directions (denoted by (110)s , and (111)s , respectively) are also constructed for both the sharp and rough interfaces. All of the ρ(r) have peaks in similar positions (even interface regions) which are shown in Figure S9 of the Supporting

28

ACS Paragon Plus Environment

Page 28 of 41

Page 29 of 41 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

Table 5: σ values (in Scm−1 ) of YSZ/ScSZ in the [100], [110] and [111] directions (denoted by (100)s , (110)s , and (111)s , respectively) for temperatures ranging from 1000 to 1800 K. Note that for each system, five randomized structural models have been used for statistical averaging and the statistical error bar (∆σ) is also reported.

T (K) 1000 1200 1400 1600 1800

(100)s

(110)s

(111)s

σ±∆σ (Scm−1 ) 0.014 ± 0.001 0.039 ± 0.002 0.073 ± 0.001 0.118 ± 0.002 0.165 ± 0.003

σ±∆σ (Scm−1 ) 0.014 ± 0.001 0.038 ± 0.001 0.072 ± 0.002 0.115 ± 0.003 0.158 ± 0.001

σ±∆σ (Scm−1 ) 0.015 ± 0.001 0.038 ± 0.002 0.071 ± 0.002 0.115 ± 0.004 0.159 ± 0.004

Information. Following the RDF analysis, the CN of the cations for YSZ/ScSZ superlattice has been separately computed (at T = 1273 K) in the dominantly YSZ, interface, and dominantly ScSZ regions, which are illustrated in Figure 10. The CN of Zr4+ is found to be almost the same for each region. We also note that the fluctuations in the CN of Y3+ and Sc3+ in (100)s YSZ and (100)s ScSZ are mainly detected near the interface region ((100)s I). In the point of view of both lattice strain and local lattice distortion, the YSZ/ScSZ heterointerfaces may offer more stable, as well as an enhancement in the ionic conductivity, σ when compared to bulk YScSZ. The calculated values of σ for YSZ/ScSZ superlattices are found in Table 5 (in contrast to that for YScSZ in Table 4).

3.2.2

Crystal direction dependence

Now, turning our attention to the possible influence of crystal directions on ionic transport in YSZ/ScSZ, we construct two other heterointerface models in the [110] and [111] directions: (110)s and (111)s where the dz values for these two models are taken as 86.82 and 88.61 ˚ A, respectively. Using the same MD analysis approach for the [100] YSZ/ScSZ heterointerface, we have calculated the < R2 >(t), the diffusion coefficient (log D) as a function of 1000/T ,

29

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

Page 30 of 41

Table 6: σ values (in Scm−1 ) for the interface region (I) of YSZ/ScSZ in the [100], [110] and [111] directions (denoted by (100)s I, (110)s I, and (111)s I, respectively) for temperatures ranging from 1000 to 1800 K. Note that for each system, five randomized structural models have been used for statistical averaging and the statistical error bar (∆σ) is also reported.

T (K) 1000 1200 1400 1600 1800

(100)s I

(110)s I

(111)s I

σ±∆σ (S cm−1 ) 0.014 ± 0.003 0.037 ± 0.005 0.073 ± 0.004 0.118 ± 0.006 0.163 ± 0.006

σ±∆σ (S cm−1 ) 0.012 ± 0.002 0.036 ± 0.003 0.068 ± 0.004 0.114 ± 0.007 0.155 ± 0.006

σ±∆σ (S cm−1 ) 0.012 ± 0.002 0.037 ± 0.003 0.070 ± 0.007 0.113 ± 0.005 0.158 ± 0.006

and the Ea values for (110)s and (111)s . The crystal direction-resolved σ values are listed in Table 5 while their region-dependent σ values are determined for various temperatures and their interface-dominant values of σ are listed together with that of (100)s in Table 6. Our results indicate that there is a rather weak dependence on crystal orientation for ionic diffusion, and this isotropic diffusion behavior corroborates well with the isotropic trend in the orthogonally resolved < R2 >(t) of oxygen ions in the bulk doped systems in Figure 6 (d).

3.2.3

Atomic interface roughness

Lastly, we will briefly comment on the likely influence of the roughness profile of the interface in YSZ/ScSZ. We note in passing that the roughness profile of the interfaces considered to this point is taken as atomically sharp, i.e. no mixing of ions across the interface is intended (and thus, the superscript “s” for sharp interface). To mimic and model a rough interface (hereafter labeled with a superscript “r”, i.e. (100)r , (110)r , and (111)r , accordingly), we have constructed a so-called mixed interface model (see Figure S2, right) and repeated our full set of MD analysis as described above for the atomically sharp YSZ/ScSZ heterostructures detailed above.

30

ACS Paragon Plus Environment

Page 31 of 41 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

Structurally, we find very little difference in the ionic pair distributions as reflected in the RDF plots for (100)r , (110)r , and (111)r (see Figure S9 of the Supporting Information). There, the primary peaks in the RDF plots are located in very similar positions for both the sharp and rough interfaces. Interestingly, by the virtue of heterointerface engineering of these doped zirconia systems at the atomic level – regardless of the atomic roughness at the interface region – seem to lift the local lattice distortions otherwise seen in the doped bulk YScSZ. When comparing their direction-dependent and region-dependent sigma values in Tables S2 and S3 of the Supporting Information, respectively, our results seem to also support that the degree of interface roughness will not have a significant effect on the overall ionic conductivity of the heterolayered YSZ/ScSZ electrolyte.

4

Summary

In summary, using molecular dynamics simulations and diffusion dynamics analysis, we have constructed various bulk and heterointerface models of doped zirconia, namely YSZ, ScSZ, YScSZ, and YSZ/ScSZ. We have performed three main analysis to examine the structural integrity of zirconia-based materials: 1) RDF analysis for the phase, 2) the CN analysis for the vacancies, and 3) Mean square displacement (< R2 >) analysis for the cations and anions. Regarding the conductivity, we have identified that an overall dopant concentration of 8 mol % is optimal to achieve very high ionic conductivity, and in the case of bulk doped ZrO2 , ScSZ yields the best results. The calculated ionic conductivity value for YScSZ is lower than YSZ and ScSZ (rather than their intermediate value) due to the structural instability of the heteroeneously-doped system. Moreover, the oxygen ion diffusion for all systems is rather isotropic about the orthogonal < 100 > directions considered. By heterointerface engineering at the atomic scale, we have introduced a superlattice model YSZ/ScSZ where the unwanted local lattice distortions are removed and desired in-plane strain effects are induced. The overall ionic conductivity in the heterolayered YSZ/ScSZ is robust with regards to crystal

31

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

orientations and atomic roughness at the interface regions, providing a significant advantage over heterogeneously-doped bulk YScSZ in terms of stability and enhanced ionic conductivity as electrolyte materials in SOFC applications.

Supporting Information Available Supporting Information. The calculated radial distribution function and ionic conductivity of heterolayered YSZ/ScSZ superlattice systems: The illustration of the distribution of homogeneously- and inhomogeneously-doped ions; The illustration of elementary zones of the considered YSZ/ScSZ superlattice with sharp and rough interfaces; During the quenching process, radial distribution function and the energy versus time relations analyses of inhomogeneously-doped YSZ, and ScSZ and heterogeneously-doped YScSZ; Radial distribution function of pristine zirconia, inhomogeneously- and heterogeneously-doped zirconia; Energy distribution of 100,000 randomly generated configurations; Radial distribution function of O–O, Zr–O, Y–O and Sc–O pairs within sharp and rough interfaces region of superlattice (100)s I, (110)s I, (111)s I, (100)r I, (110)r I, and (111)r I; The mean square displacement of cations in heterolayered YSZ/ScSZ system; The diffusion coefficient, D value of YSZ, ScSZ, and YScSZ for temperatures from 1000 to 1800 K; The ionic conductivity values of rough interface superlattices (100)r , (110)r , and (111)r (including that for their interface regions). This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement We thank Jong-Min Yun and Ji-Hwan Lee for fruitful discussions. We also gratefully acknowledge support by Samsung Research Funding Center of Samsung Electronics under Project Number SRFC-MA1501-03. Computational resources have been provided by the KISTI Supercomputing Center (KSC-2018-C3-0006) and the Australian National Computational Infrastructure (NCI). 32

ACS Paragon Plus Environment

Page 32 of 41

Page 33 of 41 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

References (1) Minh, N. Q. Ceramic Fuel Cells. J. Am. Ceram. Soc. 1993, 76, 563–588. (2) Alcock, C.; Jacob, K.; Zador, S. Atom. Energy Rev.; 1976; pp 9, 82, 119, 124, 160, 161, 172, 173, 197. (3) Christel, P.; Meunier, A.; Heller, M.; Torre, J.; Peille, C. Mechanical Properties and Short-term in Vivo Evaluation of Yttrium-oxide-partially-stabilized Zirconia. J. Biomed. Mater. Res. A 1989, 23, 45–61. (4) Kralik, B.; Chang, E. K.; Louie, S. G. Structural Properties and Quasiparticle Band Structure of Zirconia. Phys. Rev. B 1998, 57, 7027. (5) Raghavan, S.; Wang, H.; Porter, W.; Dinwiddie, R.; Mayo, M. Thermal Properties of Zirconia Co-doped with Trivalent and Pentavalent Oxides. Acta Mater. 2001, 49, 169–179. (6) Li, P.; Chen, I.-W.; Penner-Hahn, J. E. Effect of Dopants on Zirconia StabilizationAn X-ray Absorption Study: I, Trivalent dopants. J. Am. Ceram. Soc. 1994, 77, 118–128. (7) Jung, D.-H.; Lee, J.-H.; Kilic, M. E.; Soon, A. Anisotropic Vacancy-mediated Phonon Mode Softening in Sm and Gd Doped Ceria. Phys. Chem. Chem. Phys. 2018, 20, 10048–10059. (8) Kilo, M.; Argirusis, C.; Borchardt, G.; Jackson, R. A. Oxygen Diffusion in Yttria Stabilised ZirconiaExperimental Results and Molecular Dynamics Calculations. Phys. Chem. Chem. Phys. 2003, 5, 2219–2224. (9) Costa, G.; Muccillo, R. Comparative Studies on Properties of Scandia-Sabilized Zirconia Synthesized by The Polymeric Precursor and The Polyacrylamide Techniques. J. Alloy. Comp. 2010, 503, 474–479.

33

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

Page 34 of 41

(10) Laguna-Bercero, M.; Skinner, S.; Kilner, J. Performance of Solid Oxide Electrolysis Cells Based on Scandia Stabilised Zirconia. J. Power Sources 2009, 192, 126–131. (11) Li, X.; Hafskjold, B. Molecular Dynamics Simulations of Yttrium-Stabilized Zirconia. J. Phys. Condens. Matter 1995, 7, 1255. (12) Ikeda, S.; Sakurai, O.; Uematsu, K.; Mizutani, N.; Kato, M. Electrical Conductivity of Yttria-Stabilized Zirconia Single Crystals. J. Mater. Sci. 1985, 20, 4593–4600. (13) Casselton, R. Low Field DC Conduction in Yttria-Stabilized Zirconia. Phys. Status Solidi A 1970, 2, 571–585. (14) Bogicevic, A.; Wolverton, C. Nature and Strength of Defect Interactions in Cubic Stabilized Zirconia. Phys. Rev. B 2003, 67, 024106. (15) Badwal, S.; Ciacchi, F.; Milosevic, D. Scandia–Zirconia Electrolytes for Intermediate Temperature Solid Oxide Fuel Cell Operation. Solid State Ion. 2000, 136, 91–99. (16) M¨obius,

H.-H.

Sauerstoffionenleitende

dungsm¨oglichkeiten;

Festelektrolyte

und

ihre

Anwen-

Zusammensetzung und Systematik der Festelektrolyte mit

Sauerstoffionenleitung. Zeitschrift f¨ ur Chemie 1964, 4, 81–94. (17) Yamamoto, O.; Arati, Y.; Takeda, Y.; Imanishi, N.; Mizutani, Y.; Kawai, M.; Nakamura, Y. Electrical Conductivity of Stabilized Zirconia with Ytterbia and Scandia. Solid State Ion. 1995, 79, 137–142. (18) Dixon, J.; LaGrange, L.; Merten, U.; Miller, C.; Porter, J. Electrical Resistivity of Stabilized Zirconia at Elevated Temperatures. J. Electrochem. Soc. 1963, 110, 276– 280. (19) Etsell, T.; Flengas, S. N. Electrical Properties of Solid Oxide Electrolytes. Chem. Rev. 1970, 70, 339–376.

34

ACS Paragon Plus Environment

Page 35 of 41 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

(20) Ciacchi, F.; Badwal, S.; Drennan, J. The System Y2 O3 -Sc2 O3 -ZrO2 : Phase characterisation by XRD, TEM and optical microscopy. J. Eur. Ceram. Soc. 1991, 7, 185–195. (21) Badwal, S.; Ciacchi, F.; Rajendran, S.; Drennan, J. An Investigation of Conductivity, Microstructure and Stability of Electrolyte Compositions in the System 9 mol % (Sc2 O3 –Y2 O3 )–ZrO2 (Al2 O3 ). Solid State Ion. 1998, 109, 167–186. (22) Politova, T.; Irvine, J. T. Investigation of Scandia–Yttria–Zirconia System As an Electrolyte Material for Intermediate Temperature Fuel CellsInfluence of Yttria Content in System (Y2 O3 )x (Sc2 O3 )(11−x) (ZrO2 )89 . Solid State Ion. 2004, 168, 153–165. (23) Jones, R. L.; Mess, D. Improved Tetragonal Phase Stability at 1400 C with Scandia, Yttria-Stabilized Zirconia. Surf. and Coat. Technol. 1996, 86, 94–101. (24) Cammarata, A.; Martorana, A.; Duca, D. Cation Environment of BaCeO3 -based Protonic Conductors: A Computational Study. J. Phys. Chem. A 2009, 113, 6381–6390. (25) Cammarata, A.; Emanuele, A.; Martorana, A.; Duca, D. Cation Environment of BaCeO3 - Based Protonic Conductors II: New Computational Models. J. Phys. Chem. A 2011, 115, 1676–1685. (26) Qiannan, X.; HUANG, X.; ZHANG, H.; Hong, X.; ZHANG, J.; Liangshi, W. Synthesis and Characterization of High Ionic Conductivity ScSZ Core/shell Nanocomposites. J. Rare Earths 2017, 35, 567–573. (27) Schelling, P. K.; Phillpot, S. R.; Wolf, D. Mechanism of the Cubic-to-Tetragonal Phase Transition in Zirconia and Yttria-Stabilized Zirconia by Molecular-Dynamics Simulation. J. Am. Ceram. Soc. 2001, 84, 1609–1619. (28) Zacate, M. O.; Minervini, L.; Bradfield, D. J.; Grimes, R. W.; Sickafus, K. E. Defect Cluster Formation in M2 O3 -doped Cubic ZrO2 . Solid State Ion. 2000, 128, 243–254.

35

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

(29) Grimes, R. W.; Busker, G.; McCoy, M. A.; Chroneos, A.; Kilner, J. A.; Chen, S.-P. The Effect of Ion Size on Solution Mechanism and Defect Cluster Geometry. Berichte der Bunsengesellschaft f¨ ur physikalische Chemie 1997, 101, 1204–1210. (30) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (31) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291–341. (32) Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. (Berl.) 1921, 369, 253–287. (33) Morinaga, M.; Adachi, H.; Tsukada, M. Electronic Structure and Phase Stability of ZrO2. J. Phys. Chem. Solids 1983, 44, 301–306. (34) Jansen, H. J. Electronic Structure of Cubic and Tetragonal Zirconia. Phys. Rev. B 1991, 43, 7267. (35) Nos´e, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. (36) Nos´e, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (37) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A 1985, 31, 1695. (38) Yamamura, Y.; Kawasaki, S.; Sakai, H. Molecular Dynamics Analysis of Ionic Conduction Mechanism in Yttria-Stabilized Zirconia. Solid State Ion. 1999, 126, 181–189. (39) Gonz´alez-Romero, R. L.; Mel´endez, J. J.; G´omez-Garc´ıa, D.; Cumbrera, F. L.; Dom´ınguez-Rodr´ıguez, A. A Molecular Dynamics Study of Grain Boundaries in YSZ: Structure, Energetics and Diffusion of Oxygen. Solid State Ion. 2012, 219, 1–10. 36

ACS Paragon Plus Environment

Page 36 of 41

Page 37 of 41 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

(40) Sayle, T. X.; Parker, S. C.; Sayle, D. C. Ionic Conductivity in Nano-Scale CeO 2/YSZ Heterolayers. J. Mater. Chem. 2006, 16, 1067–1081. (41) Einstein, A. On The Motion of Small Particles Suspended in Liquids at Rest Required by The Molecular-kinetic Theory of Heat. Ann. Phys. 1905, 17, 549–560. (42) Murch, G. The Haven Ratio in Fast Ionic Conductors. Solid State Ionics 1982, 7, 177–198. (43) Park, K.; Olander, D. Oxygen Diffusion in Single-Crystal Tetragonal Zirconia. J. Electrochem. Soc. 1991, 138, 1154–1159. (44) Simpson, L.; Carter, R. Oxygen Exchange and Diffusion in Calcia-Stabilized Zirconia. J. Am Cerm. Soc. 1966, 49, 139–144. ¨ (45) Arrhenius, S. Uber die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker durch S¨auren. Z. Phys. Chem. 1889, 4, 226–248. (46) Goff, J.; Hayes, W.; Hull, S.; Hutchings, M.; Clausen, K. N. Defect Structure of YttriaStabilized Zirconia and its Influence on the Ionic Conductivity at Elevated Temperatures. Phys. Rev. B 1999, 59, 14202. (47) Bechepeche, A.; Treu, O.; Longo, E.; Paiva-Santos, C.; Varela, J. A. Experimental and Theoretical Aspects of the Stabilization of Zirconia. J. Mater. Sci. 1999, 34, 2751–2756. (48) Aldebert, P.; Traverse, J. P. Structure and Ionic Mobility of Zirconia at High Temperature. J. Am. Ceram. Soc. 1985, 68, 34–40. (49) Terki, R.; Feraoun, H.; Bertrand, G.; Aourag, H. First Principles Calculations of Structural, Elastic and Electronic Properties of XO2 (X= Zr, Hf and Th) in Fluorite Phase. Comput. Mater. Sci. 2005, 33, 44–52.

37

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

(50) Kim, D.-J. Lattice Parameters, Ionic Conductivities, and Solubility Limits in Fluorite– structure MO2 oxide [M = Hf4+, Zr4+, Ce4+, Th4+, U4+] Solid Solutions. J. Am. Ceram. Soc. 1989, 72, 1415–1421. (51) Fergus, J. W. The Application of Solid Fluoride Electrolytes in Chemical Sensors. Sens. Actuator B-Chem. 1997, 42, 119–130. (52) Khan, M. S.; Islam, M. S.; Bates, D. R. Cation Doping and Oxygen Diffusion in Zirconia: A Combined Atomistic Simulation and Molecular Dynamics Sudy. J. Mater. Chem. 1998, 8, 2299–2307. (53) Arachi, Y.; Sakai, H.; Yamamoto, O.; Takeda, Y.; Imanishai, N. Electrical Conductivity of the ZrO2 –Ln2 O3 (Ln= lanthanides) system. Solid State Ion. 1999, 121, 133–139. (54) Grau-Crespo, R.; Hamad, S.; Catlow, C. R. A.; De Leeuw, N. Symmetry-adapted Configurational Modelling of Fractional Site Occupancy in Solids. J. Phys. Condens. Matter 2007, 19, 256201. (55) Dong, Y.; Qi, L.; Li, J.; Chen, I.-W. A Computational Study of Yttria-stabilized Zirconia: I. Using Crystal Chemistry to Search for the Ground State on a Glassy Energy Landscape. Acta Mater. 2017, 127, 73–84. (56) Grosso, R. L.; Muccillo, E. N.; Castro, R. H. Phase Stability in Scandia-Zirconia Nanocrystals. J. Am. Chem. Soc. 2017, 100, 2199–2208. (57) Tung, K.-L.; Chang, K.-S.; Hsiung, C.-C.; Chiang, Y.-C.; Li, Y.-L. Molecular Dynamics Simulation of the Complex Dopant Effect on the Super-ionic Conduction and Microstructure of Zirconia-based Solid Electrolytes. Sep. Purif. Technol. 2010, 73, 13– 19. (58) Omar, S.; Bin Najib, W.; Chen, W.; Bonanos, N. Electrical Conductivity of 10 mol% Sc2 O3 –1 mol% M2 O3 –ZrO2 Ceramics. J. Am. Ceram. Soc. 2012, 95, 1965–1972. 38

ACS Paragon Plus Environment

Page 38 of 41

Page 39 of 41 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

(59) Dong, Y.; Qi, L.; Li, J.; Chen, I.-W. A Computational Study of Yttria-stabilized Zirconia: II. Cation Diffusion. Acta Mater. 2017, 126, 438–450. (60) Kilo, M.; Borchardt, G.; Lesage, B.; Kaıtasov, O.; Weber, S.; Scherrer, S. Cation Transport in Yttria Stabilized Cubic Zirconia: 96Zr tracer diffusion in (ZrxY1–x) O2– x/2 single crystals with 0.15? x? 0.48. J. Eur. Ceram. Soc. 2000, 20, 2069–2077. (61) 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. (62) Kamiya, M.; Shimada, E.; Ikuma, Y.; Komatsu, M.; Haneda, H. Intrinsic and Extrinsic Oxygen Diffusion and Surface Exchange Reaction in Cerium Oxide. J. Electrochem. Soc. 2000, 147, 1222–1227. (63) Brinkman, H.; Briels, W. J.; Verweij, H. Molecular Dynamics Simulations of YttriaStabilized Zirconia. Chem. Phys. Lett. 1995, 247, 386–390. (64) Araki, W.; Arai, Y. Oxygen Diffusion in Yttria-Stabilized Zirconia Subjected to Uniaxial Stress. Solid State Ion. 2010, 181, 441–446. (65) Devanathan, R.; Weber, W. J.; Singhal, S. C.; Gale, J. D. Computer Simulation of Defects and Oxygen Transport in Yttria-Stabilized Zirconia. Solid State Ion. 2006, 177, 1251–1258. (66) Suzuki, K.; Kubo, M.; Oumi, Y.; Miura, R.; Takaba, H.; Fahmi, A.; Chatterjee, A.; Teraishi, K.; Miyamoto, A. Molecular Dynamics Simulation of Enhanced Oxygen Ion Diffusion in Strained Yttria-stabilized Zirconia. Appl. Phys. Lett. 1998, 73, 1502–1504. (67) Jiang, J.; Hu, X.; Shen, W.; Ni, C.; Hertz, J. L. Improved Ionic Conductivity in Strained Yttria-stabilized Zirconia Thin Films. Appl. Phys. Lett. 2013, 102, 143901.

39

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

(68) Kushima, A.; Yildiz, B. Oxygen Ion Diffusivity in Strained Yttria Stabilized Zirconia: Where is The Fastest Strain? J. Mater. Chem. 2010, 20, 4809–4819.

40

ACS Paragon Plus Environment

Page 40 of 41

Page 41 of 41 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

Graphical TOC Entry

41

ACS Paragon Plus Environment