manuscript.org

Investigating the Energetic Ordering of Stable and Metastable \ce{TiO_2} Polymorphs Using DFT+$U$ and Hybrid Functionals

\maketitle

Introduction

Transition metal dioxides characterized as BO₂ are used in applications that involve supportive substrates in CO oxidation processes cite:Francisco2001,Bokhimi2007,Hamad2009, photocatalysis cite:Nakata2012,Indrakanti2008, electrochemical water splitting cite:Osterloh2008,Man2011, electrode coating and doping cite:Hu2006,Subasri2006,Velichenko2009, gas sensing for controlling engine fuel consumption cite:Diebold2003,Lukacevic2012, and cleaning automative exhuast cite:Reddy2005,Hamze2012, and electronic devices cite:Chiang1994,Hamad2009. In several of the applications that are listed above and others cite:Lundin1998,Dubrovinsky2001,Hugosson2002, the accurate characterization of properties, such as band gap cite:Mattesini2004_115101,Kuo2005 and the temperature and pressure dependence of phase stability cite:Olsen1999,Withers2003, of oxide materials is critical.

With respect to most photochemical applications, Anatase is the most photocatalytically active TiO₂ polymorph that naturally forms under ambient conditions cite:Indrakanti2008,Hashimoto2005,Diebold2003, though its wide band gap (E_g = 3.0-3.2 eV) limits its solar energy conversion efficiency cite:Mattesini2004_115101,Fang2007,Wu2010. Recent theoretical studies have postulated that the TiO₂ band gap can be tuned by changing its polymorphic structure from that of Anatase (tetragonal) to a cubic structure such as Fluorite cite:Mattesini2004_115101,Kuo2005 or pyrite cite:Tsai2006,Ataca2012. The thermodynamic stability of these cubic phases under particular reaction conditions has been considered experimentally cite:Haines1993,Lagarec1995,Mattesini2004_212101 and theoretically cite:Swamy2014,Mei2014,Zhou2010, though none of these studies conclusively characterize their stability.

Despite the enhanced recombination effects potentially resulting from doping, cite:Tsai2006 mixed-metal oxides composed of Fluorite-structured CeO₂ and Anatase-structured TiO₂ have also been employed to tune the band gap, additionally improving catalyst properties such as effective surface area and thermostability. cite:Francisco2001,Reddy2005,Fang2007,Bokhimi2007 Varying CeO₂ composition can tune the band gap to maximize visible light absorption in TiO₂-CeO₂ catalysts cite:Pavasupree2004 and substantially improve their photochemical activity without known negative consequences cite:Munoz-Batista2013. With particular consideration to phase stability, materials property characterization can only be fully addressed via a comprehensive methodology applicable to a broad materials set. Applications such as CO oxidation employing BO₂ supports cite:Hamad2009 and epitaxial stabilization cite:Gorbenko2002,Mehta2014 employ metal oxide candidates featuring a wide range of transition metal cations (B = Ti, V, Cr, Mn, Ru, Ir). cite:Hamad2009,Mehta2014

Density functional theory (DFT) has been previous applied as a computationally efficient cite:Chevrier2010,Akhade2012 component of a methodology capable of predictively characterizing broad sets of metal oxide properties ranging from relative energetic stability cite:Jain2011,Hautier2012,Curnan2014 to electronic structure cite:Akhade2011,Xu2014. Furthermore, this methodology should also be able to minimally predict the relative ordering of values of properties cite:Curnan2014, such as the formation energy ordering of various TiO₂ polymorphs cite:Ma2009,Mehta2014, that are close in value consistently cite:Ma2009,Mehta2014,Curnan2014. However, previous use of DFT to determine the relative energetics of TiO₂ cite:Labat2007,Ma2009,Mehta2014 structures has yielded functional-dependent results. Across the development of different basis sets and orbital overlap criteria, use of the LDA functional consistently predicts the relative energetic favorability of Rutile relative to Anatase, serving as the only functional capable of consistently demonstrating bulk Rutile to be more energetically favorable than bulk Anatase over multiple basis sets and orbital overlap criteria. Under a single basis set and orbital overlap criteria combination, the use of pure HF exchange can also achieve the relative energetic favorability of Rutile, while all other combinations of input produce HF functional results with lower relative energetic differences favoring Anatase than the corresponding energetic differences produced by non-LDA functionals. In all other tested combinations of basis set and orbital overlap criteria, spanning consideration of HF, LDA, PBE, PBE0, and B3LYP functionals, Anatase is predicted to be more energetically favorable than Rutile cite:Labat2007.

Experimentally and theoretically, the relative phase stability of TiO₂ polymorphic structures has been determined to be strongly affected by particle size cite:Hu2006,Koparde2008,Lukacevic2012, particle shape cite:Barnard2005,Tsai2006, pressure and temperature conditions during synthesis cite:Koparde2008,Li2012,Swamy2003, and the solvent surrounding that phase cite:Koparde2008,Barnard2005. Studies completed using the ReaxFF method cite:Duin2001 have begun to characterize the relationships between phase stability and these factors cite:Raju2014,Kim2012 instead of DFT due to computational cost cite:Raju2013,Kim2013_5655. Surface structural configurations cite:Kim2013_7838, atomic bond lengths and angles cite:Kim2013_5655,Kim2012, and bulk structure equations of state (EOS) cite:Raju2013,Raju2014 calculated with DFT are used to parameterize ReaxFF calculations, thus accurate electronic structure calculations are still needed to reliably characterize TiO₂ phase stability.

Though the reasons for the inconsistent energetic ordering results for TiO₂ across functionals are still contested, cite:Conesa2010 a probable source is the non-systematic error incurred by unphysical electron-electron interactions cite:Wang2006,Liechtenstein1995, which is commonly observed in 3$d$ transition metal oxide DFT calculations completed using standard local density (LDA) and generalized gradient (GGA) approximation functionals cite:Zhou2004,Jain2011. The electron-electron interaction error in these calculations can be accounted for using hybrid functionals such as PBE0 and HSE06, in which the delocalization of 3$d$ Ti electron charge density typical in standard functionals is countered by the matching localization effects induced by employing exact Hartree-Fock (HF) exchange cite:Cohen2008.

