Oxygen Vacancies in Zirconia and Their Migration: The Role of Hubbard-U Parameters in Density Functional Theory

Abstract

Cubic zirconia (c-ZrO${}_2$) is studied using Density Functional Theory with Hubbard-U corrections (DFT+U). It is shown that the determination of the U-parameters from first principles leads to values for U(Zr-4d) and U(O-2p) which are very different from standard choices. The calculated band gap with these values for U closely matches the experimental gap. Oxygen vacancies are studied using this approach, and it is found that it is possible to closely reproduce the vacancy migration energies calculated with a hybrid functional. The oxygen vacancy is associated with two excess electrons which localize in the vacancy’s cavity. In the presence of these excess electrons, the barrier for vacancy migration is very high. If instead, a charged vacancy V${}_{\mathrm{O}}$${}^{2+}$ is considered, its mobility increases considerably—an effect that is attributed to the absence of space charges localized in the cavity.

1. Introduction

Zirconia (ZrO${}_2$) is a transition metal oxide which has attracted considerable interest in recent years [1,2,3,4]. Its mechanical strength, chemical stability and high dielectric constant make it the material of choice for many applications. At ambient pressure, zirconia exists in three different phases: the monoclinic phase (m-ZrO${}_2$) is stable up to a temperature of 1478 K. At intermediate temperatures, from 1478 K to 2650 K, the tetragonal phase (t-ZrO${}_2$) is stable, while at higher temperatures, up to the melting temperature of 2983 K, zirconia exists in its cubic phase (c-ZrO${}_2$) [1]. Transitions between these phases imply a rather important change of volume, which is unsuitable for many applications with large temperature variations. For this reason, attention has been focused on the possibilities to stabilize the cubic phase at low temperatures. It was recognized first by simulations [5] and later experimentally [4] that oxygen vacancies (V${}_{\mathrm{O}}$) allow for the stabilization of c-ZrO${}_2$ at room temperature. Concentrations of V${}_{\mathrm{O}}$ as low as 3% are shown to be sufficient to stabilize the cubic phase at low temperatures.

Oxygen vacancies are of central importance not only to stabilize the desired phase of ZrO${}_2$. During oxide growth via corrosion, the oxide forms an overlayer between pure zirconium metal and an oxygen-rich environment. In such a setup, the oxide grows by forming ZrO${}_2$ at the metal–oxide interface. In this process, zirconia structures with oxygen vacancies are formed. Those V${}_{\mathrm{O}}$ migrate through the ZrO${}_2$ film towards the surface where they are annihilated by incorporation of oxygen from the environment [2]. Clearly, the existence and energetics of oxygen vacancies and their migration properties are therefore an important concept to understand growth mechanisms of zirconia, and also many of its properties in applications [6,7,8,9,10,11].

Computer simulations are one important way to elucidate the underlying elementary processes (such as vacancy migration) of a material at the atomic scale. Treating metal oxides by means of the Density Functional Theory (DFT) has therefore attracted a lot of attention in recent years. Unfortunately, the strongly localized character of electrons in transition metals’ d-orbitals leads to the so-called self-interaction error (SIE) inherent in the most commonly used density functionals. Semi-local approximations to the exchange and correlation energy fail to cancel the spurious Hartree interaction of an electron with itself. The importance of the SIE has long been recognized, and there are different approaches to cancel or mitigate this error: one possibility is the use of SIE-corrected functionals. In this approach, the electrons’ interactions with themselves are subtracted from the energy. The price to pay for this correction is that the resulting expression for the energy depends explicitly on the single-particle orbitals; thus, one has to deal with a more cumbersome minimization procedure where orbitals need to be kept orthogonal to each other by constraints. A second approach is to use hybrid functionals. In Hartree–Fock (HF) theory, the exchange energy exactly cancels the spurious self-interaction of the Hartree energy; thus, a hybrid functional, which contains a fraction of HF exchange, has less SIE. The disadvantage of hybrid functionals is their huge computational demand which often renders hybrid DFT calculations untractable practically.

