Physical description of the monoclinic phase of zirconia based on the bond-order characteristic of the Tersoff potential

Abstract

Zirconia has many important phases with Zr coordination varying from six-fold in the orthorhombic phase to eight-fold in the cubic and tetragonal phases. Development of empirical potentials to describe these zirconia phases is an important but long-standing challenge, and it is a bottleneck for theoretical investigation of large zirconia structures. Here, instead of using the standard core–shell model, we developed a new potential for zirconia by combining the long-range Coulomb interaction and bond-order Tersoff model. The bond-order characteristic of the Tersoff potential enables it to be well suited to describe the zirconia phases with different coordination numbers. In particular, the complex monoclinic phase with two inequivalent oxygen atoms, which is difficult to describe with most existing empirical potentials, is well described by this newly developed potential. This potential provides reasonable predictions of most of the static and dynamic properties of various zirconia phases. Besides its clear physical essence, this potential is at least one order of magnitude faster than core–shell based potentials in molecular dynamics simulation. This is because it does not include an ultralight shell that requires an extremely small time step. We also provide potential scripts for the widely used simulation packages GULP and LAMMPS.

Graphical abstract

Keywords

zirconia / ZrO₂ / empirical potential / molecular dynamics simulation

1 Introduction

Zirconia-based ceramics are important industrial materials with high thermal stability and high strength [1]. Yttria-stabilized zirconia acts as a thermal barrier coating material to protect the substrate components from hot gases in turbines and engines [2]. Zirconia can be used as an oxygen sensor or high quality oxygen ion channels, and artificial diamonds can be produced based on zirconia [3].

Owing to its industrial and military importance, zirconia and zirconia-based materials have attracted intense global research interest. Many properties of zirconia-based materials were first discovered by experiments rather than theory. For instance, application of yttria-stabilized zirconia as a thermal barrier coating material was achieved by experiment [2], while theoretical achievements lag behind experimental achievements. With the increase of computer power, an increasing number of ab initio calculations have been performed to study various properties of zirconia [4–10]. Ab initio calculations have high accuracy, but they are computationally expensive. Empirical potentials are desirable for investigating systems with a very large number of degrees of freedom.

As long as reliable empirical potentials remain the foundation of theoretical research, significant efforts will be devoted to developing empirical potentials for zirconia. Zirconia is an ionic oxide, so its interactions are dominated by Coulomb interactions. Classically, the long-range attractive Coulomb interaction is balanced by the short-range Born–Mayer repulsive interaction [11], which originates from Pauli repulsion from overlap of the electron density [12]. Several parameter sets of the Born–Mayer model for zirconia are available in the literature [13–15].

Combining Coulomb and Born–Mayer interactions can provide a basic description of some of the properties of cubic zirconia, but it cannot describe the important tetragonal and monoclinic phases. The stability of the tetragonal phase is closely related to the instantaneous polarization of the ions. A standard approach to describe a polarizable ion is to divide the ion into a core and a shell [16]. The core–shell model of zirconia has been parameterized in several studies [17–19]. The effect of polarizable charges can also be considered by introducing a phenomenological charge–dipole term [20]. The instability of the cubic phase and the resultant cubic–tetragonal transition can also be predicted by introducing additional Born–Mayer interactions among the oxygen ions [13].

The monoclinic phase is even more complicated than the tetragonal phase because it contains two inequivalent oxygen ions [1]. The existence of two inequivalent oxygen ions originates from charge redistribution between the oxygen ions in zirconia [21]. To correctly simultaneously describe all of the zirconia phases, the Coulomb interaction, Born–Mayer potential, core–shell model, and charge redistribution effect need to be combined [21].

Other empirical potentials to describe the atomic interactions in zirconia are also available, including the tight-binding model [22], reactive force field model [23], and neural network model [24].

Rather than using the core–shell model, we analyzed the interactions by examining the bond-order properties of zirconia. As mentioned by Smirnov et al. [21], dioxides with relatively large cations (e.g., Th, Ce, and U) are stable in fluorite-like lattices with eight-fold cation coordination, while dioxides with relatively small cations (e.g., Pb, Sn, Ti, and W) are stable in the structures with six-fold cation coordination. Differing from these dioxides, the coordination varies among the different zirconia phases. The lowest-energy state is the monoclinic phase with seven-fold cation coordination. The cubic and tetragonal phases have eight-fold cation coordination, while the orthorhombic phase has six-fold cation coordination. The energy order of these zirconia phases is monoclinic < tetragonal < cubic < orthorhombic, which explicitly shows a strong correlation between the energy and the coordination [1]. In other words, the configuration of zirconia depends on the bond order of the Zr atom. From this bond-order dependence, we believe that the atomic interactions in zirconia can be described by bond-order empirical potentials. The Tersoff model is a widely used empirical potential that possesses an explicit bond-order dependence in its functional form [25]. We propose to combine the Coulomb interaction and Tersoff potential to describe the atomic interactions of zirconia.

In this study, we propose to use the Coulomb and Tersoff (CT) potential to describe the different zirconia phases in terms of the bond order of Zr atoms, while the traditional core–shell concept is avoided. As a consequence, the CT potential is at least one order of magnitude faster than standard core–shell-based empirical potentials, while the clear physical essence is maintained. By taking advantage of the bond-order property in the Tersoff potential, the CT potential can predict the correct energy order of the various zirconia phases, with the monoclinic phase being the lowest energy phase. The instability of the cubic zirconia phase is also predicted by the CT potential. We also applied the CT potential to systematically investigate some static and dynamic properties of different zirconia phases.

2 Structure

