**3.1.2 Structural features**

#### **i. Crystallographic direction views**

Besides the different thermal responses of the lattices, a more illuminating way of exhibiting the qualitative differences between the YSZ crystals and amorphous YSZ solids can be shown by the structural features determined by the ion distribution within the lattices of these solids.

Bulk YSZ is a solid solution on a cubic fluorite lattice with yttrium and zirconium distributed on a face-centered cubic cation lattice and oxygen and vacancies distributed on a simple-cubic anion lattice. Its fully dense solid features are revealed at the top of Fig. 5. To differentiate a dense crystalline lattice with a disordered solid lattice, the Miller index that defines a crystallographic direction for the planes and directions in the crystal lattice can be very useful. By randomly distributing the Y3+ within a cubic Zr(1−*x)*Y*x*O(2−*x/*2) phase, the distribution of the Y3+ ions in both the YSZ crystals (e.g. 8-YSZc) and the amorphous solids (e.g. 8-YSZa) are found to be nearly isotropic as shown in Fig. 5. For 8-YSZc in its cubic fluorite structure, the ionic distributions in all three cubic (001), (010), and (001) directions are found to be equivalent in all cases. These directions for the amorphous cases correspond to those directions in the original crystal that was made amorphous by heating in vacuum and then repeatedly heating and cooling under hydrostatic pressure. In contrast to the YSZ amorphous systems, the YSZ crystals (e.g. 8-YSZc in Fig. 5) exhibit comparatively clear Zr4+ (cation), O2− (anion) crystalline planes in both directions shown in Fig. 5 for the crystal lattice. The dopant (Y3+) ions have a similar distribution to that of the Zr4+ ions that they replace in the crystal of YSZ (Fig. 5). For both systems (i.e. the crystal and amorphous structure), the ratio of Zr4+/O2−~ 1/2 can be found in all crystallographic planes, analogous to the ZrO2 stoichiometry which acts as a host lattice. In addition, these similar features are also found in the YSZ systems with other Y2O3 concentrations (i.e. ~ 3.0 - 12.0 mol%) and for elevated temperatures, as long as there is no grain boundary in the crystal.

#### **ii. Radial distribution function**

We analyze the bonding features of the simulated n-YSZc and n-YSZa systems. Even for ions in a disordered-network of an amorphous solid with no long range order, the short range order, however, can be obtained from the distribution of coordination number, bond length, and bond angles. Interpreting the static equilibrium structural features of these YSZ systems in another way, we calculated the radial distribution function (RDF), *g(r)* of the

thermostats (*NPT*) for further solidification and recrystallization, the MD simulation will generally yield a cubic structure (Lau et al., 2011). In the case of 8.0% Y2O3 doping in amorphous YSZ (i.e. 8-YSZa), the temperature dependence of the 8-YSZa lattice can be fitted to the expression *a(T)*=*a(*300*)*[1+*αa(T*−300*)*], and *c(T)*=*c(*300*)*[1+*αc(T*−300*)*], with *αa* and *α<sup>c</sup>* corresponding to linear thermal expansion coefficients along the *a* and *c* lattice directions. The *αa* and *αc* are found to be nearly identical for each *n-YSZa* system. Overall, the *α* for the 3- YSZa, 8-YSZa and 12-YSZa solids are found to be ~1*.*5–2*.*2×10-6 K-1, which are significantly smaller than for the corresponding YSZ crystals (Lau et al., 2011). The amorphous YSZ solids generally have a much smaller volume expansion over the temperature range shown in Fig. 4, which might be attributed to extra flexibility in spatial rearrangements for the ions

Besides the different thermal responses of the lattices, a more illuminating way of exhibiting the qualitative differences between the YSZ crystals and amorphous YSZ solids can be shown by the structural features determined by the ion distribution within the lattices of

Bulk YSZ is a solid solution on a cubic fluorite lattice with yttrium and zirconium distributed on a face-centered cubic cation lattice and oxygen and vacancies distributed on a simple-cubic anion lattice. Its fully dense solid features are revealed at the top of Fig. 5. To differentiate a dense crystalline lattice with a disordered solid lattice, the Miller index that defines a crystallographic direction for the planes and directions in the crystal lattice can be very useful. By randomly distributing the Y3+ within a cubic Zr(1−*x)*Y*x*O(2−*x/*2) phase, the distribution of the Y3+ ions in both the YSZ crystals (e.g. 8-YSZc) and the amorphous solids (e.g. 8-YSZa) are found to be nearly isotropic as shown in Fig. 5. For 8-YSZc in its cubic fluorite structure, the ionic distributions in all three cubic (001), (010), and (001) directions are found to be equivalent in all cases. These directions for the amorphous cases correspond to those directions in the original crystal that was made amorphous by heating in vacuum and then repeatedly heating and cooling under hydrostatic pressure. In contrast to the YSZ amorphous systems, the YSZ crystals (e.g. 8-YSZc in Fig. 5) exhibit comparatively clear Zr4+ (cation), O2− (anion) crystalline planes in both directions shown in Fig. 5 for the crystal lattice. The dopant (Y3+) ions have a similar distribution to that of the Zr4+ ions that they replace in the crystal of YSZ (Fig. 5). For both systems (i.e. the crystal and amorphous structure), the ratio of Zr4+/O2−~ 1/2 can be found in all crystallographic planes, analogous to the ZrO2 stoichiometry which acts as a host lattice. In addition, these similar features are also found in the YSZ systems with other Y2O3 concentrations (i.e. ~ 3.0 - 12.0 mol%) and for

elevated temperatures, as long as there is no grain boundary in the crystal.

We analyze the bonding features of the simulated n-YSZc and n-YSZa systems. Even for ions in a disordered-network of an amorphous solid with no long range order, the short range order, however, can be obtained from the distribution of coordination number, bond length, and bond angles. Interpreting the static equilibrium structural features of these YSZ systems in another way, we calculated the radial distribution function (RDF), *g(r)* of the

in disordered solids.

these solids.

**3.1.2 Structural features** 

**i. Crystallographic direction views** 

**ii. Radial distribution function** 

constituent ions, which describes the average spatial organization of cations and anions in the lattice. A Fourier transform of the radial distribution function results in the structure factor *S(q)*, which is experimentally measurable from X-ray or neutron scattering. Therefore, it is an important result from the predicted static structure of a material. Based on the prominent peaks within the short-range order (i.e. < 6 Å), the average local connectivity of the Zr4+, Y3+ and O2- ions in the system can be determined. From the RDF, the successive peaks correspond to the nearest-, the second- and the next-neighboring atomic distributions inside a YSZ system. The distinct RDFs of the YSZ crystal and YSZ amorphous solid (Fig. 6) are due to their unique local ion distributions and system densities.

Fig. 5. The 8-YSZc crystal (top) and the 8-YSZa amorphous structure (bottom). Arrows show the projection vectors in the (001) and (111) directions of these simulation cells. The right two columns of figures give the distribution of the three ions binned by atomic planes, every two angstroms, along these two crystallographic directions (Lau et al., 2011).

For 8-YSZc at ~ 300 K, the average nearest-neighbor Zr-O, Y-O, O-O, and Zr-Y bond distances are ~ 2*.*08, 2.33, 2.58, and 3*.*58 Å, respectively (Fig. 6a). These features generally are very similar to the RDF features of 14-YSZc computed using the more sophisticated ReaxFF Reactive Force Field method (van Duin et al., 2008) (Fig. 6b). For the prediction of the Zr-O distance in YSZ crystals, both interatomic potentials yield excellent agreement with EXAFS (i.e. 2.13 Å) and neutron diffraction (i.e. 2.08 Å) for the 15-YSZ crystal (van Duin et al., 2008).

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 355

EXAFS data, 2.32 Å (Ishizawa et al., 1999). Integrating the Y-O RDF to its first minimum (~3*.*22 Å), gives coordination number *N*Y = 6*.*6, slightly larger than the *N*Y = 6 of the cubic

For the O-O pair distribution, its more diffusive behavior compared the heavier cations yields a broader pair distribution over the range of O-O distances of *g(r)* that is plotted (Fig. 6a). With its O-O average distance found to be ~ 2.58 Å at T = 300 K, the average oxygen atom has *N*O = 6*.*24 (integrating up to the first O-O minimum at 2*.*98 Å), which is almost the same as in both ZrO2 and Y2O3. As the temperature rises, all the ions tend to move on average farther away from their equilibrium positions which give the thermal volume expansion of the system. All the peaks of *g(r)* in the RDF are generally lowered relative to room temperature for all the ionic pairs in a system, as a result of thermal broadening and

Similar observations are found in the ReaxFF-computed RDF of 14-YSZc (Fig. 6b). For the O-O pair distributions in 14-YSZc at T =1000 and 2000 K, the average O-O bond distance is ~ 2.67 Å at T = 1000 K, increasing to ~ 2.77 Å at T = 2000 K. Thus this suggested that for both the Buckingham potential and the ReaxFF Reactive Force Field, the description of static structural features is found to be equivalent for heavy ions (Zr4+ and Y3+) and mobile ions (O2−). As shown from Fig. 6a and Fig. 6b, for different Y2O3 compositions (from 3.0 to 14.0 mol% Y2O3), their ion pair distributions generally remain similar, as the cubic crystal lattice features are preserved. However, higher concentration of the Y3+ dopant means more Zr4+ ions were substituted; hence, the heights of the Zr-Zr *g(r )* peaks are reduced with increasing Y3+ concentration. Consistent with *N*O ~ 6*.*24 for both ZrO2 and Y2O3, the first peak height for the O-O *g(r)* profile is independent of dopant concentration, but higher Y3+ concentration requires more oxygen vacancies in the system, and therefore lowers the height of *g(r)* for

For the YSZ amorphous solid, disordered structural features are reflected in its RDF. In general, the peaks of the RDF spectrum for all ion pairs at small interatomic separations are broadened relative to the crystal. This broadening increases with separation, and structure almost disappears at the largest separations plotted. This is a clear indication of disorder, which increases with distance. Instead of having several pronounced peaks which preserve long-range order in a perfect crystal, there is no significant structure beyond ~ 8*.*0 Å for the Zr-O, Y-O, Zr-Y, and O-O pair separations (Fig. 6c). The broad profile of these ion pairs indicates that these amorphous solids have no long-range order for ions separated beyond ~ 8*.*0 Å. Based on the discrete peaks within the short-range order (i.e. *<* 8*.*0 Å), the local connectivity of these ions however can be examined. For these ion pairs, the average nearest-neighbor Zr-O, Y-O, O-O, and Zr-Y bond distances are found to be 1.98, 2.18, 2.78, and 3*.*53 Å, respectively, for the 8-YSZa solid. Despite having the same Y3+ composition as 8- YSZc in Fig. 6a, the substantial increase in volume (i.e. 8.6% volume swelling) together with loss of long-range order cause the average coordination number of Zr, *N*Zr to be ~ 5.9, substantially less than *N*Zr = 7*.*6 for the crystal represented by 8-YSZc. Similarly for *N*Y (i.e. by integrating the Y-O RDF to the minimum at ~3*.*03 Å), the *N*Y is 5.6, less than both *N*Y ~ 6*.*6 in 8-YSZc and *N*Y = 6 in the cubic Y2O3 crystal, which also might attributed to the slight expansion in the volume of amorphous YSZ. However, with increased temperature the

peaks of *g(r)* for the amorphous state broaden as they do in the perfect crystal.

Y2O3 crystal.

increased mobility.

more distant O-O pairs.

Fig. 6. (a) The RDF of 8-YSZc at 300 K based on the Buckingham (BMB) potential as defined in Eq. 1 using parameters from Table 1. (b) The RDF of 14-YSZc at 1000 K based on the ReaxFF Reactive Force Field adopted from van Duin et al. (van Duin et al., 2008). (c) The RDF of 8-YSZa at 300 K based on the BMB potential from Table 1.

Integrating the Zr-O RDF up to the first minimum ~ 3*.*27 Å past the peak at 2*.*08 Å, gives an average coordination number *N*Zr = 7*.*6, in between that of *N*zr = 7 for the monoclinic ZrO2, and *N*Zr = 8 for the tetragonal ZrO2 phase. From ReaxFF computed 14-YSZc, the RDF for the Y-O indicates an average Y-O distance to be ~ 2.32 Å (Fig. 6b), which is close to Y-O distance of 8-YSZc (i.e. ~ 2.33 Å) computed based on the BMB potential and experimental

Fig. 6. (a) The RDF of 8-YSZc at 300 K based on the Buckingham (BMB) potential as defined in Eq. 1 using parameters from Table 1. (b) The RDF of 14-YSZc at 1000 K based on the ReaxFF Reactive Force Field adopted from van Duin et al. (van Duin et al., 2008). (c) The

Integrating the Zr-O RDF up to the first minimum ~ 3*.*27 Å past the peak at 2*.*08 Å, gives an average coordination number *N*Zr = 7*.*6, in between that of *N*zr = 7 for the monoclinic ZrO2, and *N*Zr = 8 for the tetragonal ZrO2 phase. From ReaxFF computed 14-YSZc, the RDF for the Y-O indicates an average Y-O distance to be ~ 2.32 Å (Fig. 6b), which is close to Y-O distance of 8-YSZc (i.e. ~ 2.33 Å) computed based on the BMB potential and experimental

RDF of 8-YSZa at 300 K based on the BMB potential from Table 1.

For the O-O pair distribution, its more diffusive behavior compared the heavier cations yields a broader pair distribution over the range of O-O distances of *g(r)* that is plotted (Fig. 6a). With its O-O average distance found to be ~ 2.58 Å at T = 300 K, the average oxygen atom has *N*O = 6*.*24 (integrating up to the first O-O minimum at 2*.*98 Å), which is almost the same as in both ZrO2 and Y2O3. As the temperature rises, all the ions tend to move on average farther away from their equilibrium positions which give the thermal volume expansion of the system. All the peaks of *g(r)* in the RDF are generally lowered relative to room temperature for all the ionic pairs in a system, as a result of thermal broadening and increased mobility.

Similar observations are found in the ReaxFF-computed RDF of 14-YSZc (Fig. 6b). For the O-O pair distributions in 14-YSZc at T =1000 and 2000 K, the average O-O bond distance is ~ 2.67 Å at T = 1000 K, increasing to ~ 2.77 Å at T = 2000 K. Thus this suggested that for both the Buckingham potential and the ReaxFF Reactive Force Field, the description of static structural features is found to be equivalent for heavy ions (Zr4+ and Y3+) and mobile ions (O2−). As shown from Fig. 6a and Fig. 6b, for different Y2O3 compositions (from 3.0 to 14.0 mol% Y2O3), their ion pair distributions generally remain similar, as the cubic crystal lattice features are preserved. However, higher concentration of the Y3+ dopant means more Zr4+ ions were substituted; hence, the heights of the Zr-Zr *g(r )* peaks are reduced with increasing Y3+ concentration. Consistent with *N*O ~ 6*.*24 for both ZrO2 and Y2O3, the first peak height for the O-O *g(r)* profile is independent of dopant concentration, but higher Y3+ concentration requires more oxygen vacancies in the system, and therefore lowers the height of *g(r)* for more distant O-O pairs.

For the YSZ amorphous solid, disordered structural features are reflected in its RDF. In general, the peaks of the RDF spectrum for all ion pairs at small interatomic separations are broadened relative to the crystal. This broadening increases with separation, and structure almost disappears at the largest separations plotted. This is a clear indication of disorder, which increases with distance. Instead of having several pronounced peaks which preserve long-range order in a perfect crystal, there is no significant structure beyond ~ 8*.*0 Å for the Zr-O, Y-O, Zr-Y, and O-O pair separations (Fig. 6c). The broad profile of these ion pairs indicates that these amorphous solids have no long-range order for ions separated beyond ~ 8*.*0 Å. Based on the discrete peaks within the short-range order (i.e. *<* 8*.*0 Å), the local connectivity of these ions however can be examined. For these ion pairs, the average nearest-neighbor Zr-O, Y-O, O-O, and Zr-Y bond distances are found to be 1.98, 2.18, 2.78, and 3*.*53 Å, respectively, for the 8-YSZa solid. Despite having the same Y3+ composition as 8- YSZc in Fig. 6a, the substantial increase in volume (i.e. 8.6% volume swelling) together with loss of long-range order cause the average coordination number of Zr, *N*Zr to be ~ 5.9, substantially less than *N*Zr = 7*.*6 for the crystal represented by 8-YSZc. Similarly for *N*Y (i.e. by integrating the Y-O RDF to the minimum at ~3*.*03 Å), the *N*Y is 5.6, less than both *N*Y ~ 6*.*6 in 8-YSZc and *N*Y = 6 in the cubic Y2O3 crystal, which also might attributed to the slight expansion in the volume of amorphous YSZ. However, with increased temperature the peaks of *g(r)* for the amorphous state broaden as they do in the perfect crystal.

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 357

temperature range is found consistent with the reported experimental data (Kilo et al., 2003), as highlighted in Fig. 7a. From the data points of the present MD simulations for a period up to 2.5 ns for 300–1400 K, the pre-exponential diffusion constant, *D*0, and activation energy, *ΔE*act, for oxygen-ion diffusion in 8-YSZc are found to be ~5*.*83±0*.*47×10-5 cm2s-1 and 0*.*59±0*.*05 eV, respectively, in close agreement with previously reported 6*.*08×10-5 cm2s-1, and 0.60 eV theoretical results (Devanathan et al., 2006), but slightly smaller than ~ 1*.*0 eV for 10.0 mol*%* Y2O3 YSZ reported from experiment (Kilo et al., 2003). Similarly for the other YSZ crystals in the range of 3.0-12.0 mol*%* Y2O3 of this study, the predicted activation energy, *ΔE*act, for oxygen ions is ~0*.*57–0*.*68 eV, within the range ~0*.*2–1*.*0 eV of those reported in previous theoretical studies (Devanathan et al., 2006; Khan et al., 1998; Kilo et al., 2003; Li et al., 1995; van Duin et al., 2008) at various Y2O3 concentrations, including the prediction obtained from the more accurate ReaxFF Reactive Force Field (van Duin et al., 2008). From Fig. 7b our ionic

Fig. 7. The predicted oxygen ion transport properties of YSZ crystals and amorphous solids within the range ~ 300-1400 K based on the simple interatomic potential (BMB potential) from Table 1. (a) The self-diffusion coefficient *D*i (in cm2s-1) of O2− for 8-YSZc as a function of reciprocal temperature, *1/T* (in K-1) compared to the reported experimental observation for 10.0 mol% Y2O3 YSZ. (b) The temperature-dependent ionic conductivity, *σ* (in Scm-1) of YSZa and YSZ-c compared to experimental values of dc electrical conductivity measured on 8.0 mol% YSZ amorphous (i.e. 8-YSZ Exp. I) and crystalline YSZ thin films (i.e. 8-YSZ Exp. II)

Based on the MD simulations for periods up to the same 2.5 ns for 300–1000 K, the oxygenion conductivity in amorphous YSZ solids is generally found to be comparable to that of the corresponding YSZ crystals. This trend is also found to be consistent with the reported temperature-dependent dc-electrical conductivity measured on the YSZ crystalline and

conductivity is maximal at 8 mol% Y2O3.

(Lau et al., 2011).

#### **3.2 Dynamic properties**

From MD simulation, the time-dependent trajectories of the ions of a system within a time interval often provide a lot of useful information, e.g. ion diffusion patterns, ion hopping mechanisms, vibrational frequencies, transport properties (e.g. ionic and thermal), etc., that govern the diverse phenomena and the evolution of the dynamic properties of a system. For the YSZ, the most useful property is the ionic transport property that determines the electrical conductivity of an operating SOFC. The method to compute this property will be discussed in this section.

#### **3.2.1 Ion diffusion and ionic transport in YSZ solids**

#### **i. Temperature-dependent self diffusion and ionic conductivity of O2- ions**

As we know, diffusion means "spreading". This quantity can be either observable (physical) or abstract and probabilistic (stochastic) in nature. In the context of YSZ, diffusion refers to the basic dynamic properties of the collective motion of its different constituent ions in response to temperature, which are governed by the self-diffusion of the ions. In general, the self-diffusion of each type of ion is governed by thermally activated random walks in the local structure, which for oxygen ions are primarily determined by the surrounding occupied and vacant sites. At low frequencies, the existence of direct-current (dc) conductivity implies that the mean square displacement (MSD) of ions is linear, i.e. <r2(t)>dc ~ *Dt* (Sidebottom, 2009). Linear time-dependent behavior is just a reflection of random diffusion of the ions as they migrate from site to site stochastically, which is a hallmark of uncorrelated motion. To relate experimental observation and computed atomic trajectories, MSD is a key quantity of interest in quantifying the motion of ions, and it is defined as follows:

$$\text{\textbullet r}\_{i}(\text{t}) \succ = \Sigma \text{\textbullet } [\text{r}\_{i}(\text{t}) \text{-r}\_{i}(\text{0})]^{2} / \text{N} \tag{2}$$

Within the linear response regime (Dyre et al., 2009) for an isotropic medium such as a cubic lattice of YSZ, one can extract the diffusion constant *Di* for species *i* from the MSD using the Einstein relation, <*r*2*i (t)* > = 6*Dit+C*. Thus, the diffusion coefficients can be evaluated from the slope of the MSD curves versus elapsed time in a steady-state MD simulation. The direct current (dc) ionic conductivity, *σi*, of ion species *i* can be estimated using the Nernst–Einstein relation: *σi* = *Niq2 iDi/HrkBT,* where *Ni* is the charge density of charge carriers, *qi*, per unit volume, *H*R is the Haven ratio (March, 1982) (for simplicity, we chose *H*R = 1*.*0 in this work), *k*B is the Boltzmann constant, and *T* is the system temperature.

At low temperature (e.g. *T* << 300 K), all constituent ions within the YSZ lattices (i.e. both crystalline and amorphous solids) are mostly found to be 'frozen' at a metastable lattice sites with negligible hopping. However, the thermally activated ionic hopping within the lattices is essentially governed by the kinetic energy, which in turn is determined by system temperature and ionic mass. At a given temperature the heavier species (e.g. cations like Zr4+, Y3+) are less mobile than the anions (e.g. O2-, or equivalently its vacancy). As an oxygen-ion conductor at elevated temperatures, the basic feature of oxygen-ion transport within YSZ is generally well described by the Arrhenius relationship for activated processes, relationship *Di* = *D*0exp*(-ΔE*act/*kT)*. Based on the simple BMB potential (Table 1), the leastsquares linear fit to the computed self-diffusion of oxygen in the 8-YSZc system over the

From MD simulation, the time-dependent trajectories of the ions of a system within a time interval often provide a lot of useful information, e.g. ion diffusion patterns, ion hopping mechanisms, vibrational frequencies, transport properties (e.g. ionic and thermal), etc., that govern the diverse phenomena and the evolution of the dynamic properties of a system. For the YSZ, the most useful property is the ionic transport property that determines the electrical conductivity of an operating SOFC. The method to compute this property will be

As we know, diffusion means "spreading". This quantity can be either observable (physical) or abstract and probabilistic (stochastic) in nature. In the context of YSZ, diffusion refers to the basic dynamic properties of the collective motion of its different constituent ions in response to temperature, which are governed by the self-diffusion of the ions. In general, the self-diffusion of each type of ion is governed by thermally activated random walks in the local structure, which for oxygen ions are primarily determined by the surrounding occupied and vacant sites. At low frequencies, the existence of direct-current (dc) conductivity implies that the mean square displacement (MSD) of ions is linear, i.e. <r2(t)>dc ~ *Dt* (Sidebottom, 2009). Linear time-dependent behavior is just a reflection of random diffusion of the ions as they migrate from site to site stochastically, which is a hallmark of uncorrelated motion. To relate experimental observation and computed atomic trajectories, MSD is a key quantity of interest in quantifying the motion of ions, and it is

Within the linear response regime (Dyre et al., 2009) for an isotropic medium such as a cubic lattice of YSZ, one can extract the diffusion constant *Di* for species *i* from the MSD using the Einstein relation, <*r*2*i (t)* > = 6*Dit+C*. Thus, the diffusion coefficients can be evaluated from the slope of the MSD curves versus elapsed time in a steady-state MD simulation. The direct current (dc) ionic conductivity, *σi*, of ion species *i* can be estimated using the Nernst–Einstein

volume, *H*R is the Haven ratio (March, 1982) (for simplicity, we chose *H*R = 1*.*0 in this work),

At low temperature (e.g. *T* << 300 K), all constituent ions within the YSZ lattices (i.e. both crystalline and amorphous solids) are mostly found to be 'frozen' at a metastable lattice sites with negligible hopping. However, the thermally activated ionic hopping within the lattices is essentially governed by the kinetic energy, which in turn is determined by system temperature and ionic mass. At a given temperature the heavier species (e.g. cations like Zr4+, Y3+) are less mobile than the anions (e.g. O2-, or equivalently its vacancy). As an oxygen-ion conductor at elevated temperatures, the basic feature of oxygen-ion transport within YSZ is generally well described by the Arrhenius relationship for activated processes, relationship *Di* = *D*0exp*(-ΔE*act/*kT)*. Based on the simple BMB potential (Table 1), the leastsquares linear fit to the computed self-diffusion of oxygen in the 8-YSZc system over the

*k*B is the Boltzmann constant, and *T* is the system temperature.

*iDi/HrkBT,* where *Ni* is the charge density of charge carriers, *qi*, per unit

i(t)> = ΣNi[ri(t)-ri(0)]2/N (2)

**i. Temperature-dependent self diffusion and ionic conductivity of O2- ions** 

**3.2 Dynamic properties** 

discussed in this section.

defined as follows:

relation: *σi* = *Niq2*

<r2

**3.2.1 Ion diffusion and ionic transport in YSZ solids** 

temperature range is found consistent with the reported experimental data (Kilo et al., 2003), as highlighted in Fig. 7a. From the data points of the present MD simulations for a period up to 2.5 ns for 300–1400 K, the pre-exponential diffusion constant, *D*0, and activation energy, *ΔE*act, for oxygen-ion diffusion in 8-YSZc are found to be ~5*.*83±0*.*47×10-5 cm2s-1 and 0*.*59±0*.*05 eV, respectively, in close agreement with previously reported 6*.*08×10-5 cm2s-1, and 0.60 eV theoretical results (Devanathan et al., 2006), but slightly smaller than ~ 1*.*0 eV for 10.0 mol*%* Y2O3 YSZ reported from experiment (Kilo et al., 2003). Similarly for the other YSZ crystals in the range of 3.0-12.0 mol*%* Y2O3 of this study, the predicted activation energy, *ΔE*act, for oxygen ions is ~0*.*57–0*.*68 eV, within the range ~0*.*2–1*.*0 eV of those reported in previous theoretical studies (Devanathan et al., 2006; Khan et al., 1998; Kilo et al., 2003; Li et al., 1995; van Duin et al., 2008) at various Y2O3 concentrations, including the prediction obtained from the more accurate ReaxFF Reactive Force Field (van Duin et al., 2008). From Fig. 7b our ionic conductivity is maximal at 8 mol% Y2O3.

Fig. 7. The predicted oxygen ion transport properties of YSZ crystals and amorphous solids within the range ~ 300-1400 K based on the simple interatomic potential (BMB potential) from Table 1. (a) The self-diffusion coefficient *D*i (in cm2s-1) of O2− for 8-YSZc as a function of reciprocal temperature, *1/T* (in K-1) compared to the reported experimental observation for 10.0 mol% Y2O3 YSZ. (b) The temperature-dependent ionic conductivity, *σ* (in Scm-1) of YSZa and YSZ-c compared to experimental values of dc electrical conductivity measured on 8.0 mol% YSZ amorphous (i.e. 8-YSZ Exp. I) and crystalline YSZ thin films (i.e. 8-YSZ Exp. II) (Lau et al., 2011).

Based on the MD simulations for periods up to the same 2.5 ns for 300–1000 K, the oxygenion conductivity in amorphous YSZ solids is generally found to be comparable to that of the corresponding YSZ crystals. This trend is also found to be consistent with the reported temperature-dependent dc-electrical conductivity measured on the YSZ crystalline and

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 359

oxygen vacancy preferentially occupies interstitial sites adjacent to the dopant ions, implying significant dopant–vacancy attraction. However the extended X-ray absorption fine structure spectroscopy (EXAFS) and electron microscopy (Allpress et al., 1975; Catlow et al., 1986; Veal et al., 1988) studies indicates that oxygen vacancies (i.e. positively charged) have a greater preference for sitting in the first neighbor shell of Zr4+ ions than Y3+ ions. Hence, the dopant ion has a higher oxygen coordination number than the host Zr4+ ions. This argument suggests that the oxygen vacancy–Y-ion interaction is also dominated by size effects rather than a simple Coulomb interactions (note, the vacancy is positively charged and the Y3+-ion is negatively charged, relative to the host lattice, Zr4+-ion). Therefore to capture the Y2O3-dopant dependent ionic conductivity behavior, a small simulation cell might not be sufficient. In the meantime for YSZ, when an oxygen vacancy is formed, the lattice contracts in the vicinity of the vacancy. Because Zr4+ ions are smaller than Y3+ ions, Zr4+ ions can better accommodate this contraction. In this case, a simple lattice statics calculations based on a Born–Mayer interatomic potential

As shown in Fig. 8c, the simple interatomic potential from Table 1 can provide a consistent qualitative feature, a conductivity maximum, which is observed in experiments (Arachi et al., 1999; Badwal et al., 1992; Hull, 2009; Ioffe, 1978) and Kinetic Monte Carlo simulation in Fig. 8a (Krishnamurthy et al., 2004). For all three temperatures in the range 800–1200 K, the variation of conductivity, σ, with dopant concentration shows a maximum close to 8.0 mol% Y2O3 (or x ~0.16 in Zr(1−x)YxO(2-x/2)), consistent with finding the maximum close to the lower limit of the amount of doping required to stabilize the cubic YSZ phase (i.e. *c\**-YSZ) (Arachi

For YSZ amorphous solids, a similar trend of the optimal Y2O3 concentration governing ionic conductivity is also supported by a recent MD simulation (Lau et al., 2011). Of the 3.0, 8.0 and 12.0 mol% Y2O3 systems (i.e. n-YSZa) studied; it was found that σmax occurs at ~ 8.0 mol% Y2O3 similar to the maximum found for the YSZ crystals. To account for this intriguing result, one simple, yet direct explanation can be the local Y3+ aggregation within the system (Lau et al., 2011). The diffusion pathways of an ion are essentially determined by the local interactions with the surrounding ions in YSZ at a given temperature. Due to the cooperative effects of larger attractive terms (Aij and ρij in Tables 1) in the BMB potential, larger bond length (i.e. RY−O from RDF at Fig. 6) and larger Y3+ ion relative to Zr4+ ion in the lattice, steeper potential wells among Y3+–O2− relative to Zr4+–O2− are expected. Thus, this will hinder the mobility of O2− (or its vacancy), if the local Y3+ aggregation is large enough. At a given temperature, the local Y3+ aggregation can be defined qualitatively by comparing to the degree of Y3+-clustering, *η* (Table 2), based on the thermal equilibrium structures of each system. For each of the *N*Y atoms, we count the number of other Y atoms less than 3*.*7Å away. If *Nc* is the number of these Y atoms that have four or more neighboring Y atoms, then *η* = *Nc/N*×100*.*0*%*. As shown in Table 2, *η* increased as the mol% of Y2O3 increased regardless of crystal or amorphous geometry. As the system temperature goes up, *η* generally decreases due to thermal expansion of the lattice. Despite this change with temperature, the relative comparison of each system remains the same. Therefore, this suggests that YSZ ionic conductivity is a thermally activated process as shown in Fig. 7 and 8, but the subtle interplay between the vacancy concentration and the intrinsic ionic diffusion pathways that

are defined by local structures and local orderings cannot be ignored.

also support this view (Dwivedi et al., 1990).

et al., 1999; Hull, 2009; Nakamura et al., 1986).

amorphous thin films (Heiroth et al., 2008). Fitting to the Arrhenius relationship over the temperature range, a pre-exponential *D*0 for the 8-YSZa of ~3*.*56±0*.*29×10-5 cm2s-1 is obtained, whereas ~5*.*83±0*.*47×10-5 cm2s-1 was found for crystalline 8-YSZc (Lau et al., 2011). In general, the pre-exponential *D*0 of amorphous YSZ solids are found to be slightly less compared to their respective perfect crystals. This might be attributed to their distinct local structural properties, which subsequently affect the underlying random walks of the mobile O2− ions. Specifically, the Zr4+ and Y3+ ions are found to be lower in coordination number to the O2<sup>−</sup> ions in the irregular structural geometries (Sect. 3.1.2.II). Thus the effective probability of an oxygen ion finding a vacant site, relative to these reduced cation coordination numbers, in the amorphous lattices is increased relative to a highly coordinated, compact perfect cubic YSZ crystal. The effective activation energy, *ΔE*act, of these YSZ amorphous (3.0, 8.0 and 12.0 mol% Y2O3) solids is ~ 0*.*45–0*.*61 eV, comparable to YSZ crystals with the same Y2O3 concentration. For the 8-YSZa solid, the *ΔE*act is slightly smaller (i.e. 0*.*45±0*.*03 eV) compared to its crystal counterpart. The differences might be due to a slight variation in ionic transport within the lattices. For YSZ crystals, the ion diffusion is mainly dominated by thermally activated hopping processes of the mobile anions and vacancies, whereas for YSZ amorphous solids, there is an additional, non-negligible mutual diffusion that involves cations and anions moving together in the expanded lattices (Sect. 3.2.3).

#### **ii. Y2O3-dopant concentration dependent self diffusion and ionic conductivity of O2- ions**

As anticipated in Fig. 7, the high ionic conductivity in YSZ solids (both crystalline and amorphous) generally occur at rather elevated temperatures. Facilitated by increasing temperature, the isolated anion and its vacancy become mobile at *T >* 800 K. At very high temperatures 1000 K or above, all the YSZ crystals exhibit exceptionally high values of ionic conductivity (*σ*) and reach the order of 0*.*1 S cm-1 (Fig. 7), consistent with the range of reported ionic conductivity in experiment (Hull, 2004; Nakamura et al., 1986) at *T*  1000–1250 K. As an anion-deficient fluorite, the oxygen ion and its vacancy migration kinetics in YSZ crystal, however, cannot be described solely by a simple vacancy assisted diffusion model (Hull, 2004; Krishnamurthy et al., 2004). At elevated temperature, the vacancy diffusion rate is not always simply proportional to fractional vacancy concentration, *Di x* (with *x* defining Y2O3 dopant concentration) in the Zr1−*<sup>x</sup>*Y*x*O(2-*x/*2) system. The MSD of oxygen ions cannot increase proportional to molar concentration of the mol% Y2O3 because Y2O3 is an insulator. An optimal concentration of Y2O3 dopant exists and it is not predicted by a simple vacancy assisted diffusion mechanism. In many cases, maximum ionic conductivity occurs between ~ 8 and 15.0 mol% Y2O3 (Fig. 8) (Casselton, 1970; Hull, 2004; Krishnamurthy et al., 2004). Interestingly, this trend in ionic conductivity/diffusivity is also commonly observed in a wide range of oxide materials having the same fluorite structure.

Several investigations have addressed the origin of this unusual trend in ionic conductivity/diffusivity. Casselton, et al., (Casselton, 1970) attributes this conductivity anomaly to Coulomb attraction between dopant cations (i.e. Y3+) and oxygen vacancies. They propose that at higher dopant (and vacancy) concentrations, vacancies will be trapped in defect complexes, thus resulting in a decrease in the conductivity. Experimental studies of dopant–vacancy interactions were not unambiguous. X-ray and neutron scattering studies (Morinaga et al., 1979, 1980; Steele et al., 1974) indicate that the

amorphous thin films (Heiroth et al., 2008). Fitting to the Arrhenius relationship over the temperature range, a pre-exponential *D*0 for the 8-YSZa of ~3*.*56±0*.*29×10-5 cm2s-1 is obtained, whereas ~5*.*83±0*.*47×10-5 cm2s-1 was found for crystalline 8-YSZc (Lau et al., 2011). In general, the pre-exponential *D*0 of amorphous YSZ solids are found to be slightly less compared to their respective perfect crystals. This might be attributed to their distinct local structural properties, which subsequently affect the underlying random walks of the mobile O2− ions. Specifically, the Zr4+ and Y3+ ions are found to be lower in coordination number to the O2<sup>−</sup> ions in the irregular structural geometries (Sect. 3.1.2.II). Thus the effective probability of an oxygen ion finding a vacant site, relative to these reduced cation coordination numbers, in the amorphous lattices is increased relative to a highly coordinated, compact perfect cubic YSZ crystal. The effective activation energy, *ΔE*act, of these YSZ amorphous (3.0, 8.0 and 12.0 mol% Y2O3) solids is ~ 0*.*45–0*.*61 eV, comparable to YSZ crystals with the same Y2O3 concentration. For the 8-YSZa solid, the *ΔE*act is slightly smaller (i.e. 0*.*45±0*.*03 eV) compared to its crystal counterpart. The differences might be due to a slight variation in ionic transport within the lattices. For YSZ crystals, the ion diffusion is mainly dominated by thermally activated hopping processes of the mobile anions and vacancies, whereas for YSZ amorphous solids, there is an additional, non-negligible mutual diffusion that involves

**ii. Y2O3-dopant concentration dependent self diffusion and ionic conductivity of** 

As anticipated in Fig. 7, the high ionic conductivity in YSZ solids (both crystalline and amorphous) generally occur at rather elevated temperatures. Facilitated by increasing temperature, the isolated anion and its vacancy become mobile at *T >* 800 K. At very high temperatures 1000 K or above, all the YSZ crystals exhibit exceptionally high values of ionic conductivity (*σ*) and reach the order of 0*.*1 S cm-1 (Fig. 7), consistent with the range of reported ionic conductivity in experiment (Hull, 2004; Nakamura et al., 1986) at *T*  1000–1250 K. As an anion-deficient fluorite, the oxygen ion and its vacancy migration kinetics in YSZ crystal, however, cannot be described solely by a simple vacancy assisted diffusion model (Hull, 2004; Krishnamurthy et al., 2004). At elevated temperature, the vacancy diffusion rate is not always simply proportional to fractional vacancy concentration, *Di x* (with *x* defining Y2O3 dopant concentration) in the Zr1−*<sup>x</sup>*Y*x*O(2-*x/*2) system. The MSD of oxygen ions cannot increase proportional to molar concentration of the mol% Y2O3 because Y2O3 is an insulator. An optimal concentration of Y2O3 dopant exists and it is not predicted by a simple vacancy assisted diffusion mechanism. In many cases, maximum ionic conductivity occurs between ~ 8 and 15.0 mol% Y2O3 (Fig. 8) (Casselton, 1970; Hull, 2004; Krishnamurthy et al., 2004). Interestingly, this trend in ionic conductivity/diffusivity is also commonly observed in a wide range of oxide materials

Several investigations have addressed the origin of this unusual trend in ionic conductivity/diffusivity. Casselton, et al., (Casselton, 1970) attributes this conductivity anomaly to Coulomb attraction between dopant cations (i.e. Y3+) and oxygen vacancies. They propose that at higher dopant (and vacancy) concentrations, vacancies will be trapped in defect complexes, thus resulting in a decrease in the conductivity. Experimental studies of dopant–vacancy interactions were not unambiguous. X-ray and neutron scattering studies (Morinaga et al., 1979, 1980; Steele et al., 1974) indicate that the

cations and anions moving together in the expanded lattices (Sect. 3.2.3).

**O2- ions** 

having the same fluorite structure.

oxygen vacancy preferentially occupies interstitial sites adjacent to the dopant ions, implying significant dopant–vacancy attraction. However the extended X-ray absorption fine structure spectroscopy (EXAFS) and electron microscopy (Allpress et al., 1975; Catlow et al., 1986; Veal et al., 1988) studies indicates that oxygen vacancies (i.e. positively charged) have a greater preference for sitting in the first neighbor shell of Zr4+ ions than Y3+ ions. Hence, the dopant ion has a higher oxygen coordination number than the host Zr4+ ions. This argument suggests that the oxygen vacancy–Y-ion interaction is also dominated by size effects rather than a simple Coulomb interactions (note, the vacancy is positively charged and the Y3+-ion is negatively charged, relative to the host lattice, Zr4+-ion). Therefore to capture the Y2O3-dopant dependent ionic conductivity behavior, a small simulation cell might not be sufficient. In the meantime for YSZ, when an oxygen vacancy is formed, the lattice contracts in the vicinity of the vacancy. Because Zr4+ ions are smaller than Y3+ ions, Zr4+ ions can better accommodate this contraction. In this case, a simple lattice statics calculations based on a Born–Mayer interatomic potential also support this view (Dwivedi et al., 1990).

As shown in Fig. 8c, the simple interatomic potential from Table 1 can provide a consistent qualitative feature, a conductivity maximum, which is observed in experiments (Arachi et al., 1999; Badwal et al., 1992; Hull, 2009; Ioffe, 1978) and Kinetic Monte Carlo simulation in Fig. 8a (Krishnamurthy et al., 2004). For all three temperatures in the range 800–1200 K, the variation of conductivity, σ, with dopant concentration shows a maximum close to 8.0 mol% Y2O3 (or x ~0.16 in Zr(1−x)YxO(2-x/2)), consistent with finding the maximum close to the lower limit of the amount of doping required to stabilize the cubic YSZ phase (i.e. *c\**-YSZ) (Arachi et al., 1999; Hull, 2009; Nakamura et al., 1986).

For YSZ amorphous solids, a similar trend of the optimal Y2O3 concentration governing ionic conductivity is also supported by a recent MD simulation (Lau et al., 2011). Of the 3.0, 8.0 and 12.0 mol% Y2O3 systems (i.e. n-YSZa) studied; it was found that σmax occurs at ~ 8.0 mol% Y2O3 similar to the maximum found for the YSZ crystals. To account for this intriguing result, one simple, yet direct explanation can be the local Y3+ aggregation within the system (Lau et al., 2011). The diffusion pathways of an ion are essentially determined by the local interactions with the surrounding ions in YSZ at a given temperature. Due to the cooperative effects of larger attractive terms (Aij and ρij in Tables 1) in the BMB potential, larger bond length (i.e. RY−O from RDF at Fig. 6) and larger Y3+ ion relative to Zr4+ ion in the lattice, steeper potential wells among Y3+–O2− relative to Zr4+–O2− are expected. Thus, this will hinder the mobility of O2− (or its vacancy), if the local Y3+ aggregation is large enough. At a given temperature, the local Y3+ aggregation can be defined qualitatively by comparing to the degree of Y3+-clustering, *η* (Table 2), based on the thermal equilibrium structures of each system. For each of the *N*Y atoms, we count the number of other Y atoms less than 3*.*7Å away. If *Nc* is the number of these Y atoms that have four or more neighboring Y atoms, then *η* = *Nc/N*×100*.*0*%*. As shown in Table 2, *η* increased as the mol% of Y2O3 increased regardless of crystal or amorphous geometry. As the system temperature goes up, *η* generally decreases due to thermal expansion of the lattice. Despite this change with temperature, the relative comparison of each system remains the same. Therefore, this suggests that YSZ ionic conductivity is a thermally activated process as shown in Fig. 7 and 8, but the subtle interplay between the vacancy concentration and the intrinsic ionic diffusion pathways that are defined by local structures and local orderings cannot be ignored.

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 361

amorphous lattice) matrix through which the ions travel together. The charge-compensating cation sites are fixed to the matrix about which the ions are loosely bonded by the weakerthan covalent ionic bonds. Because of these ionic bonds, often the thermal energy in the solids can provide the energy needed for the ion to dissociate from its site and "hop" to an adjacent site. These phenomena lead to a so-called "activation energy" (i.e. ΔEact in Fig. 7), which defines the mean effective energy associated with the hop of an ion between adjacent sites. The rate of hopping is mostly controlled by the temperature and is given by the Boltzmann factor (or Arrhenius law) as anticipated in Fig. 7. MD simulations enable us to understand in depth this interesting mechanism. We can track the position of these ions over time, but we may need to utilize some ensemble-averaged quantities such as the meansquare displacement <*r2(t)*> (MSD) as a key to understanding these motions, and elucidating the slight differences of ionic motion between the YSZ-c (crystals) and YSZa (amorphous)

System *η* (T = 300K) *η* (T = 1000K) Crystal (mol%) 3.0 2.46 3.36 8.0 20.52 15.26 12.0 25.93 17.69 Amorphous (mol%) 3.0 1.37 0.73 8.0 16.67 12.82 12.0 22.62 24.44 Table 2. The degree of Y3+-clustering, *η* (in %), of 3.0, 8.0 and 12.0 mol% Y2O3 YSZ crystals and amorphous solids at *T* = 300 and 1000 K. (Note: only the Y3+-clusters with the size ≥ 4 ions are counted. By assuming the cluster cut-off radius [*R(*Y3+−Y3+)] to be ~ 3*.*7 Å as found by RDF, we define *η* = *Nc/N*×100*.*0*%*, where *Nc* is the number of clusters with ≥ 4 Y3+-ions

and *N* is the total number of Y3+-clusters in a system adopted from Lau et al.

As we suspected, diffusion in atomistic scale in a lattice can be mainly determined by the internal structure of a material. Thus the elements of the underlying random walks are often biased on the detailed local microstructure of the material, and the strength of the interatomic forces. In the high-temperature regime, the MSD corresponding to the diffusion of the dominant charge carriers (i.e. oxygen ions O2−) will grow linearly in time, analogous to ionic transport in a liquid (Fig. 9). This linear diffusive regime typically can also be found in other reported literature, which used interatomic potentials (Bush et al., 1994; Devanathan et al., 2006; Dwivedi et al., 1990; Fisher et al., 1998, 1999; Khan et al., 1998; Kilo et al., 2003; Lau et al., 2011; Lewis et al., 1985; Li et al., 1995; Minervini et al., 2000; Okazaki et al., 1994; Sawaguchi et al., 2000; Schelling et al., 2001; Shimojo et al., 1992; ; van Duin et al., 2008; Yamamura et al., 1999; Zacate et al., 2000) different from those of Table 1. For the much heavier cations (Zr4+ and Y3+) in YSZ crystals, the MSD generally decays and saturates over a period of time, as the kinetic energy is not sufficient to reach monotonic diffusive behavior, and therefore the path is bounded in the fixed crystal lattice (Fig. 9). This unique feature, however, is not found in YSZ amorphous solids. In the

**i. Information given by the mean-square displacement** 

structures as we found in Fig. 7.

(Lau et al., 2011).

Fig. 8. (a) Oxygen diffusivity as a function of the Y2O3 mole fraction, *y*, and yttrium ion concentration, *x*, in YSZ at different temperatures, adopted from Krishnamurthy et. al. (Krishnamurthy et al., 2004). (b) The variation of the ionic conductivity *σ*i of Zr1−*<sup>x</sup>*Y*x*O(2-*x/*2) with dopant concentration *x* at three different temperatures. The maximum in *σ*i occurs close to the lower limit of stability of the stabilized cubic form (*c\**), adopted from Hull, 2009. (c) The predicted ionic conductivity, *σ* (in Scm−1), of 3.0–12.0 mol*%* Y2O3 YSZ crystal over 800– 1200 K, with the conductivity maxima as a function of Y2O3 concentration based on the simple interatomic potential (Lau et al., 2011) of Table 1.

#### **3.2.2 Why can ionic transport in YSZ crystal and YSZ amorphous be different?**

As we have seen from previous sections, the basic static structural features and dynamics of ions in both YSZ crystals (*n*-YSZc) and amorphous solids (*n*-YSZa) can be well-described through the simple BMB potentials (Table 1). For both systems, the key-elements are almost the same: the system consists of a rigid ordered (i.e. crystalline lattice) or disordered (i.e.

Fig. 8. (a) Oxygen diffusivity as a function of the Y2O3 mole fraction, *y*, and yttrium ion concentration, *x*, in YSZ at different temperatures, adopted from Krishnamurthy et. al. (Krishnamurthy et al., 2004). (b) The variation of the ionic conductivity *σ*i of Zr1−*<sup>x</sup>*Y*x*O(2-*x/*2) with dopant concentration *x* at three different temperatures. The maximum in *σ*i occurs close to the lower limit of stability of the stabilized cubic form (*c\**), adopted from Hull, 2009. (c) The predicted ionic conductivity, *σ* (in Scm−1), of 3.0–12.0 mol*%* Y2O3 YSZ crystal over 800– 1200 K, with the conductivity maxima as a function of Y2O3 concentration based on the

**3.2.2 Why can ionic transport in YSZ crystal and YSZ amorphous be different?** 

As we have seen from previous sections, the basic static structural features and dynamics of ions in both YSZ crystals (*n*-YSZc) and amorphous solids (*n*-YSZa) can be well-described through the simple BMB potentials (Table 1). For both systems, the key-elements are almost the same: the system consists of a rigid ordered (i.e. crystalline lattice) or disordered (i.e.

simple interatomic potential (Lau et al., 2011) of Table 1.

amorphous lattice) matrix through which the ions travel together. The charge-compensating cation sites are fixed to the matrix about which the ions are loosely bonded by the weakerthan covalent ionic bonds. Because of these ionic bonds, often the thermal energy in the solids can provide the energy needed for the ion to dissociate from its site and "hop" to an adjacent site. These phenomena lead to a so-called "activation energy" (i.e. ΔEact in Fig. 7), which defines the mean effective energy associated with the hop of an ion between adjacent sites. The rate of hopping is mostly controlled by the temperature and is given by the Boltzmann factor (or Arrhenius law) as anticipated in Fig. 7. MD simulations enable us to understand in depth this interesting mechanism. We can track the position of these ions over time, but we may need to utilize some ensemble-averaged quantities such as the meansquare displacement <*r2(t)*> (MSD) as a key to understanding these motions, and elucidating the slight differences of ionic motion between the YSZ-c (crystals) and YSZa (amorphous) structures as we found in Fig. 7.


Table 2. The degree of Y3+-clustering, *η* (in %), of 3.0, 8.0 and 12.0 mol% Y2O3 YSZ crystals and amorphous solids at *T* = 300 and 1000 K. (Note: only the Y3+-clusters with the size ≥ 4 ions are counted. By assuming the cluster cut-off radius [*R(*Y3+−Y3+)] to be ~ 3*.*7 Å as found by RDF, we define *η* = *Nc/N*×100*.*0*%*, where *Nc* is the number of clusters with ≥ 4 Y3+-ions and *N* is the total number of Y3+-clusters in a system adopted from Lau et al. (Lau et al., 2011).

#### **i. Information given by the mean-square displacement**

As we suspected, diffusion in atomistic scale in a lattice can be mainly determined by the internal structure of a material. Thus the elements of the underlying random walks are often biased on the detailed local microstructure of the material, and the strength of the interatomic forces. In the high-temperature regime, the MSD corresponding to the diffusion of the dominant charge carriers (i.e. oxygen ions O2−) will grow linearly in time, analogous to ionic transport in a liquid (Fig. 9). This linear diffusive regime typically can also be found in other reported literature, which used interatomic potentials (Bush et al., 1994; Devanathan et al., 2006; Dwivedi et al., 1990; Fisher et al., 1998, 1999; Khan et al., 1998; Kilo et al., 2003; Lau et al., 2011; Lewis et al., 1985; Li et al., 1995; Minervini et al., 2000; Okazaki et al., 1994; Sawaguchi et al., 2000; Schelling et al., 2001; Shimojo et al., 1992; ; van Duin et al., 2008; Yamamura et al., 1999; Zacate et al., 2000) different from those of Table 1. For the much heavier cations (Zr4+ and Y3+) in YSZ crystals, the MSD generally decays and saturates over a period of time, as the kinetic energy is not sufficient to reach monotonic diffusive behavior, and therefore the path is bounded in the fixed crystal lattice (Fig. 9). This unique feature, however, is not found in YSZ amorphous solids. In the

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 363

amorphous structure from MD trajectories, the three constituent ions (Zr4+, Y3+, O2−) have similar slopes over a long period of simulation time, indicating that mutual diffusion driven by correlated motion between anions and cations might be substantial in the YSZ

In this case, the difference can become more evident if temperature is raised. Here it is noteworthy to point out a limit due to finite computer resources. This subtle phenomenon can also attributed to the fact that the simulated YSZ amorphous systems are mostly metastable and are continuously crystallizing at a slow rate even at elevated temperatures, analogous to crystallization of amorphous YSZ film which was reported in a recent experiment (Heiroth et al., 2008). Unlike that of a perfect crystal, the potential-energy landscape experienced by an ion in a disordered amorphous solid is irregular, containing a distribution of stationary points, wells, barrier heights, and saddle point energies. Thus, the residence sites for an ion in an amorphous system are more vulnerable to thermal fluctuations, and possibly are even metastable. Thus, a much longer simulation time might be needed to achieve thermal equilibrium for an amorphous YSZ structure compared to YSZ crystals as pointed out in a recent study (Lau et al., 2011). In real experiments, the stability of amorphous YSZ often depends on the crystallization process and is directly correlated to grain size and different local environments (Heiroth et al., 2008). This suggests that besides short-distance thermal vibrations, the dynamics of both the cation and anion can be responsible for the observed slow crystallization process and related ionic conductivity

To capture the intriguing observation as shown in MSD (Fig. 9), the physics behind the timedependent mutual diffusion and structural variations observed particularly in YSZ amorphous structure can be analyzed based on the van Hove correlation function, *G(r,t)* (Frenkel, 1996; Lau, 2011, Rapaport, 2004). By following the definition (Frenkel & Smit, 1996;

Gs(r,t) = <Σjδ(r+rj(0)-rj(t))>/Nj (3)

 Gd(r,t) = <Σi≠jδ(r+ri(0)-rj(t))>/Ni (4) which are determined by the space- and time-dependent density correlations of the constituent ions in the system. As defined, *Gs(r,t)* compares the position of a particle to its position within a time interval. In thermal equilibrium, *Gs(r,t)* depends only on the time difference, and provides information about the hopping and diffusion processes of the separated ions. Whereas *Gd(r,t)* compares the positions of a particle to the position of another particle at different time, and yields information about the correlated motions

As expected for a diffusion process, *Gs(r,t)* decreases generally with increasing distance at each time and reflects the nature of Brownian motion, which is proportional to the system temperature. For a mobile ion, the width of *Gs(r,t)* generally are broader than the stationary ion. At low temperature, the *Gs(r,t)* are generally found to be more localized, with longer residence time for the individual ions. Thus for a YSZ crystal (YSZc) and a YSZ amorphous

amorphous structure.

and the distinct part:

within a time interval.

fluctuation seen in experiments (Heiroth et al., 2008).

**ii. Information given by van hove correlation functions** 

Rapaport, 2004), *G(r,t)* can be represented separately by the self:

YSZ amorphous phase, the heavy cations are not simply bound to a rigid crystal lattice as in YSZ crystals (Fig. 9).

Fig. 9. The MSD of different ions (Zr4+, Y3+, O2-) in YSZ system (crystal and amorphous) that were obtained from MD simulation through the averaging over the ensemble of the trajectories of respective ions. The systems shown are 3-YSZc and 3-YSZa at 800 K for 2.5 ns trajectories from Lau et al. (Lau et al., 2011). Representative superposed pictures of the different motions are shown to the upper left (van Duin et al., 2008).

Instead of fixed in the rigid "stationary" lattices, substantial motions can be observed through the MSD. Comparable to the negligible motion of the cations in YSZ crystals (i.e. <*r*2> ~ 0*.*2 Å2), the displacements of Zr4+ and Y3+ ions are rather substantial (i.e. <*r2>*~ 0*.*5 Å2 at 800 K within 2.5 ns), are almost as diffusive (Fig. 9) as the mobile O2− anions. For YSZ amorphous structure from MD trajectories, the three constituent ions (Zr4+, Y3+, O2−) have similar slopes over a long period of simulation time, indicating that mutual diffusion driven by correlated motion between anions and cations might be substantial in the YSZ amorphous structure.

In this case, the difference can become more evident if temperature is raised. Here it is noteworthy to point out a limit due to finite computer resources. This subtle phenomenon can also attributed to the fact that the simulated YSZ amorphous systems are mostly metastable and are continuously crystallizing at a slow rate even at elevated temperatures, analogous to crystallization of amorphous YSZ film which was reported in a recent experiment (Heiroth et al., 2008). Unlike that of a perfect crystal, the potential-energy landscape experienced by an ion in a disordered amorphous solid is irregular, containing a distribution of stationary points, wells, barrier heights, and saddle point energies. Thus, the residence sites for an ion in an amorphous system are more vulnerable to thermal fluctuations, and possibly are even metastable. Thus, a much longer simulation time might be needed to achieve thermal equilibrium for an amorphous YSZ structure compared to YSZ crystals as pointed out in a recent study (Lau et al., 2011). In real experiments, the stability of amorphous YSZ often depends on the crystallization process and is directly correlated to grain size and different local environments (Heiroth et al., 2008). This suggests that besides short-distance thermal vibrations, the dynamics of both the cation and anion can be responsible for the observed slow crystallization process and related ionic conductivity fluctuation seen in experiments (Heiroth et al., 2008).

#### **ii. Information given by van hove correlation functions**

To capture the intriguing observation as shown in MSD (Fig. 9), the physics behind the timedependent mutual diffusion and structural variations observed particularly in YSZ amorphous structure can be analyzed based on the van Hove correlation function, *G(r,t)* (Frenkel, 1996; Lau, 2011, Rapaport, 2004). By following the definition (Frenkel & Smit, 1996; Rapaport, 2004), *G(r,t)* can be represented separately by the self:

$$\mathbf{G}\_{\mathbf{s}}(\mathbf{r}, \mathbf{t}) = \, \! \Sigma\_{\mathbf{i}} \delta(\mathbf{r} + \mathbf{r}\_{\parallel}(0) \cdot \mathbf{r}\_{\parallel}(\mathbf{t})) \rhd / \mathcal{N}\_{\mathbf{i}} \tag{3}$$

and the distinct part:

362 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

YSZ amorphous phase, the heavy cations are not simply bound to a rigid crystal lattice as

Fig. 9. The MSD of different ions (Zr4+, Y3+, O2-) in YSZ system (crystal and amorphous) that

trajectories of respective ions. The systems shown are 3-YSZc and 3-YSZa at 800 K for 2.5 ns trajectories from Lau et al. (Lau et al., 2011). Representative superposed pictures of the

Instead of fixed in the rigid "stationary" lattices, substantial motions can be observed through the MSD. Comparable to the negligible motion of the cations in YSZ crystals (i.e. <*r*2> ~ 0*.*2 Å2), the displacements of Zr4+ and Y3+ ions are rather substantial (i.e. <*r2>*~ 0*.*5 Å2 at 800 K within 2.5 ns), are almost as diffusive (Fig. 9) as the mobile O2− anions. For YSZ

were obtained from MD simulation through the averaging over the ensemble of the

different motions are shown to the upper left (van Duin et al., 2008).

in YSZ crystals (Fig. 9).

$$\mathbf{G}\_{\rm d}(\mathbf{r}, \mathbf{t}) = \sphericalangle \mathbf{G}(\mathbf{r} + \mathbf{r}\_{\rm i}(\mathbf{0}) \cdot \mathbf{r}\_{\rm \parallel}(\mathbf{t})) \gtrless / \mathbf{N}\_{\rm i} \tag{4}$$

which are determined by the space- and time-dependent density correlations of the constituent ions in the system. As defined, *Gs(r,t)* compares the position of a particle to its position within a time interval. In thermal equilibrium, *Gs(r,t)* depends only on the time difference, and provides information about the hopping and diffusion processes of the separated ions. Whereas *Gd(r,t)* compares the positions of a particle to the position of another particle at different time, and yields information about the correlated motions within a time interval.

As expected for a diffusion process, *Gs(r,t)* decreases generally with increasing distance at each time and reflects the nature of Brownian motion, which is proportional to the system temperature. For a mobile ion, the width of *Gs(r,t)* generally are broader than the stationary ion. At low temperature, the *Gs(r,t)* are generally found to be more localized, with longer residence time for the individual ions. Thus for a YSZ crystal (YSZc) and a YSZ amorphous

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 365

Fig. 10. The self-part of the van Hove correlation function 4π*r*2*Gs(r, t)* (lines with symbols) for both 3-YSZc (crystal) and 3-YSZa (amorphous) at 800 K. The functions develop from top to bottom with time, and are compared with Gaussian diffusion motion 4π*r*2*C(r,t)* with *C(r,t)*  (plane lines) determined by the steady-state diffusion constant *D* at 10 ps (a) and at 2500 ps (b) for crystal (left) and amorphous solid (center). (c) A schematic figure that illustrates a time-evolving ion's hopping (large fluctuations) and diffusion (small fluctuations) processes

in both YSZ crystal and amorphous solid that can be captured through the van Hove

amorphous solids (e.g. 3-YSZa in Fig. 10), the coexistence of the 'liquid-like' features of diffusive cations (Zr4+, Y3+) and anions (O2−) however is more obvious, with the long dispersive tails of *Gs(r,t)* spreading to large *r*. Instead of the localization of the maximum at the origin as the initial stage, substantial drifting (Δ*r ~* 0*.*4–0*.*7 Å) of the maximum for both cations and anions in YSZ amorphous solids after 2500 ps can be seen (Fig. 10b). This suggests that the residence sites in a metastable disordered amorphous solid of YSZ are

correlation function.

(YSZa) structure, their differences in MSD over a simulation time can also be highlighted in the van Hove correlation as shown in Fig. 10. To capture the time-evolution of ions' hopping and diffusion process in both the YSZ crystals and amorphous structures, a direct comparison of *Gs(r,t)* with a Gaussian form of diffusive motion can be very useful (Lau et al., 2011).

For a diffusion process, Fick's laws apply when there is a concentration gradient in a system. However as a typical thermodynamic equilibrium system, these MD calculations have no long-range gradient. To better analyze the features in *Gs (r,t)*, a comparison of *Gs(r,t)* with the diffusion transport character in Gaussian form is particularly useful. By applying the Laplace transform to Fick's second law (*D2C* =*∂C/∂t*), subject to the initial condition at time *C(x, y, z, t* = 0*)* for all *x, y, z* = 0 and the boundary condition *C(*±∞*, t)* = 0 so that a finite number of diffusive ions cannot alter the charge composition of the system, the diffusion process is subject to the diffusion solution in Gaussian form, with concentration field

$$\mathbf{C}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) = \exp(\cdot (\mathbf{x}^2 + \mathbf{y}^2 + \mathbf{z}^2) / 4 \mathbf{D} \mathbf{t}) / \langle 4 \mathbf{n} \mathbf{D} \mathbf{t} \rangle^{3/2} = \exp(\cdot \mathbf{r}^2 / 4 \mathbf{D} \mathbf{t}) / \langle 4 \mathbf{n} \mathbf{D} \mathbf{t} \rangle^{3/2} \tag{5}$$

where *r2 = x2+y2+z2*. Assuming that the time response of the van Hove self-correlation function *Gs(r, t)* of ions in the YSZ crystal and amorphous solid (i.e. 3-YSZc and 3-YSZa in Fig. 10) is reminiscent of *C(x, y, z, t)* [or *C(r, t)*], the evolution of the self-part of the van Hove function 4*πr*2*Gs(r, t)* as a function of *r* with time (i.e. 10 and 2500 ps shown in Fig. 10) at 800 K can be compared with an ideal Gaussian distribution of diffusive motion of ions, 4π*r*2*C(r, t)*  with *C(r, t)* = exp*(*−*r* 2*/*4*Dt)/(*4*π Dt)*3*/*2 where <*r*2> ~ 6*Dt* is assumed and *D* is the diffusion coefficient of the corresponding ion.

As shown by *C(r, t)*, its normal distribution is essentially governed by the ionic motion under thermal equilibrium at a given temperature that is analogous to a 'drifting wavepacket' under the thermal motion. Basically, the drifting 'wave-packet' of charged carriers is determined by the mobility of the ions in the ensemble, which can be captured by <*r2>*. As time develops, the Gaussian-shaped wave-packet widens as it travels due to the ion diffusion about the mean (Fig. 10). Compared to the Gaussian form diffusive motion, 4π*r*<sup>2</sup> *C(r, t)* in Fig. 10, the decay of the profile and variation found in 4π*r*2*Gs (r, t)* of the YSZ system is prominent as time develops. From its initial localized 'impulse-like' distribution (e.g. 4 π*r*2*Gs(r, t)* at 10 ps in Fig. 10a) analogous to a localized *δ*-function (Frenkel & Smit, 1996; Rapaport, 2004), the variation of *Gs(r, t)* for the ions is obviously less uniform compared to *C(r, t)* in Gaussian form. Here neither of the ionic motions in the two YSZ systems is fitted perfectly by a Gaussian distribution, and strong deviation from the Gaussian form is prominent as time develops (Fig. 10). In this case, the both YSZ crystals and amorphous solid share similar features in *Gs (r, t)* up to 2500 *ps*, yet are not perfectly identical. In both cases, the non-Gaussian features of the mobile O2− in *Gs(r,t)* can be attributed to a non- uniform distribution of the mobility of the carriers which varies over a wide range of velocity caused by the coexistence of hopping processes and 'liquid-like' diffusivity. As time evolves, the multiple satellite peaks which represent the hopping processes are evident in both YSZ systems for all the constituent ions. For YSZ crystal (e.g. 3- YSZc in Fig. 10), we can see only the mobile anions attain the significant diffusive features, namely a decaying peak accompanied by a tail part of the *Gs(r, t)* function that develops with time. The long tail at large *r* is related to the long-range motion found in the timedependent trajectory, i.e. by fast mobile anions with successive hopping motions. For YSZ

(YSZa) structure, their differences in MSD over a simulation time can also be highlighted in the van Hove correlation as shown in Fig. 10. To capture the time-evolution of ions' hopping and diffusion process in both the YSZ crystals and amorphous structures, a direct comparison of *Gs(r,t)* with a Gaussian form of diffusive motion can be very useful

For a diffusion process, Fick's laws apply when there is a concentration gradient in a system. However as a typical thermodynamic equilibrium system, these MD calculations have no long-range gradient. To better analyze the features in *Gs (r,t)*, a comparison of *Gs(r,t)* with the diffusion transport character in Gaussian form is particularly useful. By applying the

*C(x, y, z, t* = 0*)* for all *x, y, z* = 0 and the boundary condition *C(*±∞*, t)* = 0 so that a finite number of diffusive ions cannot alter the charge composition of the system, the diffusion

 C(x,y,z,t) = exp(-(x2+y2+z2)/4Dt)/(4πDt)3/2= exp(-r2/4Dt)/(4πDt)3/2 (5) where *r2 = x2+y2+z2*. Assuming that the time response of the van Hove self-correlation function *Gs(r, t)* of ions in the YSZ crystal and amorphous solid (i.e. 3-YSZc and 3-YSZa in Fig. 10) is reminiscent of *C(x, y, z, t)* [or *C(r, t)*], the evolution of the self-part of the van Hove function 4*πr*2*Gs(r, t)* as a function of *r* with time (i.e. 10 and 2500 ps shown in Fig. 10) at 800 K can be compared with an ideal Gaussian distribution of diffusive motion of ions, 4π*r*2*C(r, t)*  with *C(r, t)* = exp*(*−*r* 2*/*4*Dt)/(*4*π Dt)*3*/*2 where <*r*2> ~ 6*Dt* is assumed and *D* is the diffusion

As shown by *C(r, t)*, its normal distribution is essentially governed by the ionic motion under thermal equilibrium at a given temperature that is analogous to a 'drifting wavepacket' under the thermal motion. Basically, the drifting 'wave-packet' of charged carriers is determined by the mobility of the ions in the ensemble, which can be captured by <*r2>*. As time develops, the Gaussian-shaped wave-packet widens as it travels due to the ion diffusion about the mean (Fig. 10). Compared to the Gaussian form diffusive motion, 4π*r*<sup>2</sup> *C(r, t)* in Fig. 10, the decay of the profile and variation found in 4π*r*2*Gs (r, t)* of the YSZ system is prominent as time develops. From its initial localized 'impulse-like' distribution (e.g. 4 π*r*2*Gs(r, t)* at 10 ps in Fig. 10a) analogous to a localized *δ*-function (Frenkel & Smit, 1996; Rapaport, 2004), the variation of *Gs(r, t)* for the ions is obviously less uniform compared to *C(r, t)* in Gaussian form. Here neither of the ionic motions in the two YSZ systems is fitted perfectly by a Gaussian distribution, and strong deviation from the Gaussian form is prominent as time develops (Fig. 10). In this case, the both YSZ crystals and amorphous solid share similar features in *Gs (r, t)* up to 2500 *ps*, yet are not perfectly identical. In both cases, the non-Gaussian features of the mobile O2− in *Gs(r,t)* can be attributed to a non- uniform distribution of the mobility of the carriers which varies over a wide range of velocity caused by the coexistence of hopping processes and 'liquid-like' diffusivity. As time evolves, the multiple satellite peaks which represent the hopping processes are evident in both YSZ systems for all the constituent ions. For YSZ crystal (e.g. 3- YSZc in Fig. 10), we can see only the mobile anions attain the significant diffusive features, namely a decaying peak accompanied by a tail part of the *Gs(r, t)* function that develops with time. The long tail at large *r* is related to the long-range motion found in the timedependent trajectory, i.e. by fast mobile anions with successive hopping motions. For YSZ

*2C* =*∂C/∂t*), subject to the initial condition at time

process is subject to the diffusion solution in Gaussian form, with concentration field

(Lau et al., 2011).

Laplace transform to Fick's second law (*D*

coefficient of the corresponding ion.

Fig. 10. The self-part of the van Hove correlation function 4π*r*2*Gs(r, t)* (lines with symbols) for both 3-YSZc (crystal) and 3-YSZa (amorphous) at 800 K. The functions develop from top to bottom with time, and are compared with Gaussian diffusion motion 4π*r*2*C(r,t)* with *C(r,t)*  (plane lines) determined by the steady-state diffusion constant *D* at 10 ps (a) and at 2500 ps (b) for crystal (left) and amorphous solid (center). (c) A schematic figure that illustrates a time-evolving ion's hopping (large fluctuations) and diffusion (small fluctuations) processes in both YSZ crystal and amorphous solid that can be captured through the van Hove correlation function.

amorphous solids (e.g. 3-YSZa in Fig. 10), the coexistence of the 'liquid-like' features of diffusive cations (Zr4+, Y3+) and anions (O2−) however is more obvious, with the long dispersive tails of *Gs(r,t)* spreading to large *r*. Instead of the localization of the maximum at the origin as the initial stage, substantial drifting (Δ*r ~* 0*.*4–0*.*7 Å) of the maximum for both cations and anions in YSZ amorphous solids after 2500 ps can be seen (Fig. 10b). This suggests that the residence sites in a metastable disordered amorphous solid of YSZ are

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 367

amorphous solids, the mobile anions and the slowly migrating cations are strongly coupled in their motion. This mutual diffusion (cations and anions) found in the amorphous phases contributes considerably to the dynamics as seen through the timedependent van Hove correlation functions. This reduces the effective oxygen-ion conductivity. These moving ions carry charges, and thus produce an electrical response. It is expected that the intriguing features of mutual diffusion found in amorphous YSZ solids can therefore be detected by current experimental techniques at frequencies below the typical vibrational frequencies (>100 GHz). To gain further insight into the different correlated motions in various YSZ system and its interfaces, richness in morphologies, longer equilibration, better statistical methods and better atomic potentials are highly

This work was supported and by the Office of Naval Research, both directly and through

Ackermann, R J; Rauh, E G. & Alexander, C A. (1975) High Temp. Sci. Vol. 7, pp. 305.

Andersson, M.; Yuan, J. & Sundén, B. (2010) Applied Energy, Vol. 87, pp. 1461-1476.

Arachi, Y; Sakai, H; Yammoto, O; Takeda, Y. & Imanishai, N. (1999) Solid State Ion. Vol. 121,

Baller, J; Krüger J K; Birringer, R. & Prousi, C. (2000) J. Phys.:Condens. Matter Vol. 12, pp.

Bogicevic, A; Wolverton, C; Crosbie, G M. & Stechel E B. (2001) Phys. Rev. B Vol. 64, pp.

Butz, B.; Kruse, P.; Störmer, H.; Gerthsen, D.; Müller, A.; Weber, A. & Ivers-Tiffée, E. (2006)

Bush, T.S.; Gale, J.D.; Catlow, C.R.A. & Battle, P.D. (1994) J. Mater. Chem. Vol. 4, pp.

Catlow, C. R. A.; Chadwick, A. V.; Greaves, G. N. & Moroney, L. M. J. Am. Ceram. Soc. Vol.

Chen, X. J.; Khor, K A.; Chan, S H. & Yu, L G. (2002) Mater. Sci. Eng. A, Vol. 335, pp.

Cheng, Z.; Wang, J-H.; Choi, Y.; Yang, L; Lin, M.C. & Liu, M. (2011) Energy Environ. Sci.

Dash, L K; Vast, N; Baranek, P; Cheynet M C. & Reining, L. (2004) Phys. Rev. B Vol. 70, pp.

Aldebert, P. & Traverse, J P. (1985) J. Am. Ceram. Soc. Vol. 68, pp. 34. Allpress, J. G. & Rossell, H. J. J. Solid State Chem. Vol. 15, pp. 68–78.

Boysen, H; Frey, F. & Vogt, T. (1991) Acta Crystallogr. B Vol. 47, pp. 881.

Casselton, R. E. W. (1970) Phys. Status Solidi A Vol. 2, pp. 571–585.

Chu, W-F; Thangadurai, V. & Weppner, W. (2006). Ionics, Vol. 12, pp. 1-6.

Badwal, S. P. S. (1992) Solid State Ion. Vol. 52, pp. 23.

Solid State Ion. Vol. 177, pp. 3275.

Vol. DOI: 10.1039/c1ee01758f.

desirable.

**6. References** 

**5. Acknowledgements** 

pp. 133.

5403.

014106.

831.

246.

245116.

69, pp. 272–277 .

the Naval Research Laboratory.

always influenced by interactions with the migrating cations, and therefore change with time. By involving the mutual diffusion (i.e. sites previously occupied by cations can be visited by anions and vice versa), the hopping processes of ions will therefore be influenced greatly by the changing intermediate surroundings, which makes fast diffusive ions with successive hopping jumps less probable. Therefore the MSD of the mobile O2− in the amorphous solid will be less compared to the crystal (Fig. 9).