A third possibility is the introduction of so-called Hubbard-U parameters. In DFT+U calculations, a term is added to the energy which penalizes fractional occupations of localised atomic states. Originally proposed for 3d-electrons of transition metals, such Hubbard-U terms can also be added to other localized states, such as, for example, 2p-states of oxygen. By penalizing fractional occupations, DFT+U approaches correct one important feature of self-interaction, namely the tendency to delocalize electronic orbitals. Moreover, the penalty on fractional occupation numbers also restores the derivative discontinuity of the system’s total energy as a function of electron number. In many cases, the DFT+U approach allows DFT calculations in materials which would otherwise have been wrongly described as a metal or as having a small band gap, transition metal oxides being a prime example [12,13,14,15,16]. As the additional computational cost of DFT+U calculations is low, this approach has quickly led to a widespread adoption of this technique.

Zirconia is an interesting case for DFT+U: the formal oxidation state of Zr in ZrO${}_2$ is +4, which means that the 4d states have an occupation of zero: 4d${}^0$. One could therefore assume that no SIE is present given that the localized orbitals are not occupied. Furthermore, ZrO${}_2$ has a wide band gap (3.3 eV) also at the semilocal PBE level of theory, albeit smaller than the experimental band gap of 5.7 eV. It is probably for these reasons that zirconia is often treated using PBE, without any corrections for the SIE [2].

In this paper, we re-visit the case of ZrO${}_2$ and its description by DFT+U methods, with a particular focus on oxygen vacancies and their mobility, which, as discussed above, is a key property for this material.

2. Computational Methods

We use a plane–wave basis set and pseudopotentials, as implemented in the quantum-ESPRESSO suite of programs [17,18]. The density functional is approximated using a semi-local gradient approximation, as parametrized by Perdew, Burke and Ernzerhof (PBE) [19]. We use the following pseudopotentials Zr.pbe-spn-rrkjus_psl.1.0.0.UPF and O.pbe-n-rrkjus_psl.1.0.0.UPF for Zr and O, respectively. The Zr pseudopotential treats explicitly 12 electrons as valence electrons, generated in the configuration $4s^{2}4p^{6}5s^{2}4d^2$, while the oxygen pseudopotential has a valence of 6 electrons, generated in the configuration $2s^{2}2p^4$. Both pseudopotentials have originally been developed for the pslibrary version 1.0.0 [20]. Regarding the plane–wave basis set, we used a kinetic energy cutoff of 40Ry for the Kohn–Sham orbitals and of 400 Ry for the charge density. The convergence of the energies with respect to these parameters has been carefully verified and is shown in the Supplemental Information.

For calculations involving a Hubbard-U correction, we use orthogonalized atomic orbitals for the projections. Orthogonalized atomic orbitals avoid the double application of Hubbard-U corrections in those areas of space where orbitals centered on different atoms overlap.

For comparison, we also report some results obtained with hybrid functionals. We use the PBE0 functional [21,22] where the PBE exchange energy is mixed with Hartree–Fock exchange energy in a 3:1 ratio and used with the usual PBE correlation energy.

To model the c-ZrO${}_2$ system, we use a 2 × 2 × 2 supercell which contains a total of 96 atoms (the unit cell of c-ZrO${}_2$ contains 4 formula units). The lattice constant is 5.09 Å. The Brillouin zone of the supercell is sampled on a regular 4 × 4 × 4 grid of k-points.

We provide in the Supplemental Information convergence test of the energy with respect to the basis set size and to the k-point sampling in Supplementary Materials.

All systems considered here have a finite band gap. For this reason, no “electronic smearing” is needed for the occupation numbers of the Kohn–Sham levels.

The adjustable parameters for the DFT calculations are summarized in Table 1.

3. Results

3.1. Bulk c-ZrO₂

The oxidation state is a very useful concept to understand the charge distribution in an ionic system. It is based on the hypothetical assignment of an integer number of electrons to individual atoms. In our case of ZrO${}_2$, it consists of assigning the charge state −2 to oxygen and +4 to zirconium atoms. Obviously, a real crystal is never fully ionic, and an attempt to “count” electrons based on charge distributions, such as the Mulliken [23] or Bader [24] population analysis will not lead to integer oxidation numbers. However, the concept of oxidation states can be very successfully generalized to realistic systems with a mixed ionic and covalent nature. Sit and co-workers have shown [25] that the traditional oxidation states can be derived from ab initio calculations by projecting the molecular orbitals onto the relevant localized atomic states. The eigenvalues of the projected atomic occupation matrix are used to determine the oxidation number: whenever such an eigenvalue is close to unity, one adds one electron to the count for the atom.