Many zirconia phases have been experimentally observed or theoretically discussed [1]. The present work focuses on the four most studied zirconia phases shown in Fig. 1: the cubic, tetragonal, monoclinic, and orthorhombic phases. In cubic zirconia, the Zr atoms occupy the face-centered cubic lattice sites, while the O atoms are in tetrahedral positions. Each Zr atom is coordinated to eight O atoms in a symmetric manner. In tetragonal zirconia, the eight O atoms around the Zr atom are divided into two groups and relatively shifted by \mathrm{\Delta }z along the c axis, and the symmetric unit cell is elongated along this principal axis. The monoclinic phase has a more complex configuration with a monoclinic lattice [Fig. 1 (c)]. The Zr atom has seven-fold coordination. There are two inequivalent O atoms, one with three-fold coordination and another with four-fold coordination. In the orthorhombic phase, the Zr atoms have six-fold coordination, while the O atoms have three-fold coordination [Fig. 1 (d)]. The orthorhombic ZrO₂ phase only exists under high pressure [26], so this phase should have higher energy than the other zirconia phases.

3 Potential model

3.1 Born–Mayer potential

The interactions in ZrO₂ are usually described by the long-range attractive Coulomb interaction and the short-range repulsive Born–Mayer interaction [12]:

\begin{array}{}\text{(1)}& V_{ij}=\frac{Q_i{Q}_j}{r_{ij}}+A{\text{e}}^{-{\lambda }_1{r}_{ij}}\mathrm{,}\end{array}

where r_{ij} is the distance between atoms i and j. The first term is the Coulomb interaction between charges Q_i and Q_j. The second term is the Born–Mayer interaction with the two parameters A and {\lambda }_1. It has been shown that there is a constraint relation between these two parameters [21], which can be obtained as follows. The force of the repulsive term is

\begin{array}{}\text{(2)}& F=A{\lambda }_1{\text{e}}^{-{\lambda }_{1}r}\mathrm{.}\end{array}

The equilibrium Zr–O bond length of the cubic ZrO₂ phase is r_0=2.2 Å. In the equilibrium configuration, the force from the repulsive term is

\begin{array}{}\text{(3)}& F_0=A{\lambda }_1{\text{e}}^{-{\lambda }_1{r}_0}\mathrm{.}\end{array}

The repulsive force from Eq. (3) is balanced by the attractive Coulomb force. We require that the structure of the cubic ZrO₂ phase is the same for different parameters (A, {\lambda }_1), so the attractive Coulomb force does not vary. Consequently, the repulsive force in Eq. (3) does not change either owing to the equilibrium condition. As a result, Eq. (3) becomes a constraint relationship for parameters A and {\lambda }_1. This constraint relation is shown in Fig. 2(a).

For each set of parameters (A, {\lambda }_1), the energies of the cubic and orthorhombic ZrO₂ phases are compared in Fig. 2(b). It is reasonable to presume that the orthorhombic ZrO₂ phase has higher energy than the cubic ZrO₂ phase because the orthorhombic ZrO₂ phase only exists under high pressure [26]. A reasonable parameter set is {\lambda }_1=3.05 Å^–1 and A=2023.6003 eV according to the constraint relation.

The combination of the Coulomb attractive interaction and Born–Mayer repulsive interaction captures the fundamental ionic characteristic features of ZrO₂. As a consequence, several phases of ZrO₂ can be obtained by structure relaxation with this potential, including the cubic, orthorhombic, and monoclinic phases. The tetragonal phase can also be obtained by introducing an additional Born–Mayer repulsive interaction between the oxygen ions [13]. However, using only Coulomb and Born–Mayer interactions, it is not possible to obtain the correct order of the energies of the different phases because this potential does not contain the underlying mechanism of the stability for the tetragonal and monoclinic phases. Many studies have revealed that polarization of the oxygen ion is the physical mechanism driving the cubic–tetragonal phase transition [17–19], that is, polarization of the oxygen ion stabilizes the tetragonal phase of ZrO₂. The core–shell model is usually adopted to describe polarizable ions. Some studies have proposed that charge redistribution over these two inequivalent oxygen ions is the key mechanism to stabilize the monoclinic ZrO₂ phase, and a variable charge model has been developed to describe the stability of the monoclinic ZrO₂ phase [21].

3.2 Tersoff potential

Instead of using the core–shell model and variable charge model, we suggest to describe the zirconia phases by the CT potential, which combines the Coulomb interaction and Tersoff potential. The most significant characteristic of the Tersoff potential is its bond-order property, that is, the strength of each bond depends on its chemical environment. In particular, the bond strength explicitly depends on the coordination numbers of the two atoms forming the bond. As a result, the energy of each atom is dependent on its coordination. Considering that the Zr atoms have different coordination numbers in different zirconia phases, the Tersoff potential is suitable for describing the interactions in zirconia.

The CT potential takes the following form:

\begin{array}{}\text{(4)}& V_{ij}=\frac{Q_i{Q}_j}{r_{ij}}+V_{ij}^t\mathrm{,}\end{array}

where the first term is the standard Coulomb interaction. The parameters related to the Coulomb interaction are listed in Table 1. Summation of the long-range electrostatic interactions was performed by the truncation-based summation approach initially proposed by Wolf. in 1999 [27] and further developed by Fennell and Gezelter in 2006 [28]. We chose the damping parameter \alpha =0.3 Å^–1, which is used in the truncation-based summation approach. We chose the cut-off distance of r_c=10.0 Å. These two parameters (α, r_c) have been used in several previous studies [29, 30].