Considering that the relative energetic ordering of Rutile and Anatase TiO₂ polymorphs can change when switching between using standard functionals cite:Wu2010,Mei2014 and exact HF exchange cite:Fahmi1993, hybrid functionals can apparently resolve the correct relative energetics of TiO₂ polymorphs albeit with large computational expense cite:Chevrier2010. Though physical electron delocalization behavior was achieved using 20% HF exchange with Lee, Yang, and Parr correlation in these TiO₂ polymorphs cite:Finazzi2008, physically reasonable behavior was not achieved in oxides such as NiO using this amount of HF exchange cite:Finazzi2008,Morgan2010. Thus, the use of a single percentage of HF exchange cannot be generalized over many materials and cannot be applied to a predictive methodology.

Electron-electron interaction error can also be accounted for with the Hubbard $U$ model, which observes Rutile-Anatase energetic ordering changes over specific ranges of $U$ applied to the 3$d$ orbitals of Ti and gives physical delocalization behavior over the range U_Ti,3d = 3.0-4.0 eV cite:Finazzi2008. Furthermore, $U$ parameters between 3.0-3.5 eV have already been calculated for several TiO₂ polymorphs using the first-principles method of linear response theory cite:Cococcioni2005,Mattioli2008,Mattioli2010, illustrating the potential for employing this method to consistently and accurately predict phase stability across broad sets of BO₂ polymorphs cite:Cococcioni2005,Hamad2009.

This study will focus on determining the extent to which a consistent relative energetic ordering can be achieved over calculations that vary in functional and pseudopotential selection. Under typical ambient conditions, three naturally occurring TiO₂ polymorphs are expected to occur, namely Rutile, Anatase, and Brookite cite:Diebold2003,Koparde2008,Hamad2009,Wu2010,Lukacevic2012. Rutile has been experimentally shown to be the only naturally occurring stable bulk TiO₂ polymorph, cite:Swamy2002,Hu2006,Wu2010,Lukacevic2012 thus this study will investigate whether DFT can consistently confirm that the Rutile → Anatase and Rutile → Brookite formation energies are positive. Additionally, this study will confirm whether DFT can predict a consistent ordering of the metastable TiO₂ polymorphs Anatase and Brookite, namely E_Rutile \textless E_Anatase \textless E_Brookite or E_Rutile \textless E_Brookite \textless E_Anatase. Though experimental studies have shown that factors such as TiO₂ particle size and synthesis temperature can impact the relative stability of TiO₂ polymorphs that are less thermodynamically accessible than Brookite and Anatase, cite:Swamy2003,Li2012 extensive research on the bulk phase stability of TiO₂ with the application of hydrostatic pressure have shown that the Columbite TiO₂ polymorph is the next most accessible phase during compression cite:Ohsaka1979,Mammone1980,Arashi1992,Lagarec1995,Olsen1999,Sekiya2001, decompression cite:Mattesini2004_212101,Wu2010, or both cite:Arlt2000,Swamy2002 cycles. Therefore, this study will investigate whether the Rutile → Anatase and Rutile → Brookite formation energies are consistently less than that of Rutile → Columbite. Ranges of Hubbard $U$ values cite:Liechtenstein1995,Zhou2004,Wang2006,Jain2011 and fractions of exact exchange cite:Cohen2008,Finazzi2008 will be applied to DFT calculations to assess the intervals within which consistent relative energetic orderings can be achieved cite:Finazzi2008,Dompablo2011. Linear response theory will be applied to determine whether first-principles methods can calculate $U$ values that reliably fall within those intervals \textit{a priori} cite:Cococcioni2005,Kulik2006.

Methodology

Four polymorphic TiO₂ structures were investigated in this study, namely Rutile (space group $P4₂/mnm$), Anatase (space group $I4₁/amd$), Brookite (space group $Pbca$), and Columbite (space group $Pbcn$). This study will not focus on the synthetically formed Ramsdellite (TiO₂(R)) and Hollandite (TiO₂(H)) polymorphs that can be thermodynamically accessible under ambient conditions cite:Swamy2014. Initial input for the atomic positions and lattice parameters of these structures was provided from experimental data cite:Muscat2002. In all calculations involving relative energetic ordering, unit cell representations of the Rutile (6 atom), Anatase (6 atom), Columbite (12 atom), and Brookite (24 atom) polymorphs were evaluated. The Birch-Murnaghan equation of state (EOS) was employed to calculate the ground state energies of these systems cite:Murnaghan1944.

With respect to the Vienna Ab-initio Simulation Package (VASP) calculations cite:Kresse1996,Kresse1999 completed in this study, all pseudopotentials used were generated by the projector augmented wave (PAW) cite:Blochl1994 method. Relative energetic ordering was tested using the local density approximation (LDA) cite:Perdew1981 exchange-correlation functional and four generalized gradient approximation (GGA) functionals, featuring the Perdew-Burke-Ernzerhof (PBE) cite:Perdew1996_PRL, Perdew-Burke-Ernzerhof for solids (PBEsol) cite:Perdew2008, Armiento-Mattsson (AM05) cite:Armiento2005, and Perdew-Wang (PW91) cite:Perdew1992 parameterizations. VASP version 5.2.12 was used in all calculations except for those involving the PW91 functionals, which used VASP version 5.3.5 cite:Kresse1996,Kresse1999. Using the PBE functional, five combinations of Ti and O pseudopotentials were investigated in energetic ordering calculations. They are labelled Ti and O (the default pseudopotentials), Ti_pv (treating $p$ electrons as valence electrons) and O, Ti_sv (treating $s$ electrons as valence) and O, Ti_pv and O_s (soft pseudopotential for oxygen), and Ti_sv and O_s in the VASP software package. Different available PAW pseudopotentials for O, including the harder O pseudopotential (O_h), were not investigated due to the higher required energy cutoffs required by their use cite:Kresse1996,Kresse1999.

The rotationally invariant Dudarev implementation of the Hubbard \(U\) model was used to account for electron-electron interaction error in spin-polarized, paramagnetic (PM) TiO₂ polymorph calculations cite:Dudarev1998,Finazzi2008. In this implementation, the on-site Coulombic (\(U\)) and Exchange (\(J\)) terms are combined into a single effective \(U\) parameter ($U_eff$) to account for errors in exchange correlation on Ti $3d$ orbitals cite:Dudarev1998. Calculations completed using this procedure in VASP employed a 600 eV plane-wave energy, an $8 × 8 × 8$ Monkhorst-Pack cite:Monkhorst1976 $k$-point sampling mesh set with respect to each TiO₂ unit cell (containing 6, 12, or 24 atoms), and a 0.01 eV/Å force tolerance on each atom.

