**2.1 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 of a Coulomb term to describe the long-range electrostatic interactions between the ions of YSZ (i.e. Zr4+, Y3+, and O2−), 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 *ri j* with ionic charges *Qi* and *Qj* is then given by:

$$\mathbf{E}\_{i\bar{\imath}} = \mathbf{Q}\_i \mathbf{Q}\_{\bar{\jmath}} / \mathbf{r}\_{\bar{\imath}\bar{\jmath}} + \mathbf{A}\_{\bar{\imath}\bar{\imath}} \exp(-\mathbf{r}\_{\bar{\imath}\bar{\jmath}} / \rho\_{\bar{\imath}\bar{\jmath}}) - (\mathbf{C}\_{i\bar{\jmath}} / \mathbf{r}^{\phi\_{\bar{\imath}}}) \tag{1}$$

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

lattice parameters, lattice elastic properties, dielectric constants, defect formation energies (e.g. vacancies and interstitials), and phonon frequencies of cubic (*c*-ZrO2, space group *Fm*3*m*), tetragonal (*t*-ZrO2, space group *P*42*/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 P21/c ), and yttria (Baller et al., 2000; Lau et al.,

Compared to determining the geometry and energetics of the competing phases of ZrO2 polymorphs, searching for ground state atomic arrangements across a composition range of ZrO2–Y2O3 is definitely more complicated. The YSZ solid inherits the complexity of the competing phases of ZrO2 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 (Y2O3) content in the range of 8 – 40 mol% (Bogicevic et al., 2001; Ostanin et al., 2002, 2003; Predith 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% Y2O3. 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

To simplify the simulation further, the interaction of Zr4+–Y3+ 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 Y2O3 dopant concentrations of this study, a very strong first-neighbor interaction between Zr4+ and Y3+ 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., 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) In this dilute Y2O3 concentration limit of YSZ, the local environments of oxygen atoms in cubic YSZ crystals are assumed to be more like that in ZrO2 than in Y2O3, therefore the O–O potential adopted throughout the simulation can be approximated to be identical to the O–O potential of ZrO2 (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 ZrO2 phase (i.e. *c-*ZrO2) well, these potentials can also describe the tetragonal phase of ZrO2 (i.e. *t-*ZrO2) and the Y2O3 crystalline phase well, as

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

and anions of YSZ are typically chosen to be Zr4+, Y3+ and O2− ions.

pointed out in a recent paper (Lau et al., 2011).

**2.2 Atomistic models for YSZ solids** 

2009) (Y2O3, space group *Ia*3).

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* (Minervini et al., 2000; Zacate et al., 2000) were chosen for the initial guess potential functions in the fitting of interatomic potentials for ZrO2 and Y2O3 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 (Lau et al., 2009) 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 ZrO2 and Y2O3 interatomic potential that describes the entire Zr1−*<sup>x</sup>*Y*x*O2−*x/*2 system, with *x/*2 being the Y2O3 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., 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), the typical "semi-empirically" fitted properties of ZrO2 and Y2O3 crystals that are chosen for the fitting dataset can be the

of a Coulomb term to describe the long-range electrostatic interactions between the ions of YSZ (i.e. Zr4+, Y3+, and O2−), and a Born–Meyer–Buckingham (BMB) potential to describe the short-ranged interactions between the ions. The potential energy between ions *i* and *j* 

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* (Minervini et al., 2000; Zacate et al., 2000) were chosen for the initial guess potential functions in the fitting of interatomic potentials for ZrO2 and Y2O3 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 (Lau et al., 2009) 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

Stabilized YSZ in its cubic fluorite structure has cations occupying an fcc lattice and oxygen anions occupying its tetrahedral interstices. A finely tuned ZrO2 and Y2O3 interatomic potential that describes the entire Zr1−*<sup>x</sup>*Y*x*O2−*x/*2 system, with *x/*2 being the Y2O3 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., 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), the typical "semi-empirically" fitted properties of ZrO2 and Y2O3 crystals that are chosen for the fitting dataset can be the

ij) (1)

separated by a distance *ri j* with ionic charges *Qi* and *Qj* is then given by:

Ei j = Qi Qj/rij + Aij exp(−rij/ρij)−(Ci j/r6

sufficient.

lattice parameters, lattice elastic properties, dielectric constants, defect formation energies (e.g. vacancies and interstitials), and phonon frequencies of cubic (*c*-ZrO2, space group *Fm*3*m*), tetragonal (*t*-ZrO2, space group *P*42*/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 P21/c ), and yttria (Baller et al., 2000; Lau et al., 2009) (Y2O3, space group *Ia*3).

