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

Solid ionic compounds normally form insulating compounds because each ion usually has a closed-shell electronic configuration, which prevents delocalization of electrons and thus electronic conduction in the system. In an appropriate temperature regime, charge can be transported in these solid materials, however, by the motion of highly mobile ionic species (e.g. Li+, Na+, O2-, etc.) (Fig. 1 Top) (Chu et al., 2006). This phenomenon is termed “ionic conductivity”. Such ionic compounds exhibit liquid-like conductivities whilst still in the solid state, i.e. at temperatures are well below their melting points (Hull, 2004). Exploiting this unusual property, considerable progress has been made in recent years in fuel cells and batteries, which are promising key technologies to meet our rising energy and environmental needs.

(Top) Temperature dependence of the best solid ionic conductors of Ag + , Na + , Li + , and O 2− ions (adapted from Chu et al., 2006). The ionic conductivities fall within ranges which are highlighted by the shaded areas. (Middle) A simple atomistic model of SOFC that consists of three basic components: anode, cathode and electrolyte. (Bottom) The basic properties of solid electrolyte of SOFC can be simulated using standard atomistic modeling. al., 2010). In short, the overall electrochemical cell reaction is built upon the chemical potential difference between the cathode and anode sides, which is given by the Nernst equation in the simplest approach.
To maximize this superionic conductivity in a SOFC, alloying zirconia (ZrO 2 ) with various metal oxides (e.g. Y 2 O 3 , Sc 2 O 3 , La 2 O 3 , MgO, CaO, etc) is a plausible approach (Lashtabeg et al., 2006). In certain regimes, alloying can stabilize the highly symmetric cubic-fluorite phase of zirconia, which consequently facilitates the ionic conductivity through the introduction of oxygen vacancies as the host Zr 4+ cations are replaced by dopant aliovalent cations. Above all, yttria-stabilized zirconia, YSZ (i.e. the Zr 1−x Y x O 2−x/2 system, with x/2 being the Y 2 O 3 dopant concentration), is the most common choice for the electrolyte (Fig. 1 Top), due to its good oxide ion conductivity over a wide range of oxygen partial pressures, its stability under oxidizing and reducing conditions, and its good high-temperature mechanical properties (Ralph et al., 2001). This versatility ultimately arises from atomic defects (i.e. oxide ion vacancies) in the cubic zirconia crystalline lattices. The vacancies, through the coupled interactions among the vacancies and ions in these different alloys can dramatically affect the structural, thermal, mechanical and electrical properties of the system. Furthermore, the optimum ionic conductivity of each alloy varies with synthesis route and sintering conditions due to the resultant diverse local morphologies and microstructures (Badwal et al., 1992;Butz et al., 2006;Chen et al., 2002;Fukui et al., 2004;Ioffe et al., 1978;Korte et al., 2008;Zhang et al., 2007;Zhu et al., 2005). Thus a variety of different morphologies and microstructures of YSZ can be found in experiments (Badwal et al., 1992;Butz et al., 2006;Chen et al., 2002;Cheng et al., 2011;Etsell et al., 1970;Fukui et al., 2004;Ioffe et al., 1978;Korte et al., 2008;Lashtabeg et al., 2006;Zhang et a., 2007;Zhu et al., 2005). Non-cubic crystalline phases, grain boundaries, and disordered lattices of amorphous features can commonly be found. To understand how the local microstructures and system morphologies affect the ionic conductivity and degradation of the electrolyte (YSZ) upon SOFC operation, theoretical simulation can provide detailed in situ atomistic information that is difficult to obtain experimentally.
To increase the basic understanding of the scientific phenomena that underlie current experimental findings, which could most dramatically affect engineering design, atomistic modeling based on quantum mechanics (e.g. first-principles or ab initio methods), molecular dynamics and Monte Carlo simulation (Fig. 2) are particularly useful. They provide relevant predictions of crystal structure, energetics, and vibrational frequencies through detailed atomistic descriptions that complement the standard analytical continuum model of ion diffusion at the macroscopic scale (Andersson et al., 2010;Cheng et al., 2011). Furthermore, various fundamental physical processes and chemical reactions can be described via quantum mechanics. The most accurate quantum mechanical atomistic descriptions, firstprinciples or ab initio methods, however are limited to small system size (less than about 1000 atoms) and short times (less than a few nanoseconds). For long-time-scales and larger length scales in atomistic modeling, recent sophisticated studies use kinetic Monte Carlo (KMC) (Gatewood et al., 2011;Lau et al., 2008Pornprasertsuk et al., 2007;Turner et al., 2010;Wang et al., , 2011 simulations. KMC can treat the longest time scales for chemical reactions, because reactants are not followed to products in time. Instead, chemical reactions are stochastically chosen to occur at the assigned rates. This is important for SOFC modeling, because charge-separating electrochemical reactions tend to have slower rates. KMC has modeled the ionic transport in a SOFC under various operating conditions  (Gatewood et al., 2011;Lau et al., 2008Pornprasertsuk et al., 2007;Turner et al., 2010;Wang et al., , 2011. Specifically, KMC can probe SOFC performance by simultaneously capturing various reaction pathways of electrochemical and physicochemical reactions in the electrolyte and at the three-phase boundary (TPB), i.e. the interface where the gas reactants, electrolyte, and electrode meet. The electrical current through the YSZ is simulated in direct current (dc) and alternating current (ac), as electrochemical impedance spectra under various operating conditions using a minimal set of uniform chemical reaction rates on an assumed cubic YSZ lattice via KMC (Gatewood et al., 2011;Lau et al., 2008Pornprasertsuk et al., 2007;Turner et al., 2010;Wang et al., , 2011. Despite the robustness of KMC simulation, this approach is based on a rigid lattice gas model and can not predict the experimentally observed ionic conductivity maximum as a function of Y 2 O 3 dopant concentration (Hull, 2009). Beyond the kinetics driven atomic motion as implemented in minimal KMC models, the complex dynamics of real lattices and the realtime multi-particle ion-vacancy interactions at a finite temperature can be computed 'on-thefly' in simulation. However simple or complex, the ionic conductivity of any atomistic model of solid YSZ is, the conductivity can accurately be derived from standard equilibrium classical MD simulation (Frenkel & Smit, 1996). Thus to explore and understand the nature of unique ionic conductivity in the solid electrolytes of SOFC's, YSZ solids have been intensively studied using classical MD simulations (Devanathan et al., 2006;Fisher et al., 1998Fisher et al., , 1999Khan et al., 1998;Kilo et al., 2003;Lau et al., 2011;Li et al., 1995;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) in the past few years.

Basic justification of classical MD
In crystalline solids, ionic conductivity is fundamentally different from electronic conductivity. Electronic conduction in a metal, for example, occurs on a three-dimensional array of ion cores whose excess valence electrons have dissociated to form a continuous "sea of free electrons" partially filling the electronic bands around the Fermi-level. Because the electron has a small mass, its de Broglie wavelength is large and therefore quantum mechanical effects force the electrons into those bands. As ions are much heavier than electrons, their motion is far less governed by quantum mechanics. Below the typical atomic vibrational frequencies (< 100 GHz), ionic motion is best described by thermally activated hopping between (usually) charge-compensating sites (Dyre et al., 2009).
The dynamics of mobile ions in a disordered inorganic ionic conductor (e.g. amorphous YSZ) is clearly a complex multi-particle problem. Unlike a perfect crystal, the potentialenergy landscape experienced by an ion in a disordered solid is irregular and contains a distribution of depths and barrier heights. The ions' interaction with the dynamic atomic lattice network is fundamental: first, because the lattice supplies a persistent disordered potential energy landscape for the mobile anions and second, because the local fluctuations of the lattice atoms promote anionic jumps. Additional multi-particle behavior stems from the interaction among the mobile ions and vacancies. All these distinct coupled interactions contribute to the complete theoretical description of ionic conductivity in amorphous solids (Dyre et al., 2009;Lammert et al., 2010). A complete analytical microscopic theory is not available by now and, due to the complexity of the problem, would be extremely difficult to formulate without some basic understanding at the atomistic level. The direct approach to ionic conduction in the YSZ electrolyte of the SOFC is via molecular dynamics (MD) simulations. From the time-evolved atomic trajectories detailed microscopic information about the underlying mechanisms is available. Understanding them would make theoretical prediction possible.

Interatomic potentials
It is possible to model a few thousands to millions of particles in classical MD using phenomenological interatomic and intermolecular potentials. They are obtained by using the phenomenological approach of selecting a parameterized mathematical form for the interaction between atoms, and fitting its unknown parameters to various experimental or higher-level theoretical (e.g. ab initio quantum mechanics simulation) properties. In general, the flexibility, accuracy, transferability, and computational efficiency of the interatomic potentials each have to be carefully considered (Frenkel & Smit, 1996).
For the YSZ electrolyte in a SOFC, the simplest, relevant interatomic potentials are those commonly used to describe rigid ionic compounds (Lewis et al., 1985). Under the rigid ion model, the potential energy is a simple function of the distance between the ions. It consists www.intechopen.com of a Coulomb term to describe the long-range electrostatic interactions between the ions of YSZ (i.e. Zr 4+ , Y 3+ , and O 2− ), and a Born-Meyer-Buckingham (BMB) potential to describe the short-ranged interactions between the ions. The potential energy between ions i and j separated by a distance r i j with ionic charges Q i and Q j is then given by: The exponential term of the short-ranged BMB potential takes account of Pauli repulsion, whereas the r -6 term takes account of the attractive dispersion or van der Waals interaction. To calculate the long-range Coulomb interactions between the ions in the 3D periodically repeated simulation cell, Ewald summations must be used. To maximize the accuracy of these interatomic potentials and have more reliable energy landscape features, these parameters are usually fitted empirically to more accurate ab initio quantum mechanical calculations, i.e. density functional theory (DFT) or quantum chemistry methods, and relevant experimental findings on some known physical properties (Gale et al., 2003).
Of the several sets of interatomic potential parameters available in the literature (Bush et al., 1994;Dwivedi et al., 1990;Lewis et al., 1985;Minervini et al., 2000;Schelling et al., 2001;Zacate et al., 2000), the parameters proposed by Lewis et al (Lewis et al., 1985) and Minervini et al Zacate et al., 2000) were chosen for the initial guess potential functions in the fitting of interatomic potentials for ZrO 2 and Y 2 O 3 crystals. To better capture the correct dielectric constants and lattice vibrational frequencies as described in the system (Gale et al., 2003;Lindan et al., 1993), the effects of environment-dependent electronic polarizability of the ions  can be included in the fiting of the interationic potential through the shell model (Gale et al., 2003;Lau, 2011;Lindan et al., 1993). In the core-shell model, the ionic core and shell are coupled through harmonic spring force constants, k (Gale et al., 2003;Lau et al., 2011) which take into account interactions of different types: between ions, between ions and outer shell electrons, and between outer electrons. Such an approach allows one to take into account the electronic polarizability of ions that is caused by the forces acting between the ion cores and shells. However to probe the basic structural properties and ionic motion at high temperatures at a reasonable cost of computation, the electronic polarizability within the core-shell model mentioned can be ignored (i.e. all the core-shell spring constants k Zr , k Y and k O mentioned are set to zero). This is appropriate in MD simulations in which no electric field is applied, as in the case of superionic conduction at high temperature in SOFC. The explicit effect of the core-shell model was found to be small (Lindan et al., 1993), and therefore the approximate description of rigid ions (Eq. 1) without core-shell interaction in molecular dynamics should be sufficient.
Stabilized YSZ in its cubic fluorite structure has cations occupying an fcc lattice and oxygen anions occupying its tetrahedral interstices. A finely tuned ZrO 2 and Y 2 O 3 interatomic potential that describes the entire Zr 1−x Y x O 2−x/2 system, with x/2 being the Y 2 O 3 dopant concentration of YSZ has to be established. According to the reported literature (Bush et al., 1994;Devanathan et al., 2006;Dwivedi et al., 1990;Fisher et al., 1998Fisher et al., , 1999Khan 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), the typical "semi-empirically" fitted properties of ZrO 2 and Y 2 O 3 crystals that are chosen for the fitting dataset can be the www.intechopen.com lattice parameters, lattice elastic properties, dielectric constants, defect formation energies (e.g. vacancies and interstitials), and phonon frequencies of cubic (c-ZrO 2 , space group Fm3m), tetragonal (t-ZrO 2 , space group P42/nmc) (Ackermann et al., 1975;Aldebert et al., 1985;Boysen et al., 1991;Dash et al., 2004;Howard et al., 1988;Smith et al., 1965;Zhao et al., 2002;), monoclinic (m-ZrO2, space group P2 1 /c ), and yttria (Baller et al., 2000; Compared to determining the geometry and energetics of the competing phases of ZrO 2 polymorphs, searching for ground state atomic arrangements across a composition range of ZrO 2 -Y 2 O 3 is definitely more complicated. The YSZ solid inherits the complexity of the competing phases of ZrO 2 and adds the possiblility that any Zr atom can be replaced by a Y atom and half of a vacancy. The structures and lattices are not merely determined by different composition at different ambient condition, but also dictated by intrinsic long-and short-range order in the system (Bogicevic et al., 2001). Fortunately, under normal conditions, the cubic fluorite scaffold is found to be stable for yttria (Y 2 O 3 ) content in the range of 8 -40 mol% (Bogicevic et al., 2001;Ostanin et al., 2002Ostanin et al., , 2003Predith et al., 2008). For the SOFC application, the ionic conductivity (and oxygen self-diffusivity) of YSZ does not increase monotonically with increasing vacancy concentration; rather, it exhibits a maximum between 8 and 15 mol% Y 2 O 3 . These unique characteristics must also determine the basic parameter "tune-up" of the interatomic potential for the MD simulation of YSZ. To ensure neutrality of the simulation cell and that it obeys chemical stoichiometry, the cations and anions of YSZ are typically chosen to be Zr 4+ , Y 3+ and O 2− ions.
To simplify the simulation further, the interaction of Zr 4+ -Y 3+ in YSZ is assumed to be governed by the Coulomb interaction of the two ionic charges. This assumption is based on the fact that at the low Y 2 O 3 dopant concentrations of this study, a very strong first-neighbor interaction between Zr 4+ and Y 3+ in the lattice is very unlikely compared to the firstneighbor interactions between oxygen anions and vacancies (Pietrucci et al., 2008;Schelling et al., 2001;Stapper et al., 1999). This approach has been widely adopted in previous studies literature (Bush et al., 1994;Devanathan et al., 2006;Dwivedi et al., 1990;Fisher et al., 1998Fisher et al., , 1999Khan 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) In this dilute Y 2 O 3 concentration limit of YSZ, the local environments of oxygen atoms in cubic YSZ crystals are assumed to be more like that in ZrO 2 than in Y 2 O 3 , therefore the O-O potential adopted throughout the simulation can be approximated to be identical to the O-O potential of ZrO 2 (Lau et al., 2011). As one of the "semi-empirical fitted" BMB potentials that can be found in the literatures, the relevant interatomic potentials (Lau et al., 2011) is shown in Table 1. Besides describing the cubic ZrO 2 phase (i.e. c-ZrO 2 ) well, these potentials can also describe the tetragonal phase of ZrO 2 (i.e. t-ZrO 2 ) and the Y 2 O 3 crystalline phase well, as pointed out in a recent paper (Lau et al., 2011).

Atomistic models for YSZ solids
For SOFC, the challenges that delay the full commercialization include materials degradation, materials selection, materials function, and coupled interactions with other cell components. For the ion-conducting electrolyte like YSZ, the main problems can be attributed to the need for chemical, mechanical, and thermodynamic stability over a wide range of operating conditions. In addition, the optimum ionic conductivity of the electrolyte varies with different synthesis routes and sintering conditions due to the resultant diverse local morphologies, grain boundaries, and microstructures. To facilitate the engineering design, determine materials functionality, and extend experimental findings, classical MD simulation has a unique role within the multi-scale atomistic modeling to provide the relevant predictions that complement standard continuum modeling. One example is its ability to generate a large-scale atomistic structure exhibiting different morphologies and microstructures precisely with reasonable computation cost.
Interaction simulation box, an O 2− ion was removed from randomly selected anion sites for every two Y3+ dopant ions in the system. An amorphous structure, which is different from crystalline YSZ, was generated via standard structural relaxation, the "amorphization and recrystallization (A&R) strategy" proposed by Sayle et al. (Sayle et al., 1999(Sayle et al., , 2002(Sayle et al., , 2005. By adopting the standard MD techniques that are implemented in the standard classical MD codes e.g. LAMMPS and DL_POLY, a thermal well-equilibrated amorphous solid of YSZ as shown in Fig. 3b can be obtained (Lau et al., 2011).

Static and dynamic properties of YSZ bulk solids
In this section, we present a brief discussion of the principal static and dynamic properties of the YSZ bulk solids that primarily focuses on the unique characteristics of YSZ in crystalline ( Fig. 3a) and disordered bulk lattices (Fig. 3b). In particular, the accuracy of the simple atomic potential (Table 1) will be validated, and the implications of these MD predictions will be discussed. For abbreviation, all YSZ crystals will be termed n-YSZc, where n gives the mol% doping (e.g. 8-YSZc is a crystal with 8.0 mol% Y 2 O 3 dopant). The corresponding amorphous YSZ structures will be termed n-YSZa. Both the crystal and amorphous structures referred to together will be termed n-YSZ systems

Lattices
In contrast to the disordered lattice of the amorphous structure, the YSZ crystals that are typically found in the cubic fluorite structure can be fully characterized by a single lattice constant, a. From Fig. 4, the computed lattice parameters of cubic YSZ as a function of doping and temperature are generally consistent with experimental observation (Hayashi et al., 2005;Pascual et al., 1983). This further validates the interatomic potentials we employed in Table 1. Because of the larger size of the dopant Y 3+ cation, which yields larger Y-O bond distances in the lattice, the cell volume of YSZ crystals generally increase with increased doping with Y 2 O 3 . However at low temperature, this trend is not uniform. From Fig. 4, the lattice constant a of 3-YSZc is ~ 5.125 Å, is larger than 8-YSZc (a ~ 5.121 Å), but less than 12-YSZc (a ~ 5.128 Å). These discrepancies might be a signature of mixing the low-temperature tetragonal ZrO 2 ground state and high-temperature cubic structures for the 3-YSZc (or any other low mol% of Y 2 O 3 YSZ crystals in the dilute regime) at low temperature or of the limitation of the current fitted atomic potential used (Table 1), which is unable to capture accurately the phase stability of all YSZ polymorphism at low mol% of Y 2 O 3 dopant in the low-temperature regime.
Within the linear response regime, the thermal expansion of a cubic crystal lattice can be described as a(T ) = a(300)[1 + α(T − 300)], where a(300) is the lattice constant for YSZ solid at 300 K and α is the linear thermal expansion coefficient of the system. From a linear fit of MD results, the coefficients of linear thermal expansion are found to lie in the range ~ 6.0-7.0x10 -6 K -1 (Lau, 2011), close to the reported experimental (Hayashi et al., 2005;Ingel et al., 1986;Pascual et al., 1983) values of ~9.6-10.8x10 -6 K -1 and other theoretical values (Devanathan et al., 2006) (Hayashi et al., 2005;Pascual et al., 1983) (in triangles) at 300 K.
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% Y 2 O 3 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 c 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 in disordered solids.

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 Y 3+ within a cubic Zr (1−x) Y x O (2−x/2) phase, the distribution of the Y 3+ 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 Zr 4+ (cation), O 2− (anion) crystalline planes in both directions shown in Fig. 5 for the crystal lattice. The dopant (Y 3+ ) ions have a similar distribution to that of the Zr 4+ ions that they replace in the crystal of YSZ (Fig. 5). For both systems (i.e. the crystal and amorphous structure), the ratio of Zr 4+ /O 2−~ 1/2 can be found in all crystallographic planes, analogous to the ZrO 2 stoichiometry which acts as a host lattice. In addition, these similar features are also found in the YSZ systems with other Y 2 O 3 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 www.intechopen.com 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 Zr 4+ , Y 3+ and O 2-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.  (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). 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 ZrO 2 , and N Zr = 8 for the tetragonal ZrO 2 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 www.intechopen.com 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 Y 2 O 3 crystal.
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 ZrO 2 and Y 2 O 3 . 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 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 Y 3+ 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 Y 2 O 3 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. www.intechopen.com

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.

i. Temperature-dependent self diffusion and ionic conductivity of O 2-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. <r 2 (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: 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 D i for species i from the MSD using the Einstein relation, <r 2 i (t) > = 6D i t+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 = N i q 2 i D i /H r k B T, where N i is the charge density of charge carriers, q i , 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 Zr 4+ , Y 3+ ) are less mobile than the anions (e.g. O 2-, 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 D i = D 0 exp(-Δ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 www.intechopen.com 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 cm 2 s -1 and 0.59±0.05 eV, respectively, in close agreement with previously reported 6.08×10 -5 cm 2 s -1 , and 0.60 eV theoretical results (Devanathan et al., 2006), but slightly smaller than ~ 1.0 eV for 10.0 mol% Y 2 O 3 YSZ reported from experiment (Kilo et al., 2003). Similarly for the other YSZ crystals in the range of 3.0-12.0 mol% Y 2 O 3 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 Y 2 O 3 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% Y 2 O 3 . 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 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 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 cm 2 s -1 is obtained, whereas ~5.83±0.47×10 -5 cm 2 s -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 O 2− ions. Specifically, the Zr 4+ and Y 3+ ions are found to be lower in coordination number to the O 2− 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% Y 2 O 3 ) solids is ~ 0.45-0.61 eV, comparable to YSZ crystals with the same Y 2 O 3 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. Y 2 O 3 -dopant concentration dependent self diffusion and ionic conductivity of O 2-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, D i  x (with x defining Y 2 O 3 dopant concentration) in the Zr 1−x Y x O (2-x/2) system. The MSD of oxygen ions cannot increase proportional to molar concentration of the mol% Y 2 O 3 because Y 2 O 3 is an insulator. An optimal concentration of Y 2 O 3 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% Y 2 O 3 (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. Y 3+ ) 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(Morinaga et al., , 1980Steele et al., 1974) indicate that the www.intechopen.com 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 Zr 4+ ions than Y 3+ ions. Hence, the dopant ion has a higher oxygen coordination number than the host Zr 4+ 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 Y 3+ -ion is negatively charged, relative to the host lattice, Zr 4+ -ion). Therefore to capture the Y 2 O 3 -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 Zr 4+ ions are smaller than Y 3+ ions, Zr 4+ 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% Y 2 O 3 (or x ~0.16 in Zr (1−x) Y x O (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 Y 2 O 3 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% Y 2 O 3 systems (i.e. n-YSZa) studied; it was found that σ max occurs at ~ 8.0 mol% Y 2 O 3 similar to the maximum found for the YSZ crystals. To account for this intriguing result, one simple, yet direct explanation can be the local Y 3+ 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 (A ij and ij in Tables 1) in the BMB potential, larger bond length (i.e. R Y−O from RDF at Fig. 6) and larger Y 3+ ion relative to Zr 4+ ion in the lattice, steeper potential wells among Y 3+ -O 2− relative to Zr 4+ -O 2− are expected. Thus, this will hinder the mobility of O 2− (or its vacancy), if the local Y 3+ aggregation is large enough. At a given temperature, the local Y 3+ aggregation can be defined qualitatively by comparing to the degree of Y 3+ -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 N c is the number of these Y atoms that have four or more neighboring Y atoms, then η = N c /N×100.0%. As shown in Table 2, η increased as the mol% of Y 2 O 3 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 predicted ionic conductivity, σ (in Scm −1 ), of 3.0-12.0 mol% Y 2 O 3 YSZ crystal over 800-1200 K, with the conductivity maxima as a function of Y 2 O 3 concentration based on the simple interatomic potential (Lau et al., 2011) of Table 1.

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. www.intechopen.com 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. ΔE act 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 <r 2 (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. )] to be ~ 3.7 Å as found by RDF, we define η = N c /N×100.0%, where N c is the number of clusters with ≥ 4 Y 3+ -ions and N is the total number of Y 3+ -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 O 2− ) 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., 1998Fisher et al., , 1999Khan 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 (Zr 4+ and Y 3+ ) 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 YSZ amorphous phase, the heavy cations are not simply bound to a rigid crystal lattice as in YSZ crystals (Fig. 9). 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 Zr 4+ and Y 3+ ions are rather substantial (i.e. <r 2 >~ 0.5 Å 2 at 800 K within 2.5 ns), are almost as diffusive (Fig. 9) as the mobile O 2− anions. For YSZ www.intechopen.com amorphous structure from MD trajectories, the three constituent ions (Zr 4+ , Y 3+ , O 2− ) 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: G s (r,t) = <Σ j δ(r+r j (0)-r j (t))>/N j (3) and the distinct part: G d (r,t) = <Σ i≠j δ(r+r i (0)-r j (t))>/N i (4) which are determined by the space-and time-dependent density correlations of the constituent ions in the system. As defined, G s (r,t) compares the position of a particle to its position within a time interval. In thermal equilibrium, G s (r,t) depends only on the time difference, and provides information about the hopping and diffusion processes of the separated ions. Whereas G d (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, G s (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 G s (r,t) generally are broader than the stationary ion. At low temperature, the G s (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 (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 G s (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 G s (r,t), a comparison of G s (r,t) with the diffusion transport character in Gaussian form is particularly useful. By applying the Laplace transform to Fick's second law (D 2 C =∂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 C(x,y,z,t) = exp(-(x 2 +y 2 +z 2 )/4Dt)/(4 Dt) 3/2 = exp(-r 2 /4Dt)/(4 Dt) 3/2 (5) where r 2 = x 2 +y 2 +z 2 . Assuming that the time response of the van Hove self-correlation function G s (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 G s (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 /4Dt)/(4 Dt) 3/2 where <r 2 > ~ 6Dt 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 <r 2 >. 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 2 C(r, t) in Fig. 10, the decay of the profile and variation found in 4 r 2 G s (r, t) of the YSZ system is prominent as time develops. From its initial localized 'impulse-like' distribution (e.g. 4 r 2 G s (r, t) at 10 ps in Fig. 10a) analogous to a localized -function (Frenkel & Smit, 1996;Rapaport, 2004), the variation of G s (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 G s (r, t) up to 2500 ps, yet are not perfectly identical. In both cases, the non-Gaussian features of the mobile O 2− in G s (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 G s (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 www.intechopen.com Fig. 10. The self-part of the van Hove correlation function 4 r 2 G s (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 (Zr 4+ , Y 3+ ) and anions (O 2− ) however is more obvious, with the long dispersive tails of G s (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 www.intechopen.com 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 O 2− in the amorphous solid will be less compared to the crystal (Fig. 9).

Conclusion
During recent decades, it has become feasible to simulate a complicated system on a computer due to rapid progress in parallel computing. Within the scale of atomistic and molecular simulation, the application of the classical molecular dynamics (MD) simulation method covers a vast variety of systems undergoing current scientific development. The method of MD solves the classical equations of motion for an ensemble of atoms. It results in time-dependent trajectories for all atoms in a system. From these atomistic trajectories, MD can provide detailed in situ atomistic information that is difficult to obtain experimentally. As one of the robust and well-developed simulation techniques, MD simulation is an ideal scientific tool to complement experimental observations to properly characterize a complicated system. Those that involve a vast time and spatial dimension, and heterogeneous materials interfaces can be modeled using a multi-scale framework. One of the best such examples that are becoming appropriate for MD is the study of the inner workings of a solid oxide fuel cell (SOFC), which is important as an electrochemical energy conversion and clean energy storage.
SOFC is a new alternative clean energy device that converts the energy of combustion and electrochemical interactions into electricity, which utilizes the superionic conductivity (> 10 -1 Scm -1 ) of special materials at high temperature. Despite the advantages over conventional power generation technologies, there remain a number of challenges that delay the full commercialization of the SOFC and one of the challenges is to understand the basis of ionic transport in its solid electrolyte (e.g. YSZ, the Zr 1−x Y x O 2−x/2 system, with x/2 being the Y 2 O 3 dopant concentration). For YSZ, the optimum ionic conductivity can vary with different synthesis routes and sintering conditions due to the resultant diverse local morphologies, grain boundaries, and microstructures.
Within the rigid ion model approximation, we have shown a systematic study of the static and dynamic properties of YSZ crystals and amorphous solids within the typical dilute Y 2 O 3 concentration limit (i.e. 3.0 -12.0 mol% Y 2 O 3 ) in the temperature range 300 -1400 K based on a simple semi-empirical Born-Meyer-Buckingham (BMB) interatomic potential through the standard techniques in classical MD. The results suggest that the vacancy assisted ion conductivity of YSZ as a function of mol% of Y 2 O 3 at a given temperature seems to be a universal feature of YSZ electrolytes, regardless of varying local structures.
Whether it is an amorphous or a crystal, the oxygen ionic conductivity shows a maximum at 8.0mol%Y 2 O 3 , close to the lower limit of the cubic YSZ phase stability that is confirmed by experiments. For YSZ amorphous solids, their lower absolute ionic conductivity relative to YSZ crystals is consistent with the trends observed in YSZ crystalline and stabilized amorphous thin films reported from experiments. For the YSZ www.intechopen.com 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 desirable.

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