Corresponding calculations employing the PBE0 cite:Perdew1996_JChemPhys and HSE06 cite:Heyd2003,Heyd2006 hybrid functionals (VASP version 5.2.12) and PAW pseudopotentials were completed for the Rutile and Anatase TiO₂ polymorphs. Calculations completed using this procedure employed a 550 eV plane-wave energy, a $6 × 6 × 6$ Gamma point centered cite:Jain2011 $k$-point sampling mesh set with respect to each tested TiO₂ unit cell (containing 6 atoms), and a 0.02 eV/Å force tolerance on each atom. Calculation of paramagnetic, first-principles derived Hubbard \(U\) values for pertinent TiO₂ polymorph systems is accomplished via implementation of linear response theory in VASP (version 5.3.5) cite:Cococcioni2002,Zhou2004,Cococcioni2005. These linear response calculations demonstrated the effects of varying Ti pseudopotential selection (Ti, Ti_pv, or Ti_sv) and GGA functional parameterization (PBE or PBEsol), employing a 600 eV plane-wave energy cutoff and an $8 × 8 × 8$ Monkhorst-Pack cite:Monkhorst1976 $k$-point sampling mesh set with respect to each TiO₂ unit cell. $2 × 2 × 2$ supercells were used to perform linear response calculations on Rutile (48 atoms), Anatase (48 atoms), and Columbite (96 atoms) systems, while a $2 × 2 × 1$ supercell was applied in corresponding calculations on Brookite (96 atoms). These selections were made to balance the computational expense of performing linear response calculations on larger systems with the possible “delocalized background” contributions to calculated $U$ values resulting from the necessary use of finite supercells cite:Cococcioni2005.

Further information on the procedures used to complete calculations involving $U_eff$, hybrid functionals, and linear response $U$ values, as well as information on structural relaxation approaches used in particular EOS calculations, can be found in the Supporting Information document for this study.

Results and Discussion

The results of evaluating the energetic ordering of Rutile, Anatase, Columbite, and Brookite TiO$₂$ polymorphs as a function of the Hubbard $U_3d$ parameter on Ti $3d$ orbitals are calculated and presented in accordance with several conventions. Firstly, formation energies used to evaluate energetic ordering, which are calculated as the DFT total energy differences between Rutile and any other polymorph, are normalized with respect to the TiO₂ formula unit. Secondly, the ranges of $U$ values over which the relative energetic ordering of stable and metastable bulk TiO₂ polymorphs is consistent with experimental results cite:Swamy2002,Hu2006 are evaluated over $U_3d$ intervals of either 0.0-6.0 eV or 0.0-9.0 eV in 1.0 eV increments. These intervals are selected due to considerations of energetic ordering switching cite:Dompablo2011, electron delocalization cite:Finazzi2008,Morgan2010, and $U$ parameter calculation cite:Mattioli2008,Mattioli2010 completed in previous literature reports. Thirdly, despite the use of first-principles and experimentally fitted $U$ parameters to the 2$p$ orbitals of O in previous work cite:Mattioli2010, $U$ parameters were not applied to these orbitals, as the Hubbard \(U\) model can only physically characterize $d$-$d$ and $f$-$f$ orbital interactions effectively cite:Tsuneda2014. Experimental fitting of \(U\) values is also not suitable for accurate prediction of energetic properties in different materials sets cite:Hamad2009,Mehta2014 or integration into different first-principles calculations requiring physically realistic input, cite:Raju2013,Kim2013_5655 as \(U\) values themselves represent intrinsic responses to orbital occupation perturbation cite:Cococcioni2005,Kulik2010. Lastly, in order to illustratively compare the relative energetic differences of any set of polymorphs characterized by any sets of functionals or pseudopotentials studied, a graphical visualization procedure is described in the Supporting Information document of this study. However, demonstration of the impact of various factors on energetic ordering will be illustrated through several figures within this study. These figures will directly highlight comparisons between GGA (PBE) and LDA functional calculations (Figure ref:fig-PBELDAPBEs), PBE functional calculations incorporating Ti $p$-shell valence inclusive pseudopotentials (Figure ref:fig-pvpvs), PBE functional calculations incorporating Ti $s$-shell valence inclusive pseudopotentials (Figure ref:fig-svsvs), and GGA functional calculations of varied (PBEsol, PW91, AM05) parameterization (Figure ref:fig-PSPW91AM05), respectively.

Effects of pseudopotentials

As is shown in Figure ref:fig-PBELDAPBEs, the use of softer O psuedopotentials (O_s) on calculations employing the PBE functional has little effect on energetic ordering predictions, which successfully find Rutile to be the most energetically favorable polymorph above approximately $U_3d$ = 2.8 eV, Brookite to be more stable than Rutile above approximately $U_3d$ = 2.0 eV, Anatase to be less stable than Rutile and consistently more stable than Brookite between approximately $U_3d$ = 2.8-4.0 eV, and Columbite to be the least stable polymorph between approximately $U_3d$ = 0.3-4.3 eV. Therefore, in an approximate range of $U_3d$ = 2.8-4.3 eV, all experimentally expected results for energetic ordering criteria are met. In contrast, use of the LDA functional shows Rutile to be the most energetically favorable polymorph above approximately $U_3d$ = 1.5 eV, Brookite to always be more stable than Anatase though less stable than Rutile above approximately $U_3d$ = 0.9 eV, and Anatase to be less stable than Rutile and Brookite above approximately $U_3d$ = 0.4 eV. However, Columbite is never predicted to be the least energetically stable of the four phases at any tested $U_3d$ value when using the LDA functional. Therefore, LDA is not capable of entirely resolving experimentally consistent relative energetic stability upon applying any single value of $U_3d$ to all tested TiO₂ polymorphs.

