Stability of non-metal dopants to tune the photo-absorption of TiO2 at realistic temperatures and oxygen partial pressures: A hybrid DFT study

TiO2 anatase is considered to play a significant importance in energy and environmental research. However, for developing artificial photosynthesis with TiO2, the major drawback is its large bandgap of 3.2 eV. Several non-metals have been used experimentally for extending the TiO2 photo-absorption to the visible region of the spectrum. It’s therefore of paramount importance to provide theoretical guidance to experiment about the kind of defects that are thermodynamically stable at a realistic condition (e.g. Temperature (T), oxygen partial pressure (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{p}}}_{{{\bf{O}}}_{{\bf{2}}}}$$\end{document}pO2), doping). However, disentangling the relative stability of different types of defects (viz. substitution, interstitial, etc.) as a function of charge state and realistic T, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{p}}}_{{{\bf{O}}}_{{\bf{2}}}}$$\end{document}pO2 is quite challenging. We report here using state-of-the-art first-principles based methodologies, the stability and meta-stability of different non-metal dopants X (X = N, C, S, Se) at various charge states and realistic conditions. The ground state electronic structure is very accurately calculated via density functional theory with hybrid functionals, whereas the finite T and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{p}}}_{{{\bf{O}}}_{{\bf{2}}}}$$\end{document}pO2 effects are captured by ab initio atomistic thermodynamics under harmonic approximations. On comparing the defect formation energies at a given T and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{p}}}_{{{\bf{O}}}_{{\bf{2}}}}$$\end{document}pO2 (relevant to the experiment), we have found that Se interstitial defect (with two hole trapped) is energetically most favored in the p-type region, whereas N substitution (with one electron trapped) is the most abundant defect in the n-type region to provide visible region photo-absorption in TiO2. Our finding validates that the most stable defects in X doped TiO2 are not the neutral defects but the charged defects. The extra stability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\bf{S}}{\bf{e}}{\bf{O}})}_{{\bf{O}}}^{+{\bf{2}}}$$\end{document}(SeO)O+2 is carefully analyzed by comparing the individual effect of bond-making/breaking and the charge carrier trapping energies.


