**3. 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% Y2O3 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

### **3.1 Static properties**

#### **3.1.1 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 Y3+ cation, which yields larger Y-O bond distances in the lattice, the cell volume of YSZ crystals generally increase with increased doping with Y2O3. 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 ZrO2 ground state and high-temperature cubic structures for the 3-YSZc (or any other low mol% of Y2O3 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 Y2O3 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) of 6.8–8*.*0x10-6 K-1 for other Y2O3 doping. For the amorphous YSZ, the low symmetry yields the irregularity in local structure that gives distinct lattice properties (i.e. *a*, *b* and *c)*. Under the homogeneous quenching via Hoover barostats and

simulation box, an O2− 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, 2002, 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

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% Y2O3 dopant). The corresponding amorphous YSZ structures will be termed *n*-YSZa. Both the crystal and

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 Y3+ cation, which yields larger Y-O bond distances in the lattice, the cell volume of YSZ crystals generally increase with increased doping with Y2O3. 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 ZrO2 ground state and high-temperature cubic structures for the 3-YSZc (or any other low mol% of Y2O3 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 Y2O3 dopant in the

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) of 6.8–8*.*0x10-6 K-1 for other Y2O3 doping. For the amorphous YSZ, the low symmetry yields the irregularity in local structure that gives distinct lattice properties (i.e. *a*, *b* and *c)*. Under the homogeneous quenching via Hoover barostats and

shown in Fig. 3b can be obtained (Lau et al., 2011).

**3.1 Static properties** 

low-temperature regime.

**3.1.1 Lattices** 

**3. Static and dynamic properties of YSZ bulk solids** 

amorphous structures referred to together will be termed *n*-YSZ systems

Fig. 4. (Top) Lattice constant (in Å) of YSZ crystals with various Y2O3 mol% dopant concentrations as a function of temperature (Lau et al., 2011) , compared with reported experimental values (Hayashi et al., 2005; Pascual et al., 1983) (in triangles) at 300 K. (Bottom) Relative volume expansion (in %) referenced to 300 K of 3-YSZc (green dots), 8- YSZc (yellow dashes), 12-YSZc (red lines) systems and corresponding amorphous solids (lines with points) over the temperature range 300–1000 K.

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

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)

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

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).

two angstroms, along these two crystallographic directions (Lau et al., 2011).

are due to their unique local ion distributions and system densities.

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 in disordered solids.