Comparison of Brookite and Columbite results achieved with different functionals reveals that – with respect to the PBE functional phase lines of varying $U_3d$ parameter magnitude in Figures ref:fig-PBELDAPBEs, ref:fig-pvpvs, ref:fig-svsvs, and ref:fig-PSPW91AM05 – the relative formation energy of LDA functional resolved Brookite increases and Columbite decreases. Nevertheless, the largest change in LDA functional results is seen in the large upward shift of the relative formation energetics of Anatase. Figures ref:fig-pvpvs and ref:fig-svsvs, which incorporate $p$-shell and $s$-shell valence inclusive Ti pseudopotentials and O pseudopotentials of varying softness, demonstrate similar energetic behavior. In both Figures, which feature four pseudopotential combinations, Columbite is the least stable polymorph above approximately $U_3d$ = 1.0 eV and Brookite is the next least stable between approximately $U_3d$ = 3.1-7.0 eV. However, the relative energetics of the $p$-shell valence inclusive Ti pseudopotential calculations show that Anatase is unphysically more stable than Rutile at $U_3d$ values less than 4.7 eV, whereas the analogous change in energetic favorability occurs at approximately $U_3d$ = 5.8 eV in calculations featuring $s$-shell valence inclusive Ti pseudopotentials.

With respect to the value of $U_3d$ at which Rutile-Anatase energetic ordering changes, these two sets of calculations are more consistent with previously derived results cite:Dompablo2011. Overall, the results analyzed thus far indicate that, for all VASP calculations incorporating the PBE parameterization of the GGA functional in accompaniment with varied types of Ti and O pseudopotentials, there exists a pseudopotential-dependent $U_3d$ range at which the polymorph energetic ordering E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite is maintained. In Figures ref:fig-pvpvs and ref:fig-svsvs, an energetic ordering consistent with experimental results, namely E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite, is respectively maintained within the intervals 4.7-7.0 and 5.8-8.2 eV. These $U_3d$ ranges are largely maintained regardless of whether soft (O_s) or standard (O) oxygen pseudopotentials are used.

Effects of exchange-correlation functionals

As shown in Figure ref:fig-PSPW91AM05, the relative energetics of several sets of calculations involving different types of GGA functionals consistently reveal ranges of $U_3d$ over which the E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite energetic ordering is preserved. In this Figure, results for the TiO₂ Brookite polymorph featuring use of the AM05 functional were omitted due to convergence issues. In the case of PW91 and AM05 functionals, Anatase becomes more stable than Rutile at around $U_3d$ = 2.7 or 2.8 eV (respectively) and remains less favorable than Columbite until approximately 4.3 or 4.0 eV (respectively). In the case of the PBEsol functional, an upward shift in the relative energy of Anatase (w.r.t. that of Rutile) causes Anatase to become more favorable than Rutile at around $U_3d$ = 1.6 eV and less favorable than Columbite at approximately $U_3d$ = 2.1 eV. While PW91 Brookite becomes more favorable than Rutile above approximately $U_3d$ = 2.0 eV and stays less favorable than Anatase below approximately $U_3d$ = 4.0 eV, PBEsol Brookite becomes more favorable than Rutile at approximately $U_3d$ = 1.5 eV though remains less favorable than Columbite until $U_3d$ = 2.4 eV largely due to the increased formation energy values of Anatase.

Similarly, PW91 and AM05 Columbite are the least energetically favorable polymorphs within their functionals in the approximate range of $U_3d$ = 0.4-4.1 eV or 0.4-4.3 eV (respectively), while PBEsol Columbite is least favorable within the narrower $U_3d$ = 0.8-2.1 eV range. Thus, all VASP calculations incorporating GGA functionals preserve the E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite energetic ordering within ranges of $U_3d$ values on Ti, the breadth of which is most likely and prominently impeded by underestimation of the stability of Anatase relative to Rutile. The data sets featuring the PBEsol GGA (Figure ref:fig-PSPW91AM05) and LDA (Figure ref:fig-PBELDAPBEs) functionals most prominently illustrate this characteristic. Over all values of $U_3d$, these Figures also consistently illustrate a slighter increase in the relative formation energy of Brookite and a comparable decrease in the relative formation energy of Columbite. However, the magnitudes of these shifts in formation energy trends are only capable of impacting the size of the $U_3d$ range in which the E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite energetic ordering exists, instead of entirely precluding the existence of such a $U_3d$ range.

As was observed in previous research for the LDA functional cite:Labat2007 and shown above for LDA and PBEsol functionals, a consistent shift of the Rutile-Anatase relative energetic trend occurs for all values of \textit{U} tested. Previous research has indicated that functional selection has demonstrated little effect on calculated TiO₂ Rutile and Anatase band gaps cite:Dompablo2011,Hu2011 and that cell volumes for matching systems have been shown to be consistently lower when applying PBEsol and LDA functionals than PBE functionals cite:Mehta2014. Therefore, the constant energetic shifts of LDA and PBEsol Rutile-Anatase energetic trends appear to be linked to changes in equilibrium volume rather than electronic structure. In contrast, electronic structure features such as band gap have been proposed to be more directly linked to \textit{U} parameters placed on the 2\textit{p} orbitals of O anions rather than the 3\textit{d} orbitals of Ti cations in TiO₂, as the application of high magnitude \textit{U}_3d,Ti values needed to reproduce experimental band gaps does not reproduce physical vacancy defect states cite:Angelis2014,Agapito2015. With respect to previous first-principles research concerning first-principles Rutile cite:Dompablo2011,Han2011 and Anatase cite:Islam2011,Han2013 electronic structure calculations, the changing of Rutile-Anatase energetic ordering across these energetic trends appears to be most strongly related to the narrowing and upward contraction of Ti 3\textit{d} conduction band densities of state with the increase of \textit{U}_3d,Ti. A comparison between bulk Rutile cite:Han2011 and Anatase cite:Han2013 band structures reveals that, with respect to their corresponding GGA electronic structures, the contribution of comparable values of \textit{U}_3d,Ti to both polymorphs disproportionately changes the differences between their conduction band (CB) minima and valence band (VB) maxima. Thus, this disproportionality between changes in the CB-VB differences of Rutile and Anatase with the incrementation of \textit{U}, which features greater increases in the CB-VB difference of Rutile with increasing \textit{U} resulting from a more pronounced contraction of the Ti 3\textit{d} states occupying its CB minima, serves to possibly indicate a physical link between changes in Rutile-Anatase relative energetic ordering and the corresponding electronic structures of Rutile and Anatase cite:Dompablo2011.

Linear response U