The second term in Eq. (4) is the Tersoff potential. This potential was proposed by Tersoff in 1986 [25], modified in 1988 [31], and then generalized to multi-component systems in 1989 [32]. There are some minor differences in the notations of different versions. The present work uses the following functional form of the Tersoff potential:

\begin{array}{}\text{(5)}& V_{ij}^t=f_C\left(r_{ij}\right)\left[f_R\left(r_{ij}\right)+b_{ij}{f}_A\left(r_{ij}\right)\right]\mathrm{.}\end{array}

The cut-off function is

The repulsive and attractive terms are

\begin{array}{}\text{(7)}& f_R\left(r_{ij}\right)=A{\text{e}}^{-{\lambda }_1{r}_{ij}};\end{array}

\begin{array}{}\text{(8)}& f_A\left(r_{ij}\right)=-B{\text{e}}^{-{\lambda }_2{r}_{ij}}\mathrm{.}\end{array}

It should be noted that the repulsive term in the Tersoff potential is exactly the same as the Born–Mayer potential. Hence, the values of the parameters A and {\lambda }_1 in the Tersoff potential are set to the same as those of the Born–Mayer potential, that is, λ₁=3.05 Å^–1 and A=2023.6003 eV.

Following the Morse potential [33], the {\lambda }_2 parameter in the attractive term can be set to {\lambda }_2={\lambda }_1/2 (i.e., {\lambda }_2=1.525 Å). To determine the energy parameter B, it should be noted that both the Coulomb interaction and the f_A term in the Tersoff potential are attractive interactions. The ionic bond model is most suitable for ionic crystals consisting of atoms from columns I and VII in the periodic table, such as NaCl. This is because the energy of the valence electron in a metallic atom is much higher than that in a chlorine-like atom, so valence electron transfer is nearly complete, while electron coupling is only a small perturbation. Zirconia consists of a transition metal from column II and oxygen from column VI, both of which are closer to the center of the periodic table than Na and Cl in NaCl. Thus, the energy difference between the valence electrons of the Zr and O atoms is smaller, so valence electron transfer is slightly weakened. As a result, the covalent component will increase, but the ionicity still dominates, which will be reflected by reduction of the ionic charges to the effective charges [34]. Indeed, the charges of the Zr cations and O anions are slightly reduced from their normal values to +3.8 and –1.9, which indicates that a small portion of the ionic interaction is replaced by a covalent interaction. The variation of the charge is chosen to be the same as that in the charge redistribution model, in which approximately 5% of the charges are transferred between oxygen ions [21]. The B parameter was fitted to the structure of cubic zirconia. The fitted parameter (B=17.3376 eV) is around two orders of magnitude smaller than the A parameter, which further confirms that the interaction is mostly ionic.

The characteristic feature of the Tersoff potential is the bond-order term

\begin{array}{}\text{(9)}& b_{ij}={\left(1+{\beta }^n{\zeta }_{ij}^n\right)}^{-\frac{1}{2n}}\mathrm{.}\end{array}

The effective coordination {\zeta }_{ij} includes the local environment effect through the following expression:

\begin{array}{}\text{(10)}& {\zeta }_{ij}={\displaystyle \sum _{k\ne i\mathrm{,}j}}{f}_c(r_{ik})g({\theta }_{ijk}){\text{e}}^{{\lambda }_3^m{(r_{ij}-r_{ik})}^m}\mathrm{,}\end{array}

where the summation {\displaystyle {\sum }_k} is over other bonds i-k around atom i. The coordination of atom i is closely related to {\zeta }_{ij}, which is regarded as the effective coordination of atom i. The three-body term is

\begin{array}{}\text{(11)}& g({\theta }_{ijk})=1+\frac{c^2}{d^2}-\frac{c^2}{d^2+{\left(h-\cos {\theta }_{ijk}\right)}^2}\mathrm{.}\end{array}

Before presenting the parameters of the Tersoff potential, we will discuss the suitability of the Tersoff potential for describing the monoclinic phase by taking advantage of its bond-order characteristic. Let us consider a simple situation with g({\theta }_{ijk})=1 and {\lambda }_3=0.0 in Eq. (10), so the actual coordination of atom i is {\zeta }_{ij}+1. Because of this close relationship between {\zeta }_{ij} and the coordination number, {\zeta }_{ij} is usually called the effective coordination number. The energy minimum position r_m is determined by

\begin{array}{}\text{(12)}& \frac{\partial V_{ij}^t}{\partial r}{|}_{r=r_m}=\mathrm{0,}\end{array}

which gives

\begin{array}{}\text{(13)}& r_m=\frac{1}{{\lambda }_1-{\lambda }_2}\ln \frac{A{\lambda }_1}{bB{\lambda }_2}\mathrm{.}\end{array}

The corresponding energy minimum is

\begin{array}{}\text{(14)}& V_{ij}^m=A{\text{e}}^{-{\lambda }_1{r}_m}-bB{\text{e}}^{-{\lambda }_2{r}_m}\mathrm{.}\end{array}

Assuming that all of the bonds around atom i have the same energy, the total energy of atom i can be determined by

\begin{array}{}\text{(14)}& \begin{array}{cccc}{V}_i(\zeta )\hfill & =V_m\times \left(\zeta +1\right)\hfill & \hfill & \left(\text{15}\right)\hfill \\ \hfill & =\left(A{\text{e}}^{-{\lambda }_1{r}_m}-bB{\text{e}}^{-{\lambda }_2{r}_m}\right)\times \left(\zeta +1\right)\mathrm{.}\hfill & \hfill & \left(\text{16}\right)\hfill \end{array}\end{array}

