Ferroelastic switching of doped zirconia : Modeling and understanding from first principles

The properties of materials at high temperatures are often determined by complex thermodynamic mechanisms. One of the most prominent examples is the stabilization of tetragonal and cubic zirconia, which we investigate using density functional theory. The results show that the minimum energy path for the tetragonal-to-cubic phase transformation differs significantly from the paths discussed in the literature so far. This provides insight into the properties of compositions codoped with yttria and titania, an approach that has recently been proposed for the design of thermal barrier coatings.

Zirconia-based materials offer many appealing properties, including a low thermal conductivity and a high ionic conductivity, as well as the potential for remarkable toughness.Hence, they are employed in a wide variety of applications, e.g., as catalyst supports, ionic conductors in sensors and solid oxide fuel cells, and thermal barrier coatings (TBC) for gas turbines in propulsion and power generation.In the latter case, a 150-1000 μm thick zirconia-based coating is applied to the turbine's superalloy components to protect them from the extreme temperatures generated during combustion.This improves the durability of the components and also enables an increase of the operation temperature and fuel efficiency [1].
Pure ZrO 2 is not suitable for most applications because of the monoclinic-tetragonal phase transition at approximately 1500 K, which involves a disruptive volume change (∼4%) that degrades its mechanical integrity [1].Thus, most ZrO 2 materials are doped to control or suppress the transformation, as it is the case for state-of-the-art TBCs based on singlephase tetragonal zirconia stabilized with 8 ± 1 mol % YO 1.5 (8YSZ).The remarkable durability of 8YSZ is ascribed to its superior toughness [2] based on a ferroelastic domain switching mechanism [1,3].However, 8YSZ is metastable as a single tetragonal phase and eventually decomposes into Y-rich cubic and Y-lean tetragonal domains, the latter susceptible to the disruptive monoclinic transformation [4].Accordingly, important goals in the development of advanced TBCs are phase stability at the higher temperatures (>1400 K) and improved toughness to mitigate a host of potential damage mechanisms [1].For this purpose, it is essential to understand the fundamental origins of the mechanisms underlying the stabilization and toughening of the tetragonal phase.
A common approach to stabilize the higher temperature forms of ZrO 2 involves doping with trivalent cations, most commonly Y 3+ , with the concomitant introduction of charged oxygen vacancies (F 2+ ) [5,6].For TBCs, the useful composi-tional range is very narrow (8 ± 1 mol % YO 1.5 ): Underdoping renders the structure transformable to monoclinic on cooling, while overdoping reduces the tetragonality and hence the potential for ferroelastic toughening [1,3].Indeed, excess doping eventually leads to stabilization of the cubic form with a much lower toughness [2].In contrast thereto, recent work involving doping with Y 3+ and Ti 4+ has led to tetragonal alloys with superior phase stability and improved toughness [3].The present paper sheds light on the effects of doping on the electronic and atomistic mechanisms underpinning the tetragonal phase stability and the energetic barriers for domain switching.
We performed density functional theory (DFT) calculations of pure and doped zirconia with the full potential, numeric atomic orbital code FHI-aims [7].To correctly describe the breakdown of short-range periodicity upon doping, the calculations were performed for 2 × 2 × 2 supercells (96 atoms).Equilibrium configurations are determined by relaxing both the atomic positions and the unit-cell shape; energetic barriers and respective minimum energy paths (MEPs) are computed by means of the climbing-image-nudged-elasticband (CI-NEB) method [8] and its generalization for solids (CI-G-SSNEB) [9] that additionally accounts for unit-cell degrees of freedom.The employed numerical settings [10] ensure a numerical error smaller than 5 meV/atom for energy differences.
Since electronic exchange and correlation effects play an important role for the relative stability of different zirconia polymorphs [6], it is important to validate our approach by comparing calculations performed at different levels of theory.For this purpose, we here discuss the influence of the exchangecorrelation functional on the fully optimized geometries and total energy differences for monoclinic (baddeleyite structure, P2 1 /c), tetragonal (P4 2 /nmc), and cubic (fluorite structure, Fm3m) zirconia, as summarized in Table I.In spite of the fact that all approaches yield the same qualitative trends, we find that the various local (LDA [11]) and semilocal (PBEsol [12], PBE [13], RPBE [14]) DFT exchange-correlation functionals predict energy differences E m,t and E t,c that differ almost by a factor of three.In particular, we find that the LDA functional, which has been used extensively in studies of ZrO 2 [15,16], severely underestimates E m,t and E t,c , whereas the PBEsol results lie closest to the outcome of higher level, As shown in Fig. 1, the zirconium ions in tetragonal, pure ZrO 2 form a slightly elongated fcc crystal (c/a t > 1) with eight oxygen anions that are pairwise displaced (dz > 0) from the tetrahedral site in an antiparallel fashion along the c axis.It is important to realize that cubic zirconia (upper ball-and-stick model in Fig. 2), which corresponds to the special case of a tetragonal configuration with c/a c = 1 and dx,dy,dz = 0, does not correspond to an energetic minimum, but to a saddle point in the energetic landscape, as demonstrated before [15,16].As the potential-energy curve in Fig. 2(a) shows, the slightest perturbation leads to a spontaneous relaxation into one of the minima, i.e., in one of the six (a pair for each cartesian direction) symmetry-equivalent structures with nonvanishing tetragonal displacements.The cubic structure thus corresponds to the transition state along the minimum energy path that connects the tetragonal structures in their primitive unit cell [also see Fig. 1(a)].
The potential-energy curve shown in Fig. 2(a) has often been used to rationalize the dynamics in tetragonal zirconia, in particular the second order, tetragonal-to-cubic phase transition that occurs at high temperatures 2500 K (dynamical stabilization) in pure ZrO 2 [15,16].At the LDA level of theory used in these studies, the free energy surface becomes flat at these temperatures, which implies that the oxygen atoms can overcome the barrier E t,c and that the thermodynamic average of dz becomes zero [15].Similarly, various theoretical stud-  ies [5,18] have inspected changes in this potential-energy curve and in E t,c as a function of doping to explain the decrease in the tetragonal-to-cubic phase transition temperature (∼1700 K in 8YSZ [1]) and the stabilization of cubic zirconia at high levels of doping (>16 mol % YO 1.5 ).
Compared to these earlier studies (cf. the LDA-based model [5,15] described in the previous paragraph) our results imply a significantly, even qualitatively, different dynamics of pure and doped zirconia, since the higher value of E t,c at the PBEsol level of theory cannot be easily overcome at temperatures 3000 K: As the potential-energy surfaces (PESs) in Figs.3(a) and 3(a ) reveal, the actual MEP does not go along the cubic (dx,dy,dz = 0) configuration if tetragonal displacements dx and dy are allowed.Rather, a realignment of the geometry along a different cartesian direction takes place.Due to this realignment of the tetragonal displacement (and of the tetragonality in the CI-G-SSNEB calculation), the dz = 0 configuration becomes a (local) minimum of the MEP and the respective energy barrier E b is halved with respect to E t,c [cf.Fig. 2(b)].In particular, this peculiar MEP provides critical insight in the mechanism that governs ferroelastic switching, i.e., the collective reorientation of both the anions and the lattice vectors that hitherto lacked an atomistic description and explanation [2].
We find direct evidence of these ferroelastic switches in ab initio molecular dynamics (aiMD) calculations performed for a 96-atom supercell: Thermal fluctuations are significant, but we can track the time evolution of the average tetragonal displacements dx , dy and dz by backfolding the geometries of the supercell to the primitive unit cell and averaging over all oxygen atoms.Both in the tetragonal [cf. at low temperatures (900 and 300 K), since dx , dy and dz just oscillate around their respective equilibrium values.At increased temperatures (1500 and 600 K), spontaneous changes in the average tetragonal displacements occur, i.e., from dz to − dz in Fig. 4(b) and from dz to dy in Fig. 4(e).For even higher temperatures (2400 and 1200 K), the oscillations grow so strong that the time averages of the tetragonal displacements dx , dy and dz vanish, as shown in Figs.4(c) and 4(f).
In all these cases (a)-(f), however, the respective radii dr = dx 2 + dy 2 + dz 2 , i.e., the distances from the highsymmetry cubic configuration averaged over all oxygens, are finite and the probability P ( dr ) that the cubic Fm3m structure ( dr = 0) occurs is vanishingly small.This indicates that the free energy surface F ( dr ) = −k B T ln (P ( dr )) is still very corrugated, since the cubic structure corresponds to a local maximum of F ( dr ) at all temperatures below the melting point.In spite of that, the spontaneous realignment of the tetragonal displacements results in an apparent cubic structure at high temperatures, since the time average of dx , dy , dz vanishes [20].Under these thermodynamic conditions the system has no preferential axis and hence exhibits cubic fluorite Fm3m symmetry on average.
Due to the fact that earlier insightful studies of pure zirconia by Finnis et al. relied on a LDA-based tight-binding Hamiltonian [15], the relevance of the MEP discussed above and the implications to the understanding of ferroelastic switching could not be realized before.Similarly, studies of doped zirconia [5,18,21,22] did not investigate the relevant transition state associated to this MEP.We illustrate the importance of these effects by comparing the PES and MEP of (a) pure zirconia with those of (b) zirconia with a F 2+ vacancy (VSZ) [5,23], (c) 6.25 mol % YO 1.5 doped ZrO 2 (6YSZ) and (d) 6.25 mol % YO 1.5 plus 3.125 mol % TiO 2 doped ZrO 2 (3Ti6YSZ).By fully optimizing all possible tetragonal and monoclinic defect configurations, we identified various (nearly-)degenerate geometries.All associated MEPs show the same qualitative features.Thus, in the following we limit the discussion to the most stable geometry for each compositions, i.e., the one with titanium in the first and yttria in the second coordination shell of the vacancy, as found in previous studies [18,21,22].As a generalized coordinate for plotting of the MEPs in Fig. 5, we use again dz , i.e., the tetragonal displacement averaged over all oxygen atoms, since dz itself is no longer equal for all oxygen atoms due to the local relaxations induced by the dopants.For the exact same reason, no meaningful two-dimensional representation of the PES can be plotted by using the same value of dx, dy and dz for all oxygens.To still give qualitative insight on the topology of the PES, we first determined the respective equilibrium geometries r ±dz at ±dz and then used r 0 = r +dz + r −dz , r z = r +dz − r −dz and the respectively rotated r y to map out the PES as a function of the average tetragonal displacements dy and dz in Figs.In the case of VSZ, the MEP in Fig. 5 shows that the vacancy reduces dz (c/a t = 1.017) and lowers E m,t and E b and thereby stabilizes the tetragonal (with respect to the monoclinic) structure and reduces the energy barrier to domain switching.Since the vacancy also reduces E t,c , previous studies focusing on the erroneous MEP were able to find a similar qualitative trend [5].The topology of the PES, however, remains unchanged by the F 2+ center [cf.Fig. 3(b)].Conversely, the additional presence of yttria anions in 6YSZ breaks the symmetry, as a result of which the degeneracy of the ± dz states is lifted.Also, the yttria further lowers E m,t but hardly affects E b and hence further stabilizes the tetragonal structure (c/a t = 1.014), but does not affect the barrier for domain switching with respect to VSZ.In 3Ti6YSZ, the additional titanium induces strong local tetragonal displacements in its first coordination shell dz > 0.06 that lead to an increase in c/a t = 1.015.Since similar relaxations occur in the monoclinic structure, Ti doping hardly affects E m,t with respect to 6YSZ; the respective MEP and the topology of the PES, however, change significantly: The addition of titania increases the energetic difference between the different switched configurations, since TiYSZ strongly favors geometries with a large vacancy-titania distance (+ dz ), so that also the barrier for ferroelastic switches increases dramatically to a value that is even higher than in pure ZrO 2 (cf.Fig. 5).Given that VSZ, 6YSZ, and 3Ti6YSZ contain the same amount of vacancies, this reveals that the dynamics of doped ZrO 2 cannot be rationalized solely in terms of oxygen vacancies as suggested by previous studies [5].
In summary, our calculations show that the minimum energy path associated with the tetragonal-to-cubic phase transformation in pure zirconia does not go across the high-symmetry cubic configuration, but rather involves a spontaneous and collective realignment of the tetragonal geometry along a different cartesian direction.In thermodynamic average, this ferroelastic switching dynamics leads to an apparent cubic Fm3m structure that is, however, never realized as a stable or metastable structure in pure ZrO 2 .Recent advancements in pulsed, high-intensity electron and light sources [24,25] should enable an experimental examination of our findings (and of the computed barriers), given that the tetragonal displacements of the anions can be detected by vibrational spectroscopy [26] or electron diffraction [27].Similar mechanisms may also be active in other materials, given that soft vibrational modes and tetragonal-cubic phase transitions are not a unique feature of zirconia.Even more importantly, the discussed dynamics supports the models that ascribe the durability of TBCs to a remarkable toughening mechanism based on ferroelastic switching: The reorientation of tetragonal domains under the stress field of a crack yields a residual stress field that shields the crack tip and increases the work of fracture [2,3].Our calculations show that doping with YO 1.5 stabilizes the tetragonal structure, but also lifts the degeneracy between the various tetragonal orientations, so that ferroelastic switches may become a viable mechanism to absorb energy during crack propagation.Codoping with titania induces strong local displacements that further increase the energetic difference between differently oriented configurations and that increase the energetic barrier for ferroelastic switches, which in turn would translate into a beneficial increment in fracture energy.These results are consistent with the experimental finding that codoping of 8YSZ with titania gives rise to an increase in toughness [1,3].Accordingly, our calculations provide insights in the underlying atomistic mechanism that go well beyond the geometrical considerations used in ferroelastic toughening models so far [2].In particular, we show that the energetics of domain switching can be tailored by cation codoping and thus lay the foundation for further research along these lines, e.g., in multiscale models.

FIG. 1
FIG. 1. (Color online) (a) Primitive six-atom unit cell of tetragonal ZrO 2 with a tetragonal displacement dz.Displacements dx,dy orthogonal to the c axis are impossible, since the periodic images of the oxygen atoms (orange) are interlinked in this plane via the lattice vectors.Such displacements become feasible in the 12-atom unit cell, as shown in the ball-and-stick model (b).
FIG. 2. (Color online) Calculated (DFT-PBEsol) MEPs (in eV per 12-atom unit cell) between two symmetry-equivalent tetragonal configurations in pure ZrO 2 .The black lines ( ) denote CI-NEB calculations performed in the tetragonal equilibrium lattice (cf.Table I), whereas the red lines ( ) denote CI-G-SSNEB calculations, in which also the unit-cell shape is optimized along the path.Symbols mark the "images" used in the actual calculations.Panel (a) shows the MEP in the primitive tetragonal unit cell, panel (b) the MEP in nonprimitive supercells that allow for tetragonal displacements dx and dy.The respective initial, intermediate (dz = 0), and final geometries are shown, where blue arrows mark the tetragonal displacement in the first oxygen layer.
FIG. 3. (Color online) DFT-PBEsol potential-energy surfaces (E − E t in eV per 12-atom ZrO 2 unit cell, interpolated from an equally spaced 17 × 17 grid) as a function of the (average) tetragonal displacements dy, dz ( dy , dz ) for (a) pure Zr 4 O 8 with dx = 0, (b) VSZ (Zr 32 O 63 F 2+ ), (c) 6YSZ (Y 2 Zr 30 O 63 ), and (d) 3Ti6YSZ (TiY 2 Zr 29 O 63 ).In the doped cases (b)-(d), in which dz is no longer equal for all oxygen atoms due to local relaxations, we show specific two-dimensional cuts of the full PESs (see text).The green arrow in panel (a) denotes the MEP shown in Fig. 2(a), whereas the dotted orange lines hint at the MEPs shown in Figs.2(b) and 5.The relaxation of the lattice vectors was accounted for only in the plots (a )-(d ).
3(b)-3(d).For Figs. 3(a )-3(d ), a relaxation of the lattice vectors was performed at each point while keeping the fractional coordinates fixed.

FIG. 5 .
FIG. 5. (Color online) Computed (DFT-PBEsol) MEPs (eV per 12-atom ZrO 2 unit-cell) as a function of the average tetragonal displacement dz for (a) pure Zr 4 O 8 , (b) VSZ, (c) 6YSZ, and (d) 3Ti6YSZ.The left plot shows MEPs in the respective tetragonal equilibrium lattices (CI-NEB), whereas the relaxation of the lattice vectors was taken into account in the right plot (CI-G-SSNEB).Open symbols mark the "images" of the NEB, full symbols denote the transition states including their energy barriers E b and dashed lines denote E m,t .

TABLE I .
Equilibrium zero Kelvin lattice constants a c and a t , tetragonalities c/a t , tetragonal displacements dz (cf.Fig.1), and energy differences per 12-atom cell E t,c = E c − E t and E m,t = E t − E m for fully optimized cubic (c), tetragonal (t), and monoclinic (m) zirconia at different levels of theory.
[17]utationally much more involved calculations with the hybrid functional HSE06[17].Since PBEsol yields almost the same geometries as HSE06, we limit ourselves to present PBEsol calculations in what follows.