Applying this concept of occupation matrices to ZrO${}_2$ at the PBE level of theory and projecting the Kohn–Sham orbitals onto the zirconium d-electrons, one finds the following eigenvalues of the 5 × 5 occupation matrix:

$\begin{array}{ccccc}0.135& 0.135& 0.262& 0.262& 0.262.\end{array}$

Here, one clearly notes the crystal field splitting in threefold $t_{2g}$ and twofold $e_g$ symmetries, but in the context of oxidation state, none of the eigenvalues is close to unity, and therefore, as expected, one can assign a zero occupation to the zirconium d-states: $d^0$. This fact is the main reason why PBE calculations (without Hubbard-U or other corrections for self-interaction) have often been carried out successfully in this material [2].

The main shortcoming of PBE-level calculations in bulk ZrO${}_2$ is the strong under-estimation of the band gap. While the experimental gap is 5.7 eV, a standard PBE approach leads to a gap of 3.3 eV in c-ZrO${}_2$. It is interesting to note that the hybrid functional HSE06 [26] also leads in c-ZrO${}_2$ to a band gap of 4.72 eV [27], again underestimated with respect to experiments. As shown in Figure 1a, the states at the valence band maximum have mainly an O-2p character, while the conduction band minimum is mainly made of Zr-4d states. The fact that even standard PBE calculations lead to a relatively large band gap is a justification for using uncorrected semilocal functionals in this material, unless one is explicitly interested in optical properties. Many other transition metal oxides do not allow such an approach as they are wrongly described as metals at this level of theory.

For these reasons, it is a common practice to add a Hubbard-U correction to PBE calculations, with the value of the U-parameter fixed by the requirement to reproduce the experimental band gap. One such case is shown in panel (b) of Figure 1, where such a correction is added to the Zr-4d and the O-2p states, i.e., to those states which mainly make up the frontier orbitals. Clearly, the PBE+U approach is able to reproduce the experimental band gap.

At this point, the question arises whether the band gap alone is enough to determine suitable values for the Hubbard-U parameters. To answer this question, we have calculated the band gap as a function of the two mentioned parameters, the Us for the O-2p and for the Zr-4d sates. The result is shown in Figure 2. The black line in the figure corresponds to the experimental band gap. It is evident that there is a whole range of parameter combinations where the band gap is reproduced well: from very high values of U(Zr-4d) (around 9.5 eV) when no correction to the O-2p states is applied, to lower values for U(Zr-4d) when a correction is also applied to the oxygen atomic orbitals.

The ambiguity of the choice of U which is manifest in Figure 2 can be resolved in various ways. One possibility is to simply choose a second property to be fitted. Such an approach would, however, still depend on the choice of properties. Nevertheless, such an empirical fitting of U is still a very common practice. A second approach, free of any such arbitrariness and more in line with the concept of a parameter-free calculation from the first principles, is to calculate the values of the Hubbard parameters. This idea was first proposed by Cococcioni and de Gironcoli [28] and was recently implemented in a different but equivalent formulation in quantum-ESPRESSO [29]. The basic procedure is to use a linear-response formalism to determine values of the parameters U. Applied to the case of c-ZrO${}_2$, one finds the values U(Zr-4d) = 2.85 eV and U(O-2p) = 9.16 eV. The DOS of the system using these values is shown in panel (b) of Figure 1, where one can clearly see that the band gap is reproduced well by this non-empirical choice of the parameters. This is also evident in Figure 2, where these two values of U are close to the black line corresponding to the experimental band gap. It should be noted that such a relatively low value of U for the zirconium d-states and a high value for the oxygen p-states is quite uncommon. For example, by semi-empirically fitting these parameters to properties, the values of U(Zr-4d) = 8.0 eV and U(O-2p) = 4.35 eV have been used in Refs. [1,30] for m-ZrO${}_2$. In addition, this choice of parameters lies close to the black line in Figure 2, meaning that it is a suitable choice for U to reproduce the band gap in c-ZrO${}_2$.

3.2. Oxygen Vacancies in c-ZrO₂

As discussed in the Introduction, oxygen vacancies play a crucial role for the stability of c-ZrO${}_2$ at room temperature and for the growth mechanisms of the material. For this reason, we have studied the energetics of the migration of V${}_{\mathrm{O}}$. The oxidation number of oxygen in ZrO${}_2$ is −2. Therefore, removing a neutral oxygen atom from the crystal will lead to two excess electrons in the defected crystal. In order to examine the role of those excess electrons, we also considered the case where one neutral oxygen atom and two electrons have been removed from the crystal, leading to a charged vacancy V${}_{\mathrm{O}}$${}^{2+}$.