The subscripts i and j have been omitted. Note that b is also a function of ζ, as can be observed in Eq. (9). For a given set of parameters (n, β), the energy V_i(\zeta ) is an explicit function of the effective coordination number ζ. The function V_i(\zeta ) is shown in Fig. 3, which shows that the minimum of V_i(\zeta ) with respective to ζ can be tuned by varying the parameters n and β. The energy V_i(\zeta ) is minimum at \zeta =6 (i.e., actual coordination of 7) for n=5 and \beta =0.24. Thus, the bond-order characteristic of the Tersoff potential is able to ensure that the monoclinic phase (with seven-fold coordination of Zr) is the lowest energy configuration among all of the zirconia phases. We have thus shown that the stability of monoclinic zirconia among the different phases (with different Zr coordination) can be well described by the Tersoff potential using its bond-order property. The parameters n and β for the Zr atom are primarily determined according to Fig. 3. It should be noted that in zirconia phases with inequivalent bond lengths, the bonds around an atom will have different energies, so the assumption in Eq. (16) does not hold in this situation. Consequently, the final n and β parameters slightly deviate from the ideal values in Fig. 3.

We emphasize that Tersoff potential contains both repulsive and attractive terms. The repulsive term in the Tersoff potential is exactly the same as the usual Born–Mayer potential, both of which describe Pauli repulsion owing to overlap of the electron density. The novelty of the CT potential lies in the bond-order dependent attractive term in the Tersoff potential. The bond-order attractive term is only a small fraction of the whole CT potential, but this is the key component that provides good descriptions of both the tetragonal and monoclinic phases.

Here, we summarize some of the major steps in the parameterization process of the CT potential. Normally, the parameters of empirical potentials are numerically fitted to some quantities by minimizing the square root error. In contrast to this numerical fitting procedure, the present CT potential was parameterized in a physical manner, similar to that used in developing the variable charge model [21]. More specifically, the effective charges of the cation (+3.8) and anion (–1.9) were determined by considering the amount of charge redistribution in the variable charge model [21]. The damping factor and the cut-off in the Coulomb summation are typical numerical parameters that have been successfully used in previous studies. Parameters A and {\lambda }_1 were determined by the Born–Mayer model because the repulsive term in the Tersoff potential is the same as the Born–Mayer model. The parameter {\lambda }_2 was determined following the Morse model [33], which is a typical treatment of the Tersoff potential [25]. Parameters n and β are the two most important parameters of the Tersoff potential because these two parameters determine the bond-order property of the Tersoff potential. We fitted n and β in a physical manner, as shown in Fig. 3, where n and β were fitted to ensure that the monoclinic phase was the lowest energy zirconia structure. The final three parameters c, d, and h were fitted to the correct energy order of the different phases of zirconia.

The parameters of the Tersoff potential of ZrO₂ are listed in Table 2. The Tersoff potential files for GULP [35] and LAMMPS [36] are available from the personal website of the corresponding author (jiangjinwu.org).

4 Application of the CT potential

We applied the CT potential to predict some of the typical properties of the different zirconia phases, and then compared the results with available experimental or ab initio results.

4.1 Monoclinic phase

The CT potential includes both the long-range Coulomb interaction and bond-order Tersoff potential. The Tersoff potential has been implemented in most lattice dynamics and molecule dynamics (MD) simulation packages, such as GULP [35] and LAMMPS [36]. In the present work, the GULP [35] package was used to calculate the static properties of the different zirconia phases, including the structural properties, energy barrier, phonon dispersion, Young’s modulus, and Poisson’s ratio.

The structural properties of monoclinic zirconia calculated with the present CT potential are in good agreement with experimental or ab initio results (Table 3), which indicates that the CT potential can successfully describe the monoclinic phase. Furthermore, the monoclinic phase has the lowest energy among all of the zirconia phases (Table 6). The energy barrier between the tetragonal and monoclinic phases \mathrm{\Delta }{E}_{tm}=80.27 meV/molecule is in good agreement with experimental [3] and ab initio results [6, 37]. The volume of the monoclinic zirconia is clearly larger than those of the cubic and tetragonal phases, which agrees well with previous studies. These results confirm that the CT potential correctly predicts that the monoclinic zirconia phase is the most stable phase among the different zirconia phases.

4.2 Tetragonal phase

While the success of the CT potential in describing the monoclinic phase is expected considering the bond-order feature of the Tersoff potential, the Tersoff potential does not have any specific features to guarantee transition of the cubic phase to the more stable tetragonal phase. However, the present CT potential also predicts this transition. To illustrate this fact, the phonon dispersion of cubic zirconia is shown in Fig. 4. A primitive unit cell containing one Zr atom and two O atoms was used for the calculation. The phonon dispersion was calculated with GULP [35]. These phonon branches show similar behavior to the results from ab initio calculations [5, 7–9]. The most significant feature in the phonon dispersion is the softening of one optical branch, which eventually becomes imaginary for wave vectors around the X point in the Brillouin zone. Note that the convention is to show an imaginary value as a negative value for the frequency in phonon dispersion, like in Fig. 4.

The vibrational morphology corresponding to the imaginary mode at the X point (i.e., the {\text{X}}_2^{-} mode) is shown at the bottom of Fig. 4. This figure was plotted with the XCRYSDEN package [38]. The Zr atom does not vibrate, while the eight surrounding oxygen atoms are divided into two groups. These two groups of oxygen atoms vibrate in opposite directions. The cubic phase will be distorted into the tetragonal phase by deformation following the vibrational morphology of the {\text{X}}_2^{-} mode. The {\text{X}}_2^{-} mode is the origin of the transition from cubic zirconia to tetragonal zirconia.