In order to evaluate the extent to which relative TiO₂ polymorph energetics can be predicted, the linear response method cite:Cococcioni2005,Kulik2006 is employed to calculate the first-principles values of the Hubbard $U$ applied to all Ti cations in each system. In past research, a significant range of first-principles calculated $U$ values have been resolved for Rutile TiO₂ polymorph systems, namely those studied to evaluate different properties of interest. For example, a $U_3d,Ti$ value of 0.15 eV was achieved and applied to calculations improving TiO₂ band gap estimation relative to experimentally calculated values cite:Agapito2015, while a corresponding value of 4.95 eV was achieved and implemented in calculations comparing large sets of surface adsorption energetics to matching oxygen evolution reaction (OER) activity trends cite:Xu2015. Furthermore, application of the self-consistent extension cite:Kulik2006 of the linear response approach to comparatively model intrinsic defect behavior in TiO₂ Rutile and Anatase achieved respective first-principles resolved $U$ values of 3.23 and 3.25 eV for each system cite:Mattioli2010.

As shown in Table ref:table-Linear_Response_U, the $U$ values achieved in this study are not only comparable to several previously calculated values, but also illustrate possible reasons for differences in those previously calculated $U$ values. When applying the PBE functional parameterization and the $p$-valence inclusive Ti pseudopotential to the linear response calculation of $U_3d,Ti$ for Rutile in VASP, a $U$ value of 4.773 ± 0.267 eV is resolved. This result, which is consistent with that achieved by Xu et al. (4.95 eV) within uncertainty cite:Xu2015, was obtained in Quantum Espresso (QE) cite:Giannozzi2009 while using the same functional as that applied in Xu et al. and a semi-core $p$-valence state inclusive, Ultrasoft, GBRV Ti pseudopotential featuring a non-linear core correction cite:Garrity2014. Similarly, applying the PBE functional parameterization and a standard pseudopotential in VASP calculations yields a result of 3.102 ± 0.137 eV, which is equivalent to the QE resolved result shown in Mattioli et al. (3.23 eV) within uncertainty cite:Mattioli2010. This result from Mattioli et al. was achieved while also using the PBE functional and a highly similar Ultrasoft pseudopotential, which differs from its GRBV analogue predominately in its lack of non-linear core correction cite:Giannozzi2009,Kulik2008. Despite the difference in the software package used to calculate $U$ value results in this paper (VASP) and past research (QE), strong similarity between the effects of implementing $p$-valence inclusive Ti pseudopotentials in VASP and including a non-linear core correction in a semi-core valence state inclusive Ti pseudopotential in QE is evident. Thus, when considering the effects of pseudopotential selection on TiO₂ polymorph linear response $U$ calculation independent of other factors, the use of pseudopotentials that include core states in their valence configurations apparently and systematically increases the values of calculated $U$ results.

The extent to which this increase in calculated $U$ values affects all considered polymorphs is additionally illustrated in Table ref:table-Linear_Response_U, as Rutile $U$ values largely uniformly increase from 3.102 eV to 6.030 or 4.773 eV in the cases of $s$-valence and $p$-valence inclusive pseudopotentials, respectively. Anatase $U$ values also largely uniformly increase from 2.929 eV to either 5.321 or 4.295 eV for $p$-valence and $s$-valence inclusive pseudopotentials, respectively. Consideration of the Rutile-Anatase relative formation energy $U$ range results derived in past research cite:Dompablo2011 or formerly in this study, in conjunction with the established link between pseudopotential selection and $U$ value, reveals that both $U$ values and their associated formation energies increase proportionately with changes in pseudopotential selection. In the case of changing from standard to $p$-valence or $s$-valence inclusive Ti pseudopotentials in Rutile-based calculations, changes in $U$ values and experimentally consistent $U$ ranges are equivalent and energetic ordering predictions are conserved within uncertainty as pseudopotentials are changed. However, this does not occur in the case of $p$-valence and $s$-valence inclusive pseudopotentials for Anatase-based calculations. In this case, $p$-valence and $s$-valence Anatase-based calculations would require $U$ values and associated uncertainties that respectively spanned ranges encompassing 4.7 and 5.8 eV, though Table ref:table-Linear_Response_U indicates that this is only true for Rutile-based calculations. A possible explanation for this inequivalence of $U$ value and $U$ range changes when assessing the Anatase polymorph is proposed in the Supporting Information document.

In contrast to changes in pseudopotential selection, changing the functional from PBE to PBEsol while maintaining standard pseudopotential selection does not cause a proportional, joint shift of calculated $U$ values and experimentally consistent $U$ ranges. This results from the non-systematic upward shift observed by the Rutile-Anatase formation energy curve, as is shown in Figure ref:fig-PSPW91AM05. Though both Rutile and Anatase calculated $U$ values decrease significantly when considering PBEsol functional results (w.r.t. PBE results) within uncertainty, neither polymorph $U$ value spans the narrow 1.5-2.1 eV $U$ range within which experimental energetic ordering consistency is achieved. Results employing the PBE functional and standard pseudopotentials are available for all four studied polymorphs. Within uncertainty, each of their $U$ values is explicitly above 2.8 eV and below 4.3 eV. Thus, in accordance with the $U$ ranges visualized in Figure ref:fig-PBELDAPBEs, these calculated $U$ results indicate that the energetic ordering of all four polymorph systems is necessarily experimentally consistent via the relationship E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite. Additionally, within uncertainty, all four polymorphs share a common set of $U$ values through which they can be directly energetically compared.

Polymorph	Functional	Pseudopotentials	$U$ (eV)
Rutile	PBE	Ti, O	3.102 ± 0.137
Anatase	PBE	Ti, O	2.929 ± 0.133
Columbite	PBE	Ti, O	2.983 ± 0.137
Brookite	PBE	Ti, O	2.934 ± 0.128
Rutile	PBEsol	Ti, O	2.727 ± 0.120
Anatase	PBEsol	Ti, O	2.558 ± 0.117
Rutile	PBE	Ti_pv, O	4.773 ± 0.267
Anatase	PBE	Ti_pv, O	4.295 ± 0.243
Rutile	PBE	Ti_sv, O	6.030 ± 0.392
Anatase	PBE	Ti_sv, O	5.321 ± 0.370

Impact of hybrid functionals and fraction of exact exchange