Using the values for the Hubbard-U parameters found by linear response (see Section 3.1 above), in Figure 3 we show a plot of the difference in charge density between the V${}_{\mathrm{O}}$ and V${}_{\mathrm{O}}$${}^{2+}$ systems. As the contour plot of the charge density nicely shows, the excess electrons are localized mainly in the space left by the removed oxygen atom. This electronic cloud has a strong influence on the migration of the vacancy, as we will see below. The role of the two excess electrons is also clearly visible in the calculated DOS (see Figure 4). In the case of V${}_{\mathrm{O}}$, a filled state within the gap appears (panel (a) in Figure 4) that corresponds to the two electrons “left behind” after the removal of a neutral O atom. In the case of V${}_{\mathrm{O}}$${}^{2+}$ instead (panel (b) in Figure 4), no such deep gap state is present, and the system’s band gap is close to the gap of the pristine c-ZrO${}_2$ crystal.

We used the Nudged Elastic Band (NEB) method [31] to study the migration of V${}_{\mathrm{O}}$ to its nearest neighbour oxygen site. Due to the symmetry of c-ZrO${}_2$, all oxygen neighbours of a given oxygen atom are equivalent, and it is therefore sufficient to consider one single pair of nearest-neighbour O-atoms.

Figure 5 shows the diffusion barriers for the case of V${}_{\mathrm{O}}$ (panel (a)) and V${}_{\mathrm{O}}$ (panel (b)). For the NEB calculations, we used two different sets of values for the Hubbard-U parameters. The values of U obtained by a linear response calculation lead to the green curves, while the alternative set of parameters leads to the red curves. As discussed in Section 3.1, no empirical information or fitting is used for the values leading to the green curves, while the red curves are a choice that can be found in the literature [1,30] and which has been fitted. Both possible choices of U parameters lead to the correct band gap in agreement with the experiment.

It is important to note that also in a system such as ZrO${}_2$, where the formal occupation of the 4d-states is zero, the choice of U has such an important influence on energies in the electronic ground state, as evidenced in Figure 5. As it is mostly the conduction band manifold that has a strong Zr-4d character (see Figure 1 and Figure 4), one can clearly expect that the Hubbard corrections will influence the band gap or optical properties. However, in conjunction with the Hubbard-U for the oxygen 2p-states, there is also a strong dependence on the energies along the vacancy migration path on the U-parameters.

In the absence of experimental information on the energy of the transition state, we have resorted to a hybrid calculation (using PBE0) to gauge the quality of the obtained NEB curves. In PBE0, the self-interaction error is partially corrected thanks to the incorporation of a fraction of Hartree–Fock exchange—and this without resorting to Hubbard-U corrections. In Figure 5, it can be seen that the U values from linear response lead to curves which are in closer agreement with the PBE0 curves than the empirically chosen values for U. Apart from the band gap, this is a second test of the linear response parameter set which provides a consistently good description of both band gaps and ground-state energies.

Similar to the scan of U values for the band gap in Figure 2, we have also performed a scan of the barrier for ${\mathrm{V}}_{\mathrm{O}}$ migration. It is shown in Figure 6. For this figure, we used the geometries of the initial state and of the transition state obtained by NEB with the U-parameters from linear response. The diffusion barrier shown in Figure 6 is the difference in total energy between these two geometries, calculated with varying values of U. The black line indicates the barrier height from a PBE0 calculation, which we might take as a reference. Comparing the black lines in Figure 2 and Figure 6, one can see that the calculated values (U(Zr-4d) = 2.85 eV and U(O-2p) = 9.16 eV) are indeed close to the intersection of the two lines.

Concerning the absolute value of the height of the diffusion barrier, it is interesting to see that V${}_{\mathrm{O}}$ has a very high barrier of 1.69 eV, which effectively prevents thermal activation and migration of an oxygen vacancy in c-ZrO${}_2$. In contrast to this, the charged vacancy V${}_{\mathrm{O}}$${}^{2+}$ has no such barrier: at the PBE0 level of theory, we find that the vacancy needs to overcome only 0.22 eV to migrate to a neighbouring site. In this case, the vacancy is thus very mobile, even at room temperature. We attribute this effect to the localized charge cloud depicted in Figure 3. When the vacancy is close to the position of an oxygen atom of the pristine lattice, there is a cavity available where the excess electrons can localize. Close to the transition state, an oxygen atom finds itself in the middle between two equilibrium O positions. No cavity exists in this case, and the excess electrons cannot localize in the same way. This impediment of migration does not exist for the case of V${}_{\mathrm{O}}$${}^{2+}$ where no corresponding cloud of electrons needs to move together with the vacancy.