The imaginary mode at the X point can be attributed to instantaneous polarization of the Zr cation or O anions, which can be described by the core–shell model. Adding additional Born–Mayer interactions among the oxygen ions can also provide correct prediction of this imaginary mode. Here, we have provided a third solution for the {\text{X}}_2^{-} imaginary mode, that is, by the bond-order Tersoff potential.

To further explore the relationship between the cubic and tetragonal phases, we investigated evolution of the structure from the cubic phase to the tetragonal phase. We introduce the parameter η for evolution of the structure according to the expression {\overrightarrow{R}}_{\eta }=\frac{1-\eta }{2}{\overrightarrow{R}}_c+\frac{1+\eta }{2}{\overrightarrow{R}}_t, where {\overrightarrow{R}}_c and {\overrightarrow{R}}_t represent the structures of the cubic and tetragonal phases, respectively. The structure with \eta =-1 is the cubic phase, while the structure with \eta =+1 corresponds to the tetragonal phase. For an arbitrary η, the structure {\overrightarrow{R}}_{\eta } is constructed based on the cubic and tetragonal configurations. Evolution of the energy of the structure by varying η is shown in Fig. 5. The cubic phase is located at a local maximum energy position, while the tetragonal phase is at a local minimum energy position. As a result, the cubic phase is unstable and will transform to the tetragonal phase.

4.3 Static and mechanical properties

The structural and energy properties of the zirconia phases are listed in Table 6. The magnitude of the frequency of the {\text{X}}_2^{-} mode is comparable to previous calculations, which again verifies the instability of the cubic phase. The quantity \mathrm{\Delta }z shown in Fig. 1 quantifies the magnitude of cubic–tetragonal distortion. \mathrm{\Delta }z from the present work is smaller than those from previous studies [3, 6, 21, 37], which indicates that cubic–tetragonal distortion is underestimated by the present CT potential. For the same reason, the energy barrier \mathrm{\Delta }{E}_{ct} between the cubic and tetragonal phases is underestimated by the present CT potential. These results show that the present CT potential can only provide a qualitative description of distortion of the cubic phase to the tetragonal phase. The actual magnitude of the distortion is underestimated with the CT potential.

The Young’s modulus and Poisson’s ratio for the different zirconia phases are compared in Table 5. We applied the present CT potential to predict the properties of orthorhombic zirconia (Table 4). The energy of orthorhombic zirconia is much higher than the other zirconia phases, so the orthorhombic phase is relatively unstable. This is consistent with the fact that orthorhombic ZrO₂ only exists under high pressure [26].

4.4 Dynamic properties

We then used the CT potential for MD simulations. The MD simulations were performed using the publicly available simulation code LAMMPS [36]. The OVITO package was used for visualization of the MD snapshots [39]. The standard Newton equations of motion were integrated in time using the velocity Verlet algorithm with a time step of 1.0 fs. The structure contained 2×2×2 symmetric unit cells. A larger structure of 4×4×4 symmetric unit cells was also simulated and similar results were obtained. Periodic boundary conditions were applied in all three directions.

We simulated the heating process of monoclinic zirconia. The change of the volume during this process is shown in Fig. 6(a). There was an abrupt decrease in the volume at around 747 K, which indicates a possible monoclinic–tetragonal phase transition considering the smaller volume of the tetragonal phase. The Zr coordination number increased from seven to eight after this phase transition [Fig. 6(b)], which further confirms the transition from the monoclinic phase (with seven-fold Zr coordination) to the tetragonal phase (with eight-fold Zr coordination). Note that the decrease of the coordination in the high temperature range is due to strong thermal vibration at high temperature, where some neighboring oxygen atoms vibrate to a distance far from the Zr atom.

To explore the structural transition in more detail, the atomic displacement caused by the phase transition is shown in Fig. 7. The color bar indicates the coordination number of each atom. It clearly shows the structure transition from the monoclinic phase to the tetragonal phase. The arrow on each atom represents its displacement induced by the phase transition. The oxygen atoms O^Ⅰ with three-fold coordination were seriously reconstructed by large displacements and eventually deformed to four-fold coordination.

It should be noted that the critical temperature for the monoclinic–tetragonal transition predicted by the present CT potential is lower than the experimental value [1]. Furthermore, there was no obvious tetragonal–cubic phase transition at high temperature, which occurs at approximately 2377 K in the experiment [1]. This is because cubic–tetragonal distortion is underestimated by the present CT potential, so cubic–tetragonal distortion cannot take effect at high temperature.

5 Comparison with existing potentials

In this section, the CT potential developed in the present work is compared to some existing empirical models for zirconia.

5.1 CT potential versus the core–shell model: Simulation efficiency