The use of different functionals has been shown to impact predictions of possible energetic ordering in stable, metastable, and unstable TiO₂ polymorphs, particularly affecting the Rutile-Anatase relative formation energy trend while varying $U_3d$. In past research, varying the $U_3d$ on Ti $3d$ orbitals as a tunable parameter over ranges comparable to 2.8-4.3 eV, namely the range shown capable of achieving experimentally consistent energetic ordering previously, has been demonstrated to affect the physicality of electron localization on those Ti $3d$ orbitals cite:Finazzi2008,Morgan2010. Hybrid functionals have been shown capable of affecting the relative energetic ordering of Rutile and Anatase to the possible extent of reversing their ordering cite:Fahmi1993. Proposals of the relationship between variation in the fraction of exact exchange contributed to a hybrid functional and the relative localization of valence ($3d$) electrons on Ti cations cite:Finazzi2008,Morgan2010 and in general applications cite:Cohen2008 have been made in past research. Considering that achieving physical electron localization behavior on Ti cations has been demonstrated possible by tuning either $U_3d$ in Hubbard $U$ calculations or the fraction of exact exchange in hybrid functional calculations, cite:Finazzi2008,Morgan2010 a Rutile-Anatase relative formation energy consistent with experimental expectations should be achievable within a range of either $U_3d$ values or fractions of exact exchange.

In accordance with the calculation of $U_3d$ intervals containing an experimentally consistent Rutile-Anatase energetic ordering, Figure ref:fig-PBE0HSE06 illustrates the fractions of exact exchange contributed to the PBE0 and HSE06 hybrid functionals in order to achieve E_Rutile \textless E_Anatase. When reviewed over several fractions of exact exchange (\textit{a} = 0.00, 0.25, 0.50, 0.75, 0.825, 0.875, 0.95, and 1.00) cite:Perdew1996_JChemPhys,Heyd2003, the PBE0 and HSE06 hybrid functionals both initially observe a monotonic upward trend when increasing the fraction of exact exchange, starting at a shared relative Rutile-Anatase formation energy of approximately -0.081 eV/TiO₂. Energetic ordering changes from E_Anatase \textless E_Rutile to E_Rutile \textless E_Anatase occur for the PBE0 hybrid functional at approximately \textit{a} = 0.72 and occur for the HSE06 hybrid functional at approximately \textit{a} = 0.78. As further detailed in the Supporting Information document, these monotonic increases in Rutile-Anatase relative formation energy appear to occur until the limit of HF exact exchange is reached.

For the HSE06 and PBE0 functionals, formation energy peaks are separated by approximately 0.015 eV/TiO₂, with PBE0 having a higher peak located between 0.040 and 0.045 eV/TiO₂. A more exact interpolation of the points at which the Rutile-Anatase relative formation energy peaks and reverts back to E_Anatase \textless E_Rutile is impaired by both the narrowness of the range of \textit{a} over which the monotonic decline in formation energy is observed and the low number of data points resolved within that range. Despite the achievement of an experimentally consistent reversal in the Rutile-Anatase relative energetic ordering at higher fractions of HF exact exchange, improvement of the exact exchange fraction in Lee, Yang, and Parr (LYP) parameterized hybrid functionals from 20% (B3LYP or Becke, three parameter, Lee-Yang-Parr) to 50% (H&HLYP or Half and Half, Lee-Yang-Parr) led to the calculation of physically unrealistic electronic structures in Rutile and Anatase TiO₂ cite:Finazzi2008,Morgan2010. Therefore, even though the improvement of exact exchange fraction in hybrid functionals can lead to experimentally consistent energetics in BO₂ systems, the physicality of the electronic structure yielding those structures cannot be guaranteed using solely hybrid functional energetic calculations.

Despite the evidence of consistent change in Ti $3d$ electron localization and delocalization behavior with the incrementation of $U_3d$ value or exact exchange fraction (\textit{a}) cite:Finazzi2008,Cohen2008, the Rutile-Anatase relative formation energy trend is clearly not monotonic over all values of the exact exchange fraction. Therefore, electron localization on Ti $3d$ orbitals cannot be the only factor impacting relative formation energetics in Rutile and Anatase polymorphs. The short-range, long-range (SR-LR) separation parameter ($ω$ = 0.2) of the HSE06 hybrid functional cite:Heyd2003, which distinguishes it from the related PBE0 functional, defines the long-range distance after which HF exchange is replaced with PBE exchange in these calculations. This replacement of exchange has clear effects on the relative Rutile-Anatase energetics when considering that the range in which HSE06 achieves E_Rutile \textless E_Anatase is significantly narrower than that of PBE0, while the interpolated peak of the monotonic increase in HSE06 is also significantly lower than that of PBE0. Though the removal of long-range HF exchange clearly affects energetics, this contribution to energetics is insufficient to prevent the Rutile-Anatase energetic ordering change, as both PBE0 and HSE06 feature a Rutile-Anatase energetic ordering change.

Further analysis of data involving several features linked to the non-monotonic trend formed by varying the fraction of HF exchange are developed in the Supporting Information document for this study. For both polymorphs in both hybrid functionals tested, a summary of these results reveals that the inversely proportional relationship between Rutile and Anatase cell volumes and the Rutile-Anatase formation energy is very strong, illustrating that an analogue of the relationship between atomic structure and relative energetics found when varying $U$ in Hubbard $U$ calculations is also present when varying the HF exchange fraction in hybrid functional calculations cite:Kulik2008,Kulik2011. In the case of Rutile, cell volume first decreases in proportion with the initial monotonic trend in relative formation energy then increases upon reversal of the trend, yielding system volumes at the HF exact exchange limit (\textit{a} = 1.0) that vary across tested functionals. In the case of Anatase, the same increasing and decreasing trends in cell volume are present, though the system volume at the HF exact exchange limit more consistently matches that observed when \textit{a} = 0.25 for both PBE0 and HSE06 functionals.