The ease of migration of V${}_{\mathrm{O}}$${}^{2+}$ with respect to V${}_{\mathrm{O}}$ supports the picture of zirconium corrosion detailed in Ref. [2]: the formation of V${}_{\mathrm{O}}$ is associated with a V${}_{\mathrm{O}}$${}^{2+}$/2e${}^{-}$ separation, leading to an anodic property of the metal–oxide interface and similarly to a cathodic property at the interface of the oxide and the oxidizing surroundings [2].

4. Conclusions

We have applied the DFT+U methodology to the case of c-ZrO${}_2$. Due to the d${}^0$ oxidation state of Zr in this compound, it can be expected that the standard application of a Hubbard-U correction to zirconium 4d-states only influences properties, such as the band gap, which involve empty states. We show that if Hubbard corrections are employed to both Zr-4d and O-2p states, then this correction has a strong influence also on ground-state energies. This is particularly apparent in the case of oxygen vacancies, for which the diffusion barrier depends strongly on the chosen values for U.

There exists a wide range of possible values for U(Zr-4d) and U(O-2p) which lead to a band gap prediction in close agreement with the experiment. A linear response calculation of the two parameters leads to U(Zr-4d) = 2.85 eV and U(O-2p) = 9.16 eV, two values which are very different from the usual choice of a high Hubbard-U for zirconium and a relatively small correction for oxygen. In spite of this unusual result, we show that this calculated choice of parameters leads to excellent agreement of the band gap with the experiment and also to a vacancy migration barrier that closely matches the results from a hybrid DFT calculation.

It is found that the creation of an oxygen vacancy by removal of a neutral O atom from the lattice leads to two excess electrons which localize in the created cavity. This charge is an impediment for the thermally activated migration of the vacancy, which has a barrier of 1.69 eV. In contrast to this situation, a charged oxygen vacancy V${}_{\mathrm{O}}$${}^{2+}$, where the two electrons are removed, has a much lower barrier of migration. These findings support the picture of a redox reaction associated with the creation of vacancies and thus the corrosion process of zirconium to form ZrO${}_2$.

Conflicts of Interest

The author declares no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. (a) The density of states (DOS) of c-ZrO[Forumla omitted. See PDF.] calculated at the PBE level of theory; (b) DOS of the same system, calculated with U(Zr-4d) = 2.85 eV and U(O-2p) = 9.16 eV. See text for a discussion of these values which have been calculated by linear response. The Fermi energy is set to zero in both panels.

Figure 2. Calculated band gap of c-ZrO[Forumla omitted. See PDF.] as a function of the Hubbard-U parameters. The black line indicates the experimental band gap of 5.7 eV.

Figure 3. Localization of the excess electrons upon formation of an oxygen vacancy. Zr atoms are shown in grey, oxygen atoms in red. The isosurface of the excess electron density is shown in light blue and is localized mainly around the oxygen vacancy in the middle of the graph.

Figure 4. DOS for c-ZrO[Forumla omitted. See PDF.] with an oxygen vacancy: (a) for the case of V[Forumla omitted. See PDF.]; (b) for the case of V[Forumla omitted. See PDF.][Forumla omitted. See PDF.]. The Fermi level has been set to 0 eV.

Figure 5. Nudged Elastic Band (NEB) calculation of the transition path for a vacancy migration: (a) for the case of V[Forumla omitted. See PDF.]; (b) for the case of V[Forumla omitted. See PDF.][Forumla omitted. See PDF.].

Figure 6. Calculated V[Forumla omitted. See PDF.] diffusion barrier as a function of the Hubbard-U parameters. The black line indicates the barrier height from a hybrid DFT calculation (using PBE0).

Table 1

Used parameters for all calculations.

ecutwfc	ecutrho	n${}_{\mathit{k}}$	Lattice Parameter
40 Ry	400 Ry	4	5.09 Å

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cryst13040574/s1, Convergence tests for ZrO${}_2$