To quantify the efficiency of the CT potential, we performed comparative MD simulations with the CT potential and the widely used core–shell model. A series of zirconia systems were investigated by MD simulations with either the CT potential or the Coulomb and core–shell model. The simulation parameters of the CT potential have been described in previous sections. The parameters of the core–shell model were taken from a previous study [21]. In the core–shell model, each ion is constructed by a core and shell pair, which are connected by a spring. The ion’s charge and mass are divided into the core and shell. The relative motion of the shell with respective to the core mimics the polarization of the cation. Considering that the polarization process is very fast, a very small mass must be set for the shell so that the shell can follow the core’s motion in a quasi-adiabatic manner. In other words, the shell moves very fast so that it can rapidly respond to the motion of the core. To accurately simulate the shell’s fast motion, the time step \mathrm{\Delta }t used for numerical integration of Newton’s equations of motion must be relatively small. As a consequence, the computational cost of a MD simulation with the core–shell model is very expensive. The time step should be properly chosen in MD simulations. A small time step can reduce numerical errors in integration of Newton’s equations of motion, but it requires a long simulation time. In contrast, a large time step can reduce simulation time, but it will cause large numerical errors and may even lead to artificial collapse of the simulated structure. We introduce the critical value \mathrm{\Delta }{t}_c as the upper limit of the time step, above which the simulated structure will artificially collapse. The \mathrm{\Delta }{t}_c values for the time step in MD simulations with the CT potential and core–shell model are compared in Fig. 8. In this set of simulations, for a given time step, we ran MD simulations for 2×10⁴ steps and checked whether the structure collapsed. Zirconia does not melt for temperatures below 2500 K, so such collapse (if it occurs) is an artificial effect induced by a large time step. We can thus determine the critical value of the time step. We found that the critical time step decreases with increasing temperature because of stronger thermal vibration at higher temperature. A distinct feature is that the critical time step for the core–shell model is approximately one order of magnitude smaller than \mathrm{\Delta }{t}_c for the CT potential. For example, at 2000 K, the critical time step for the core–shell model is 0.35 fs, which is around 20 times smaller than the value of 7.25 fs for the CT potential. As a result of the extremely small time step, a much longer simulation time is needed for MD simulations using the core–shell model.

We directly compared the actual simulation times for the CT potential and core–shell model. We simulated several zirconia structures with total atom numbers ranging from 10² to 10⁴. All structures were simulated for 10 ps at room temperature. According to Fig. 8, the time step was chosen to be 1.0 fs and 0.05 fs for MD simulations with the CT potential and core–shell model, respectively. The simulation time was proportional to the size of the structure (Fig. 9). This shows that the CT potential is approximately one or two orders of magnitude faster than the core–shell model. The higher efficiency of the CT potential is mainly because an extremely small time step is not required, but it is required for the core–shell model.

5.2 CT potential versus the Born–Mayer potential: Surface reconstruction

Free surfaces have important properties for ionic crystals like zirconia. We applied the CT potential and Born–Mayer potential to study the typical free surfaces of zirconia. The parameters of the Born–Mayer potential were taken from a previous study [13]. The surface energy is a useful quantity that determines many physical properties of a free surface. The surface energy measures the energy excess caused by the appearance of a free surface. It can be calculated as follows:

\begin{array}{}\text{(17)}& \gamma =\frac{E_s-E_b}{A}\mathrm{,}\end{array}

where A is the total area of the two surfaces created, and E_b and E_s are the energy of bulk zirconia and the energy of zirconia with free surfaces, respectively.

We performed two calculations of each zirconia system with given dimensions of l_x\times l_y\times l_z. In the first calculation, bulk zirconia was mimicked by applying periodic boundary conditions in all directions. The bulk zirconia was relaxed by minimizing its potential energy with the conjugate gradient algorithm. The total potential energy gives the value of E_b. In the second calculation, two free surfaces perpendicular to the z axis were generated by applying the free boundary condition along the z direction, while periodic boundary conditions were applied in the x and y directions. The structure was first relaxed by minimizing the potential energy with the conjugate gradient algorithm. The free surface was not sufficiently relaxed by the energy minimization process because of the symmetry constraints in the conjugate gradient algorithm. The free surface was thus further optimized by MD simulation at a very low temperature of 4.2 K with the NPT (constant number of particles, pressure, and temperature) ensemble. The randomness in the MD simulation was able to assist relaxation of the free surface to its lowest energy configuration. An extremely low temperature was used so that the temperature-induced variation in the potential energy was negligible. The MD simulation was performed for 100 ps. The final potential energy gives E_s in Eq.17.

The thickness dependence of the surface energies calculated with the CT potential is shown in Fig. 10. The monoclinic phase is the most stable structure of zirconia at low temperature. We thus simulated monoclinic zirconia with the CT potential. The dimensions of zirconia were l_x=51.5 Å and l_y=52.0 Å, and the thickness l_z was increased from 5.2 Å to 131.2 Å. Note that l_x, l_y, and l_z are the sizes along the directions of the three lattice vectors in the monoclinic phase. The surface energies of the (001) and (100) surfaces are not sensitive to the thickness. Reconstruction of the (001) and (100) surfaces is weak (Fig. 11) because both surfaces are non-polarized. Approximately half of the neighboring atoms are removed for atoms within the surface region because of the appearance of the free surface. The energies of the atoms within the surface region will increase, which is the main origin of the surface energies of the (001) and (100) surfaces. The (010) surface energy decreases with increasing thickness and converges to a constant value at large thickness. This size dependence is closely related to surface reconstruction, as shown in Fig. 11(c), which is because of polarization of this surface. Surface reconstruction is more serious for thinner zirconia because of the larger surface to volume ratio. At large thickness, reconstruction is localized at the surface region, so the surface energy becomes constant at large thickness.

An overall feature of Fig. 10 is that the surface energy calculated with the CT potential depends on the thickness l_z, but it converges to a constant value at large thickness for all of the surfaces. This indicates the suitability of the CT potential in describing the free surfaces of monoclinic zirconia. These converged values are compared with the results from first-principles calculations in Table 7. The results obtained with the CT potential are generally larger than those obtained by first-principles calculations. However, the CT potential predicts that the (001) surface is the lowest energy surface, which is the same as first-principles calculations.