When reviewing contributions to the expansion and contraction of the Rutile and Anatase unit cells with the addition of HF exchange across both hybrid functionals, a component of cell volume that changes highly proportionately with Rutile-Anatase formation energetics is the \textit{c/a} ratio of Rutile. A similar level of proportionality can be observed in comparing the relationship between the \textit{c/a} ratio of Anatase and HF exchange fraction with the relationship between corresponding Anatase cell volumes and HF exchange fractions. In these relationships, which are depicted in the Supporting Information document, all volumetric information is normalized with respect to the volumetric datum possessing the lowest HF exchange fraction in each plot. This normalization, which is accomplished for each system and each material property plotted, facilitates the direct comparison of different systems and reveals the apparent proportional relationships shared by different systems across the same property. In both proportional relationships, the non-monotonicity originally observed in the Rutile-Anatase formation energy trend as $a$ → 1 is observed, namely as a discontinuity of decreasing magnitude in the \textit{c/a} ratio of Rutile and one of increasing magnitude in the \textit{c/a} ratio of Anatase, the cell volume of Rutile, and the cell volume of Anatase. Therefore, the discontinuity observed in formation energy can be strongly linked to variations in multiple structural properties independent of long-range exchange interaction screening, inferring that the physicality of the energetic discontinuity can be more directly determined via comparisons between the structural features yielded from experimental results and first-principles calculations.

In connection with the discontinuity shared by Rutile-Anatase formation energetics and the \textit{c/a} ratio trends of Rutile and Anatase that both vary with HF exchange fraction, continuous portions of those trends share common characteristics. These shared characteristics include the direct relationship between the \textit{c/a} ratio of Rutile and the Rutile-Anatase formation energy as a function of exact HF exchange, in addition to the inverse relationship between the corresponding \textit{c/a} ratio of Anatase and that formation energy. These relationships directly relate the Rutile-Anatase formation energy with relative changes in \textit{c}-axis length for Rutile and \textit{a}-axis length for Anatase. When considering previous conclusions stating that electron localization on Ti 3\textit{d} orbitals cannot be the only factor influencing Rutile-Anatase formation energetics and that neglecting to consider long-range exact exchange interactions in HSE06 was not sufficient to change relative energetic ordering, the possible link between short-range interactions and changes in the lengths of the shorter axes of both tetragonal polymorphs with corresponding changes in Rutile-Anatase relative energetic ordering can be developed. Previous research indicates that van der Waals dispersion interactions, which are accounted for by short-range energetic contributions proportional to \textit{r}^-6 (\textit{r} represents interatomic distances), are expected to occur in TiO₂ Rutile, Anatase, and other polymorphs cite:Conesa2010,Gerosa2015. Given the poor overlap between TiO₂ Rutile and Anatase experimental cell volumes and cell volumes with experimentally consistent energetic ordering shown in the Supporting Information, cite:Mehta2014 the introduction of both short-range dispersion and electron-electron interaction corrections to simultaneously achieve experimentally consistent energetic and structural information of improved accuracy is strongly supported by analysis completed in this research.

Given the connection between structure and Rutile-Anatase formation energetics in TiO₂ polymorphs, the calculation of formation energetics in similarly structured BO₂ polymorphs using hybrid functionals can also be affected by discontinuities, especially when considering the small magnitude of relative energetic changes resulting from variation in \textit{U} or \textit{a} for TiO₂ polymorphs. When measuring the maximum error that can result from using a PBE or HF non-hybridized functional rather than a hybrid functional, the HF exchange fraction of which is set with respect to the interpolated maximum of each formation energy trend ($a$ → 1), the relative energetic differences between the results produced by either the HSE06 or PBE0 functional with $a$ → 1 and corresponding results produced via a non-hybridized PBE or HF functional range from 0.068 to 0.127 eV/TiO₂. As shown in the Supporting Information and previous research, cite:Mehta2014 epitaxial stabilization of similarly composed polymorphs occurs within an energetic window of 0.1-0.2 eV. Considering that the maximum error reported previously observes the same order of magnitude of the energetic window implemented in epitaxial stabilization applications, the \textit{a priori} use of hybrid functionals as a predictive tool capable of materials selection in these applications is questionable.

Conclusions

In this study, the relative energetic ordering of Rutile, Anatase, Brookite, and Columbite TiO₂ polymorphs has been assessed using DFT+$U$ and hybrid functional methodologies. The relative formation energies of Anatase, Brookite, and Columbite with respect to Rutile were evaluated over variation in the \textit{U}_\textit{3d} parameter on Ti in Hubbard $U$ calculations and the fraction of exact exchange in PBE0 and HSE06 calculations. Past research not incorporating Hubbard $U$ or hybrid functional methodologies of these TiO₂ polymorphs has indicated that their energetic ordering varies with functional selection cite:Mehta2014. However, when incorporating a Hubbard $U$ methodology, energetic trends resulting from variation in $U$ are able to clearly resolve an energetic ordering consistent with experiment that persists regardless of the GGA-based functionals or pseudopotentials employed. In addition to using the LDA functional, this evaluation was completed using the GGA-based functionals PBE, PBEsol, AM05, and PW91. Also, multiple permutations of the $s$-valence inclusive Ti (Ti_sv), $p$-valence inclusive Ti (Ti_pv), soft O (O_s), standard Ti (Ti), and standard O (O) pseudopotentials were considered in this evaluation. Relative energetic orderings consistent with experimental expectations included E_Rutile \textless E_Anatase \textless E_Brookite \textless E_Columbite, E_Rutile \textless E_Brookite \textless E_Anatase \textless E_Columbite, or E_Rutile \textless E_Anatase \textless E_Columbite cite:Swamy2002,Wu2010.

When applying the PBE, AM05, and PW91 functionals in accompaniment with standard or O_s pseudopotentials, an interval of $U$ values between approximately 2.8-4.0 eV demonstrated relative energetic ordering consistent with experiment for all polymorphs studied. Inclusion of the Ti_pv or Ti_sv pseudopotentials on PBE calculations preserved an interval of energetic ordering consistent with experiment that was qualitatively similar to those achieved using standard pseudopotentials. However, in the cases of Ti_pv and Ti_sv inclusive calculations, these energetic orderings differed from their PBE analogues largely by being shifted rightward with respect to them, as these qualitatively consistent energetic intervals occurred at \(U\)=4.7-7.0 and 5.8-8.2 eV, respectively. As a predominant result of upward shifts of the Rutile-Anatase energetic trends within them, the PBEsol functional with standard pseudopotentials observed an energetic ordering consistent with experiment within the $U$ interval of 1.5-2.1 eV, whereas the LDA functional did not produce an experimentally consistent energetic ordering at any value of $U$.