Results and Discussion
Most stable defect configurations. We  For this purpose, we have adopted an iterative strategy 46 : starting with the Ti 16 O 32 supercell, we first identified the most favorable sites for the O-□ by scanning over all possible sites. We find that, since they belong to the same symmetry point, all the atomic sites are equivalent. Hence, we can substitute non-metal dopants at any O-site in the system. For interstitial and complex interstitial configurations, we simply consider one oxygen atom in the Ti 16 O 32 framework and scan its all surrounding positions. The configuration having minimum energy is taken as a final structure for studying defects. Note that for all geometry relaxations we have used PBE functional as the difference in geometry optimized by PBE and HSE06 are insignificant. Therefore, on top of the final relaxed geometry, we have performed HSE06 calculations to determine the correct energetics associated with those optimized structures.
Validation of DFT functionals. To ensure that our findings are not an artifact of the chosen treatment of the xc-functional, we have first thoroughly benchmarked the xc-functionals. Note that using PBE, the band gap of TiO 2 is quite underestimated from the experimental value [PBE (2.12 eV), Expt. (3.2 eV)]. On varying the "Hubbard parameter U" one can reproduce the experimental band gap with PBE + U approach. With U = 6.3 eV, we could reproduce the band gap to be 3.2 eV. However, the energetics of PBE + U to estimate defect formation energy of a single oxygen vacancy at various charge states doesn't match with respect to advanced hybrid xc-functional HSE06 (see Supplementary Information; Figs S1 and S2). Note that, the default α = 0.25 value of HSE06 is unable to reproduce the exact experimental band gap. We find that a proportion of 22% Hartree Fock exchange with 78% PBE exchange produces the experimental band-gap (3.2 eV) for HSE06. In view of this, the accuracy of HSE06 with α = 0.22 needs also to be validated. We thus further validate HSE06 xc-functional (with α = 0.22) with respect to experimentally observed optical transition of various (X) O (X = N, C, S, Se) (i.e. X substituted at oxygen vacancy defects) in TiO 2 . The optical transitions associated with the process ν are shown in (Fig. 2) by forming a configuration co-ordinate diagram. To do this, we have taken 3 different geometries of X O (optimized for charges = 0, +1, +2). Then, we have taken intermediate structure between geometry 1 and geometry 2 (simply by averaging the atomic positions of the geometry 1 and 2). The same is then repeated for the geometry 2 and geometry 3. In this way, we have achieved the 5 intermediate relaxed structure. We then linearly interpolate the relaxed structures in the ground and the constrained state by using five intermediate structures and calculate the formation energies using HSE06 xc-functional (with α = 0.22). A parabolic fit of the formation energies of these five intermediate geometries is thus obtained. The same is then repeated for the other charge states (+2 is not shown here). Now the possible transitions are illustrated as shown in Fig. 2. Knowing the optical transitions, the corresponding optical spectra can thus be computed. The details of this state-of-the-art theoretical spectroscopy techniques are explained by Rinke et al. 47 . The process involved absorption of photon and emission of electron. The Peak positions are determined by assuming that the atoms in the vicinity of the defect do not have time to relax during the transition from charge state q to charge state q′ = q ± 1(Franck-Condon principle). Thus, the difference in formation energies in charge states q and q′ gives the absorption and emission energy. The difference in the value of absorption and emission peaks arises due to difference in geometry relaxation for both the charge states q and q′. The amount of energy lost during relaxation to the new equilibrium position is the Franck-Condon shift. We have found that X O leads to sub-bandgap optical transitions centred at 2.8 eV (442 nm) for N, 2.04 eV (607 nm) for C, 2.32 eV (534 nm) for S and 1.92 eV (645 nm) for Se respectively. These results are in good agreement with the experimental observations in which sizeable absorption occurs at visible region (400-700 nm) [48][49][50][51] . Note that, we have also tried the same for Phosphorous substitution at O-site in TiO 2 anatase. Unfortunately, we have got the response in infra-red region (1.48 eV) (see Supplementary Information; Fig. S6). In view of this, we have not considered it for further study. Therefore, this validates that our HSE06 xc-functional with (α = 0.22) is sufficiently accurate to ensure correct energetics at various charge defects in TiO 2 . www.nature.com/scientificreports www.nature.com/scientificreports/ Formation energy of charged defects. We have now considered three types of defects viz. (X) O , (XO) O and (X 2 ) O with X = N, C, S, Se to understand the most stable defect with specific charge state q at a given T and p O 2 . The explicit dependence on T and p O 2 is introduced by employing the concept of ab initio atomistic thermodynamics 39 . We address the stability of the substitutional and interstitial configurations of N, C, S, and Se dopants in TiO 2 by computing formation energies as a function of chemical potentials of the respective species viz. μ O , μ N , μ C , μ S , μ Se as discussed in the following sections. Here the doping is considered as varying the chemical potential of electrons μ e . In the following sections, we will use E [(X) ] N-related dopant. The formation energy of (N) O in charge state q is given by: [see details of this methodology in ref. 52 ]. where m is the mass,  is the spin multiplicity, and σ is the symmetry number.
. Solving this we get, under O-poor condition Δμ Ti = −1.51 eV and Δμ O = −4.11 eV. Defect formation energy depends on the chemical potential of Ti, O and the doped impurities. Note that to avoid the precipitation of the host elements, the chemical potential μ X must have a constraint that Δμ X ≤ 0. We find that under O-poor conditions (i.e. Ti-rich) Δμ X = 0 and Δμ O = −4.1 eV. Under O-rich conditions Δμ O = 0 and μ N is estimated by the formation of N 2 molecule, i.e., Δμ N = 0; while the limit of the other impurities X (excluding N) (X = C, S, Se) are determined by the formation of the corresponding oxides i.e., The chemical potentials Δμ X for X = C, S, Se are, therefore, determined by formation enthalpy of CO 2 (ΔH f = −4.07 eV), SO 3 (ΔH f = −4.101 eV) and SeO 2 (ΔH f = −2.336 eV) respectively 54 . Note that we have chosen only those oxides which are used as precursor in experiments. With this background, we report now the most stable charged impurities for various non-metals X.  www.nature.com/scientificreports www.nature.com/scientificreports/ The intersection of the formation energies [e.g. 0/−1, +1/−1, +2/−1 etc.] for the different charge states correspond to the thermodynamic transition levels. Note that N has one electron less than O in the valence shell. Therefore, (N) O in TiO 2 is expected to act as an acceptor level. It can be stable either in (N) O 0 or − (N) O charge states. At a lower μ e i.e., near valence band maximum (VBM) (N) O 0 is stable and for a larger value of μ e (i.e., in the upper part near the conduction band minimum (CBm)) − (N) O is stable (see Fig. 3). Since, Nitrogen will act as an electron acceptor because it has one electron less than oxygen, it is more feasible at higher values of μ e i.e., n-type doping. Our key conclusion from calculated formation energies indicate that (N) O  The same can be calculated directly from DFT total energy from a supercell of C 8 . The equations related to the formation energy of C related defects in charge state q are given in Supplementary Information. In Fig. 4 Fig. 4) for all the three conditions. Our main concern is near CBm as it is an acceptor dopant.
Note that for interstitial C in TiO 2 i.e., (CO www.nature.com/scientificreports www.nature.com/scientificreports/ configurations using either S or SO 3 as a precursor. In either case we have verified that the results are same. To use SO 3 as a precursor, we have used the following formula to estimate E(S): The same can also be calculated directly from DFT total energy from a supercell of S 8 (see equations in Supplementary Information). In  www.nature.com/scientificreports www.nature.com/scientificreports/ (less Ti atoms), the calculated formation energy shows that interstitial S (SO) O is more prevalent defect than other configurations.
Se-related dopant. The formation energy formula for Selenium related defects for all the three configurations using SeO 2 as a precursor is given in Supplementary Information. The calculated Formation energy is shown in Fig. 6  Note that, from all the above explanation it is cleared that, in order to stabilize the system in case of doping, we need to supply external charges into the system. This means if we assume of having reservoir of electrons (from where the dopant can trap one electron), the dopant will be stabilized at the defect site. This is as good as claiming the presence of having several oxygen vacancies in the lattice with uncompensated electrons and those electrons can be trapped by the external dopant. Thus we can adapt either of the approaches i.e. (1) addition of external charge to the dopant or (2) explicit presence of oxygen vacancy for the neutral dopants. In the literatures, www.nature.com/scientificreports www.nature.com/scientificreports/ people have adopted the methodology in which charge compensation are done to the (neutral) dopant by creating O-vacancies 24,25,[59][60][61][62] . However, here we have adopted a different approach, where the charge compensation are performed by adding or removing the extra electrons/holes to the dopant (charged defects). Both the approaches yield to the same conclusion as they are doing effectively the same charge compensation at the defect site and therefore our results are well in agreement with the previous findings of the Pacchioni's group 24,25,[59][60][61][62]  is carefully further analyzed quantitatively. For this purpose, we have decomposed the overall defect formation energy into two parts: (i) bond-making/breaking energy (Δε) and (ii) charge carrier trapping energy (Δζ). The combined effect of these two competing factors decides which defect configurations will get more stabilized over other configurations. Now Δε is given by: Note that μ e is varied from VBM (i.e. highest possible p-type doping) to the CBm (i.e. highest possible n-type doping). Here to understand the stability at a given μ e , we have fixed our μ e at VBM of pristine TiO 2 at charge state www.nature.com/scientificreports www.nature.com/scientificreports/ 0. The VBM of pristine Ti 16 O 32 is taken as a reference for calculating charge carrier trapping energy. Therefore, the charge carrier trapping energy i.e., Δζ is given by (see Fig. 8): It should be mentioned here that Eq. 12 is explicitly for p-type doping. Note that, for n-type doping one needs to set a different μ e as reference level in the n-type region (viz. CBm(TiO 2 ) 0 ). From Table 1    . A similar explanation can also be drawn for the stable compositions at other environmental conditions viz. O-poor and O-rich as shown in (see Supplementary Information; Fig. S5).

conclusion
In summary, we have presented an exhaustive study to understand the most stable non-metal related defect in bulk form of anatase TiO 2 as a function of charge at a realistic condition. As a first step, we have validated our DFT functionals, where we find that involvement of hybrid functional HSE06 is essential and semi-local functional and its improved variants are not sufficient even to address this problem qualitatively. We further notice that the α parameter of HSE06 needs to be adjusted to 0.22 so as to correctly predict the experimental bandgap of 3. The enhanced stability is also carefully analyzed quantitatively by including the individual effect of bond-making/breaking and the charge carrier trapping energies.

computational Details
We have performed the DFT calculations with PAW pseudopotential method 63 as implemented in Vienna ab initio simulation package (VASP) 64 . The supercell, consisting 48 atoms, is constructed by 2 × 2 × 1 replication of the tetragonal TiO 2 unit cell (space group number:141, I41/amd). The supercell size ensures enough spatial separation between the periodic images of the doped impurities under periodic boundary conditions. We have explicitly checked this by making a single O-vacancy to be fully localized inside the cell. To ensure the supercell size convergence, test calculations using 96-atom supercell have also been performed for the case of O-vacancy. The results are consistent with that obtained from a 48-atom supercell. However, a drawback of limited system sizes is that the dopant concentrations are artificially higher (perhaps two to three times) than the experimental case 11 . In order to ensure that our findings are not just an artifact of DFT functionals, we have used a number of exchange and correlation (xc) functionals viz. generalized gradient approximation (GGA) with PBE 65 , GGA + U 66 (where U is Hubbard parameter) and hybrid functional HSE06 67,68 . The results given here are from HSE06 xc-functional. The performance of the other xc-functionals are discussed in detail in Supplementary Information; Figs S1 and S2. All the structures are fully relaxed (atomic position) upto 0.001 eV/Å force minimization using conjugate gradient minimization with 4 × 4 × 2 K-mesh. For electronic structure energy calculations, the brillouin zone is sampled with a 8 × 8 × 4 Monkhorst-Pack 69 K-mesh with 0.01 meV energy tolerance. In all our calculations, the plane wave energy cut-off is set to 600 eV.