The thickness dependence of the surface energies calculated with the Born–Mayer potential is shown in Fig. 12. We used the Born–Mayer potential to calculate the surface energies of tetragonal zirconia, which is the most stable phase of zirconia using the Born–Mayer potential. The oxygen-terminated surface is energetically more favorable than the Zr or double oxygen layer (O–O) terminated surfaces [37, 40]. For the (100) and (001) surfaces, the ideal surfaces are polar and can be turned into neutral surfaces by removing half of the oxygen atoms of the outermost atomic layer [41]. The surface energies of the (111), (110), and (100) free surfaces converge with increasing thickness l_z. The (111) and (110) surfaces are non-polarized, so their surface energies do not depend on the thickness. There is a sharp increase in the surface energy of the (100) free surface around the critical thickness of 80 Å, below which the surface energy is sensitive to the thickness of zirconia. The underlying mechanism for this critical thickness is a transition from the tetragonal phase to the orthorhombic phase in thin zirconia, as shown in Fig. 13(a). For thick zirconia, the energy gained from surface reconstruction is not able to induce the tetragonal–orthorhombic phase transition, so reconstruction is localized in the surface region [Fig. 13(b)]. The surface energies of the (101) and (001) surfaces do not converge at large thickness. This is because reconstruction of the two surfaces is so strong that a tetragonal–orthorhombic transition is induced (Fig. 14). This is different from previous first-principles calculations, which predicted convergence of these two surface energies. The surface energies obtained with the Born–Mayer potential are compared to the results from first-principles calculations in Table 8. The surface energies obtained with the Born–Mayer potential are generally larger than those obtained by first-principles calculations.

6 Conclusion

Before concluding, we will make some general remarks about the positive and negative features of the present CT potential to help readers decide whether the CT potential is suitable for their research. The most significant feature of the CT potential is substitution of the bond-order Tersoff potential for the core–shell model. Both positive and negative features of the CT potential directly result from this substitution. We list these positive and negative features below.

Positive features. (i) The CT potential predicts the correct energy order of the four zirconia phases, with the monoclinic phase as the lowest energy structure. A clear monoclinic–tetragonal phase transition is observed in MD simulation. (ii) The CT potential is at least one order of magnitude faster than the core–shell model in MD simulation. This is because, to mimic an adiabatic response of the shells to the cores, the shells must have a very small mass and thus require a very small time step in MD simulation. (iii) The Tersoff potential is widely used in the computational community and it has been implemented in most simulation packages, so the CT potential can be conveniently used.

Negative features. (i) The magnitude of the cubic–tetragonal distortion is weaker than experiment, and consequently the tetragonal–cubic phase transition is not observed in MD simulation. (ii) The polarization effect is effectively treated by the Tersoff potential without introducing shells, so the effect of an external electric field on polarization of the shells cannot be simulated by the CT potential.

In summary, we have developed a potential by combining the Coulomb interaction and Tersoff potential to describe the atomic interactions of zirconia. The bond-order property of the Tersoff potential means that this potential is suitable for describing the well-known zirconia phases. In particular, this potential predicts that monoclinic zirconia is the most stable phase in the low-temperature region. The cubic phase is not stable and spontaneously transforms to the tetragonal phase. The orthorhombic phase has the highest energy and only exists under high pressure. These predictions agree reasonably well with experiments or ab initio calculations. We also used this potential to predict various static and dynamic properties of the zirconia phases. The potential scripts for GULP and LAMMPS are available from the website of the corresponding author (jiangjinwu.org).

References

[1] E. H. Kisi and C. J. Howard, Crystal structures of zirconia phases and their inter-relation, Key Eng. Mater. 153-154, 1 (1998)

[2] N. P. Padture, M. Gell, and E. H. Jordan, Thermal barrier coatings for gas-turbine engine applications, Science 296(5566), 280 (2002)

[3] E. H. Kisi and C. J. Howard, in: Zirconia Engineering Ceramics: Old Challenges-new Ideas, Netikon-Zurich: Trans Tech, 1998

[4] H. J. F. Jansen and J. A. Gardner, Total energy calculations for ZrO₂, Physica B 150(1–2), 10 (1988)

[5] K. Parlinski, Z. Q. Li, and Y. Kawazoe, First-principles determination of the soft mode in cubic ZrO₂, Phys. Rev. Lett. 78(21), 4063 (1997)

[6] G. Jomard, T. Petit, A. Pasturel, L. Magaud, G. Kresse, and J. Hafner, First-principles calculations to describe zirconia pseudopolymorphs, Phys. Rev. B 59(6), 4044 (1999)

[7] A. Kuwabara, T. Tohei, T. Yamamoto, and I. Tanaka, Ab initio lattice dynamics and phase transformations of ZrO₂, Phys. Rev. B 71(6), 064301 (2005)

[8] P. Souvatzis and S. P. Rudin, Dynamical stabilization of cubic ZrO₂ by phonon-phonon interactions: Ab initio calculations, Phys. Rev. B 78(18), 184304 (2008)

[9] H. Wu, Y. Duan, K. Liu, D. Lv, L. Qin, L. Shi, and G. Tang, First-principles study of phase transition and band structure of ZrO₂ under pressure, J. Alloys Compd. 645, 352 (2015)

[10] C. W. Li, H. L. Smith, T. Lan, J. L. Niedziela, J. A. Munoz, J. B. Keith, L. Mauger, D. L. Abernathy, and B. Fultz, Phonon anharmonicity of monoclinic zirconia and yttrium-stabilized zirconia, Phys. Rev. B 91(14), 144302 (2015)

[11] G. V. Lewis and C. R. A. Catlow, Potential models for ionic oxides, J. Phys. C 18(6), 1149 (1985)

[12] M. Born and K. Huang, Dynamical Theory of Crystal Lattices, Oxford: Oxford University Press, 1954