References

1. Khattab, E.S.R.; Abd El Rehim, S.S.; Hassan, W.M.I.; El-Shazly, T.S. Band Structure Engineering and Optical Properties of Pristine and Doped Monoclinic Zirconia (m-ZrO₂): Density Functional Theory Theoretical Prospective. ACS Omega; 2021; 6, pp. 30061-30068. [DOI: https://dx.doi.org/10.1021/acsomega.1c04756] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34778677]

2. Lindgren, M.; Panas, I. Oxygen Vacancy Formation, Mobility, and Hydrogen Pick-up during Oxidation of Zirconium by Water. Oxid. Met.; 2017; 87, pp. 355-365. [DOI: https://dx.doi.org/10.1007/s11085-016-9695-z]

3. Zhu, H.; Li, J.; Chen, K.; Yi, X.; Cheng, S.; Gan, F. Nature of charge transport and p-electron ferromagnetism in nitrogen-doped ZrO₂: An ab initio perspective. Sci. Rep.; 2015; 5, 8586. [DOI: https://dx.doi.org/10.1038/srep08586] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25717033]

4. Raza, M.; Cornil, D.; Cornil, J.; Lucas, S.; Snyders, R.; Konstantinidis, S. Oxygen vacancy stabilized zirconia (OVSZ); a joint experimental and theoretical study. Scr. Mater.; 2016; 124, pp. 26-29. [DOI: https://dx.doi.org/10.1016/j.scriptamat.2016.06.025]

5. Fabris, S. A stabilization mechanism of zirconia based on oxygen vacancies only. Acta Mater.; 2002; 50, pp. 5171-5178. [DOI: https://dx.doi.org/10.1016/S1359-6454(02)00385-3]

6. Oka, M.; Kamisaka, H.; Fukumura, T.; Hasegawa, T. Density functional theory-based ab initio molecular dynamics simulation of ionic conduction in N-/F-doped ZrO₂ under epitaxial strain. Comput. Mater. Sci.; 2018; 154, pp. 91-96. [DOI: https://dx.doi.org/10.1016/j.commatsci.2018.07.038]

7. Lee, Y.L.; Duan, Y.; Morgan, D.; Sorescu, D.; Abernathy, H.; Hackett, G. Density functional theory modeling of cation diffusion in bulk tetragonal zirconia. Ceram. Trans.; 2019; 266, pp. 95-110. [DOI: https://dx.doi.org/10.1002/9781119631460.ch10]

8. Schultze, T.K.; Arnold, J.P.; Grieshammer, S. Ab Initio Investigation of Migration Mechanisms in La Apatites. ACS Appl. Energy Mater.; 2019; 2, pp. 4708-4717. [DOI: https://dx.doi.org/10.1021/acsaem.9b00226]

9. Koga, H.; Hayashi, A.; Ato, Y.; Tada, K.; Hosokawa, S.; Tanaka, T.; Okumura, M. Facile NO-CO elimination over zirconia-coated Cu (110) surfaces: Further evidence from DFT + U calculations. Appl. Surf. Sci.; 2020; 508, 145252. [DOI: https://dx.doi.org/10.1016/j.apsusc.2020.145252]

10. Mueller, M.P.; Pingen, K.; Hardtdegen, A.; Aussen, S.; Kindsmueller, A.; Hoffmann-Eifert, S.; De Souza, R.A. Cation diffusion in polycrystalline thin films of monoclinic HfO₂ deposited by atomic layer deposition. APL Mater.; 2020; 8, 081104. [DOI: https://dx.doi.org/10.1063/5.0013965]

11. Lee, Y.L.; Duan, Y.; Sorescu, D.C.; Morgan, D.; Abernathy, H.; Kalapos, T.; Hackett, G. Density functional theory modeling of cation diffusion in tetragonal bulk ZrO₂: Effects of humidity and hydrogen defect complexes on cation transport. Phys. Rev. Res.; 2021; 3, 013121. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.013121]

12. Kumar, N.; Seriani, N.; Gebauer, R. DFT insights into electrocatalytic CO₂ reduction to methanol on α-Fe₂O₃ (0001) surfaces. Phys. Chem. Chem. Phys.; 2020; 22, pp. 10819-10827. [DOI: https://dx.doi.org/10.1039/C9CP06453B] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32400771]

13. Ulman, K.; Poli, E.; Seriani, N.; Piccinin, S.; Gebauer, R. Understanding the electrochemical double layer at the hematite/water interface: A first principles molecular dynamics study. J. Chem. Phys.; 2019; 150, 041707. [DOI: https://dx.doi.org/10.1063/1.5047930] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30709242]

14. Ahamed, I.; Ulman, K.; Seriani, N.; Gebauer, R.; Kashyap, A. Magnetoelectric e-Fe₂O₃: DFT study of a potential candidate for electrode material in photoelectrochemical cells. J. Chem. Phys.; 2018; 148, 214707. [DOI: https://dx.doi.org/10.1063/1.5025779] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29884044]

15. Nguyen, M.T.; Gebauer, R. Graphene Supported on Hematite Surfaces: A Density Functional Study. J. Phys. Chem. C; 2014; 118, pp. 8455-8461. [DOI: https://dx.doi.org/10.1021/jp501171t]

16. Walker, B.G.; Hendy, S.C.; Gebauer, R.; Tilley, R.D. Application of Lanczos-based time-dependent density-functional theory approach to semiconductor nanoparticle quantum dots. Eur. Phys. J. B; 2008; 66, pp. 7-15. [DOI: https://dx.doi.org/10.1140/epjb/e2008-00388-1]

17. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G.L.; Cococcioni, M.; Dabo, I. et al. QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter; 2009; 21, 395502. [DOI: https://dx.doi.org/10.1088/0953-8984/21/39/395502]

18. Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Buongiorno Nardelli, M.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys. Condens. Matter; 2017; 29, 465901. [DOI: https://dx.doi.org/10.1088/1361-648X/aa8f79]

19. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett.; 1996; 77, pp. 3865-3868. [DOI: https://dx.doi.org/10.1103/PhysRevLett.77.3865]

20. Dal Corso, A. Pseudopotentials periodic table: From H to Pu. Comput. Mater. Sci.; 2014; 95, pp. 337-350. [DOI: https://dx.doi.org/10.1016/j.commatsci.2014.07.043]

21. Perdew, J.P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys.; 1996; 105, pp. 9982-9985. [DOI: https://dx.doi.org/10.1063/1.472933]

22. Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys.; 1999; 110, pp. 6158-6170. [DOI: https://dx.doi.org/10.1063/1.478522]

23. Mulliken, R.S. Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I. J. Chem. Phys.; 1955; 23, pp. 1833-1840. [DOI: https://dx.doi.org/10.1063/1.1740588]

24. Bader, R.F.W. Atoms in Molecules—A Quantum Theory; Oxford University Press: Oxford, UK, 1990.

25. Sit, P.H.L.; Car, R.; Cohen, M.H.; Selloni, A. Simple, unambiguous theoretical approach to oxidation state determination via first-principles calculations. Inorg. Chem.; 2011; 50, pp. 10259-10267. [DOI: https://dx.doi.org/10.1021/ic2013107]

26. Krukau, A.V.; Vydrov, O.A.; Izmaylov, A.F.; Scuseria, G.E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys.; 2006; 125, 224106. [DOI: https://dx.doi.org/10.1063/1.2404663] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17176133]

27. Yang, Y.L.; Fan, X.L.; Liu, C.; Ran, R.X. First principles study of structural and electronic properties of cubic phase of ZrO₂ and HfO₂. Phys. B Condens. Matter; 2014; 434, pp. 7-13. [DOI: https://dx.doi.org/10.1016/j.physb.2013.10.037]

28. Cococcioni, M.; de Gironcoli, S. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B; 2005; 71, 035105. [DOI: https://dx.doi.org/10.1103/PhysRevB.71.035105]

29. Timrov, I.; Marzari, N.; Cococcioni, M. Hubbard parameters from density-functional perturbation theory. Phys. Rev. B; 2018; 98, 085127. [DOI: https://dx.doi.org/10.1103/PhysRevB.98.085127]

30. Li, J.; Meng, S.; Niu, J.; Lu, H. Electronic structures and optical properties of monoclinic ZrO₂ studied by first-principles local density approximation + U approach. J. Adv. Ceram.; 2017; 6, pp. 43-49. [DOI: https://dx.doi.org/10.1007/s40145-016-0216-y]

31. Henkelman, G.; Uberuaga, B.P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys.; 2000; 113, pp. 9901-9904. [DOI: https://dx.doi.org/10.1063/1.1329672]