First-principles derived values of $U$, which have been achieved from linear response theory calculations performed in this study and past research cite:Mattioli2010,Xu2015, were applied to determine whether experimentally consistent relative energetic ordering can be predicted in TiO₂ polymorphs. Application of either the first-principles derived $U$ values from this work or corresponding results from other research cite:Mattioli2010,Xu2015 indicates that experimentally consistent results can at least be achieved with the PBE, AM05, and PW91 functionals accompanied by standard pseudopotentials. In past research cite:Finazzi2008,Mattioli2010, an approximate $U$ value of 3 eV has been assigned to Rutile and Anatase polymorph Ti cations to improve the physicality of electronic structure properties relative to experimental results. Application of this $U$ value, or corresponding values calculated via linear response, to the Rutile-Anatase formation energy reveals a result within the range of 0.00 and 0.01 eV/f.u., which is in strong quantitative agreement with past experimental results cite:Dompablo2011,Rao1961.

The Rutile-Anatase relative formation energy, which was most responsible for inconsistencies in relative energetic ordering predictions across functionals, was more extensively analyzed using hybrid functional calculations, namely by calculating formation energy as a function of exact exchange using both the PBE0 and HSE06 hybrid functionals. For PBE0 and HSE06 functionals, experimentally consistent E_Rutile \textless E_Anatase energetic ordering was achieved at fractions of exact exchange of approximately 0.72 and 0.78, respectively, indicating that controlling the charge localization and delocalization in Ti $3d$ orbitals can produce experimentally consistent energetics. However, in the limit of 100% exact HF exchange, Rutile-Anatase energetic ordering was shown to revert to an experimentally inconsistent energetic ordering. For both TiO₂ Rutile and Anatase, structural features observed strong inverse or direct relationships with formation energetics. Given these connections between structure and energetics and the relative energetic error associated with the inability to select a single fraction of exact exchange for both Rutile and Anatase predictively, the arbitrary selection of an exact exchange fraction has greater implications on the prediction of relative energetics in generalized BO₂ systems using hybrid functionals. In combination with the absence of a methodology for predicting a suitable fraction of exact exchange for particular sets of TiO₂ polymorphs, the errors in relative energetic ordering resulting from the use of hybrid functionals on TiO₂ polymorphs suggest care should be taken in their use in predicting relative energetics in less studied BO₂ systems. Minimally, these predictive calculations could require the use of a high fraction of exact exchange not equal to the HF limit to accommodate relative formation energetics similar to the Rutile-Anatase energy, though exact magnitudes of exact exchange cannot be prescribed based on this study of TiO₂ polymorphs.

In each case of varied pseudopotential and functional selection evaluated in this study for Hubbard \textit{U} and hybrid functional inclusive calculations, the quantitative differences between the relative energetics of distinct TiO₂ polymorphs, namely energetics contained by the prescribed ranges of linear response calculated \textit{U} or experimentally consistent energetic ordering, are generally within a range of 0.1 eV of one another. Despite the low magnitude of these energetic differences, the addition of a constant \textit{U} = 3 eV to the 3\textit{d} orbitals of Ti cations has been demonstrated to predict experimentally consistent energetic ordering within the measurement uncertainty associated with \textit{U} calculations. Though relative energetic results resolved in this study are most directly related to internal energy cite:Xu2015 or formation enthalpy values achieved while neglecting pressure involved contributions, cite:Wang2006,Dompablo2011 the consideration of energetic contributions beyond those resolved in first-principles ionic and electronic relaxation calculations is required when applying the \textit{B}O₂ relative energetics yielded in this study to corresponding Gibbs free energy calculations cite:Kitchin2008. Previous information has shown that the consideration of entropy contributions over a wide range (0-1300 K) of temperatures does not affect Rutile-Anatase energetic ordering, as the magnitude of the standard entropy of formation of Rutile is strictly greater than that of Anatase over the evaluated temperature range cite:Smith2009. Past research evaluating the zero-point energy vibrational contributions of TiO₂ Rutile and Anatase, which was calculated at \textit{U} = 0 eV, reveals an approximate difference of 0.02 eV between the two contributions that favors stabilizing Anatase cite:Shirley2010,Moellmann2012. Nevertheless, the calculation of vibrational contributions at finite temperatures using phonon densities of state (also at \textit{U} = 0 eV) reveals, in accompaniment with corresponding experimental data, that a TiO₂ Rutile-Anatase phase transformation occurs between 1200 and 1340 K cite:Mei2014. Though both experimental and computational results imply that vibrational contributions could affect relative energetic ordering in calculated Gibbs free energy results at higher temperatures, these vibrational effects are not anticipated to stabilize Anatase relative to Rutile in DFT calculations and under experimental conditions not involving elevated temperature and pressure conditions.

Recent efforts have been made in the determination of distinct exact exchange fractions suitable for the calculation of properties in individual systems by comparing GW self-energy calculation contributions to those of hybrid functionals, namely via relating the inverted dielectric matrix of the screened non-local exchange term of the GW calculation to the fraction of exact HF exchange contributed to hybrid functional calculations cite:Alkauskas2011,Gerosa2015. Despite the largely experimentally consistent results achieved when extending this approach to compute a pair of related reaction energies involving two systems (the TiO₂ Rutile polymorph and the Ti₂O₃ Corundum polymorph) that each possessed a distinct, calculated fraction of exact exchange, the application of this approach to the relative formation energies involving Rutile, Anatase, and Brookite TiO₂ polymorph pairs yielded an experimentally inconsistent energetic ordering (E_Anatase \textless E_Brookite \textless E_Rutile) cite:Gerosa2015. Nevertheless, appropriate application of the Hubbard $U$ method and self-consistent linear response theory enables consistent, precise prediction of experimental energetic orderings of TiO₂ polymorphs and related systems, encouraging the use of Hubbard $U$ methodologies to predict relative energetic ordering in less studied BO₂ systems.

\begin{acknowledgement} The authors of this paper thank Dr. Giuseppe Mattioli for his contributions to discussions involving the use of the Hubbard $U$ method in this article. The computing resource used to complete this paper was provided by Carnegie Mellon University and its Department of Chemical Engineering. JRK gratefully acknowledges partial support from the DOE Office of Science Early Career Research program (DE-SC0004031). \end{acknowledgement}

\begin{suppinfo} A complete database of the results from this work with examples of using the data to generate the figures are provided in the supporting information. \end{suppinfo}

bibliography:shorttitles.bib,references.bib

Table of Contents Graphic