[13] P. K. Schelling, S. R. Phillpot, and D. Wolf, Mechanism of the cubic-to-tetragonal phase transition in zirconia and yttria-stabilized zirconia by molecular-dynamics simulation, J. Am. Ceram. Soc. 84(7), 1609 (2001)

[14] M. Kilo, C. Argirusis, G. Borchardt, and R. A. Jackson, Oxygen diffusion in yttria stabilised zirconia — experimental results and molecular dynamics calculations, Phys. Chem. Chem. Phys. 5(11), 2219 (2003)

[15] C. Yang, K. Trachenko, S. Hull, I. T. Todorov, and M. T. Dove, Emergence of microstructure and oxygen diffusion in yttrium-stabilized cubic zirconia, Phys. Rev. B 97(18), 184107 (2018)

[16] J. B. G. Dick and A. W. Overhauser, Theory of the dielectric constants of alkali halide crystals, Phys. Rev. 112(1), 90 (1958)

[17] A. Dwivedi and A. N. Cormack, A computer simulation study of the defect structure of calcia-stabilized zirconia, Philos. Mag. A 61(1), 1 (1990)

[18] M. Wilson, U. Schonberger, and M. W. Finnis, Transferable atomistic model to describe the energetics of zirconia, Phys. Rev. B 54(13), 9147 (1996)

[19] K. C. Lau and B. I. Dunlap, Molecular dynamics simulation of yttria-stabilized zirconia (YSZ) crystalline and amorphous solids, J. Phys.: Condens. Matter 23(3), 035401 (2011)

[20] F. Shimojo, T. Okabe, F. Tachibana, M. Kobayashi, and H. Okazaki, Molecular dynamics studies of yttria stabilized zirconia (I): Structure and oxygen diffusion, J. Phys. Soc. Jpn. 61(8), 2848 (1992)

[21] M. Smirnov, A. Mirgorodsky, and R. Guinebretiere, Phenomenological theory of lattice dynamics and polymorphism of ZrO₂, Phys. Rev. B 68(10), 104106 (2003)

[22] S. Fabris, A. T. Paxton, and M. W. Finnis, Relative energetics and structural properties of zirconia using a self-consistent tight-binding model, Phys. Rev. B 61(10), 6617 (2000)

[23] A. C. T. van Duin, B. V. Merinov, S. S. Jang, and W. A. Goddard, The ReaxFF reactive force field for solid oxide fuel cell systems with application to oxygen ion transport in yttria-stabilized zirconia, J. Phys. Chem. A 112, 3133 (2008)

[24] C. Wang, A. Tharval, and J. R. Kitchin, A density functional theory parameterised neural network model of zirconia, Mol. Simul. 44(8), 623 (2018) https://doi.org/10.1080/08927022.2017.1420185

[25] J. Tersoff, New empirical model for the structural properties of silicon, Phys. Rev. Lett. 56(6), 632 (1986)

[26] O. Ohtaka, H. Fukui, T. Kunisada, T. Fujisawa, K. Funakoshi, W. Utsumi, T. Irifune, K. Kuroda, and T. Kikegawa, Phase relations and equations of state of ZrO₂ under high temperature and high pressure, Phys. Rev. B 63(17), 174108 (2001)

[27] D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r–1 summation, J. Chem. Phys. 110(17), 8254 (1999)

[28] C. J. Fennell and J. D. Gezelter, Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics, J. Chem. Phys. 124(23), 234104 (2006)

[29] S. Dai, M. Gharbi, P. Sharma, and H. S. Park, Surface piezoelectricity: Size effects in nanostructures and the emergence of piezoelectricity in non-piezoelectric materials, J. Appl. Phys. 110(10), 104305 (2011)

[30] R. Agrawal, B. Peng, E. E. Gdoutos, and H. D. Espinosa, Elasticity size effects in ZnO nanowires: A combined experimental-computational approach, Nano Lett. 8(11), 3668 (2008)

[31] J. Tersoff, New empirical approach for the structure and energy of covalent systems, Phys. Rev. B 37(12), 6991 (1988)

[32] J. Tersoff, Modeling solid-state chemistry: Interatomic potentials for multicomponent systems, Phys. Rev. B 39(8), 5566 (1989)

[33] P. M. Morse, Diatomic molecules according to the wave mechanics (II): Vibrational levels, Phys. Rev. 34(1), 57 (1929)

[34] W. A. Harrison, Elemetary Electronic Structure, Singapore: World Scientific, 2004

[35] J. D. Gale, GULP: A computer program for the symmetry-adapted simulation of solids, J. Chem. Soc. Faraday Trans. 93(4), 629 (1997)

[36] S. J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117(1), 1 (1995)

[37] A. Christensen and E. A. Carter, First-principles study of the surfaces of zirconia, Phys. Rev. B 58(12), 8050 (1998)

[38] A. Kokalj, Computer graphics and graphical user interfaces as tools in simulations of matter at the atomic scale, Comput. Mater. Sci. 28(2), 155 (2003)

[39] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO — the Open Visualization Tool, Model. Simul. Mater. Sci. Eng. 18(1), 015012 (2010)

[40] A. Eichler and G. Kresse, First-principles calculations for the surface termination of pure and yttria-doped zirconia surfaces, Phys. Rev. B 69(4), 045402 (2004)

[41] G. Ballabio, M. Bernasconi, F. Pietrucci, and S. Serra, Ab initio study of yttria-stabilized cubic zirconia surfaces, Phys. Rev. B 70(7), 075417 (2004)

Rights & permissions

Higher Education Press 2021