Compared to determining the geometry and energetics of the competing phases of ZrO2 polymorphs, searching for ground state atomic arrangements across a composition range of ZrO2–Y2O3 is definitely more complicated. The YSZ solid inherits the complexity of the competing phases of ZrO2 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 (Y2O3) content in the range of 8 – 40 mol% (Bogicevic et al., 2001; Ostanin et al., 2002, 2003; Predith 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% Y2O3. 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 Zr4+, Y3+ and O2− ions.

To simplify the simulation further, the interaction of Zr4+–Y3+ 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 Y2O3 dopant concentrations of this study, a very strong first-neighbor interaction between Zr4+ and Y3+ 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., 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) In this dilute Y2O3 concentration limit of YSZ, the local environments of oxygen atoms in cubic YSZ crystals are assumed to be more like that in ZrO2 than in Y2O3, therefore the O–O potential adopted throughout the simulation can be approximated to be identical to the O–O potential of ZrO2 (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 ZrO2 phase (i.e. *c-*ZrO2) well, these potentials can also describe the tetragonal phase of ZrO2 (i.e. *t-*ZrO2) and the Y2O3 crystalline phase well, as pointed out in a recent paper (Lau et al., 2011).

#### **2.2 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

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

Fig. 3. Various morphologies of YSZ solids (i.e. Zr4+ is blue, Y3+ is yellow and O2- is red ion) that can be modeled using classical MD simulation (a) crystal, (b) amorphous structure, (c) grain boundary, (d) crystalline thin film structure, and (e) the finite nanocrystalline

structure.

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.


Table 1. Short-range interatomic BMB potential for the YSZ solids and their predicted properties compared to *ab initio* DFT planewave calculations within the Local Density Approximation (LDA) (Lau et al., 2009). Here, the Ω-YSZ (i.e. Zr5Y2O13 solid) consists of 17 mol% yttria, whereas the Δ-YSZ (i.e. Zr3Y4O12) consists of 40 mol% yttria. The *ωmin* and *ωmax* refer to the lowest and highest vibration frequency of phonons for the corresponding YSZ system.

Within the rigid-ion approximation validated in Table 1, various types of YSZ solids that arise from different experimental synthesis routes and sintering conditions can be modeled using the appropriate spatial boundary conditions under realistic thermodynamic conditions (Sect. 2.1). Those include the perfect crystal, amorphous structures, grain boundaries, thin-film surfaces, and nanocrystals of Fig. 3. For the YSZ systems shown in that figure, the number of Zr4+, Y3+, and O2− ions depends on the chosen Y2O3 dopant concentration. Initially, the cubic fluorite lattice of ZrO2 was constructed, and yttrium Y3+ ions were generated by random replacements of zirconium Zr4+ ions according to the predefined chemical composition. To keep stoichiometry and the neutrality of the

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

Interaction *Ai j* (eV) *ρi j* (Å) *Ci j* (eV Å6)

System Property This Work DFT Ω-YSZ *a* (Å) 6.26 6.40

Δ-YSZ *a* (Å) 6.23 6.37

Table 1. Short-range interatomic BMB potential for the YSZ solids and their predicted properties compared to *ab initio* DFT planewave calculations within the Local Density Approximation (LDA) (Lau et al., 2009). Here, the Ω-YSZ (i.e. Zr5Y2O13 solid) consists of 17 mol% yttria, whereas the Δ-YSZ (i.e. Zr3Y4O12) consists of 40 mol% yttria. The *ωmin* and *ωmax* refer to the lowest and highest vibration frequency of phonons for the corresponding YSZ

Within the rigid-ion approximation validated in Table 1, various types of YSZ solids that arise from different experimental synthesis routes and sintering conditions can be modeled using the appropriate spatial boundary conditions under realistic thermodynamic conditions (Sect. 2.1). Those include the perfect crystal, amorphous structures, grain boundaries, thin-film surfaces, and nanocrystals of Fig. 3. For the YSZ systems shown in that figure, the number of Zr4+, Y3+, and O2− ions depends on the chosen Y2O3 dopant concentration. Initially, the cubic fluorite lattice of ZrO2 was constructed, and yttrium Y3+ ions were generated by random replacements of zirconium Zr4+ ions according to the predefined chemical composition. To keep stoichiometry and the neutrality of the

1*.*29285x103 3*.*58388x10<sup>−</sup>1 1*.*93646x101

1*.*642724x103 3*.*53197x10<sup>−</sup>1 1*.*04180x102

1*.*30989x104 2*.*19670x10<sup>−</sup>1 4*.*92998x101

*b* (Å) 6.30 6.20 *c* (Å) 6.29 6.30 *α* (deg) 100.71 99.96 *β* (deg) 100.13 100.75 *γ* (deg) 99.69 98.57 *ωmin* (cm-1) 101 114 *ωmax* (cm-1) 696 719

*b* (Å) 6.29 6.28 *c* (Å) 6.37 6.48 *α* (deg) 81.23 81.53 *β* (deg) 99.49 99.59 *γ* (deg) 80.59 79.98 *ωmin* (cm-1) 146 115 *ωmax* (cm-1) 678 707

microstructures precisely with reasonable computation cost.

Zr4+–O2<sup>−</sup>

Y3+–O2<sup>−</sup>

O2<sup>−</sup> –O2<sup>−</sup>

system.

Fig. 3. Various morphologies of YSZ solids (i.e. Zr4+ is blue, Y3+ is yellow and O2- is red ion) that can be modeled using classical MD simulation (a) crystal, (b) amorphous structure, (c) grain boundary, (d) crystalline thin film structure, and (e) the finite nanocrystalline structure.

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

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.

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 shown in Fig. 3b can be obtained (Lau et al., 2011).
