Nonmetal Doping at Octahedral Vacancy Sites in ... - ACS Publications

Jun 29, 2007 - Nicholas C. Wilson,† Ian E. Grey,*,† and Salvy P. Russo‡. CSIRO Minerals, BayView AVenue, Clayton, Victoria 3169, Australia, and ...
7 downloads 0 Views 201KB Size
Download PDF
J. Phys. Chem. C 2007, 111, 10915-10922

10915

Nonmetal Doping at Octahedral Vacancy Sites in Rutile: A Quantum Mechanical Study Nicholas C. Wilson,† Ian E. Grey,*,† and Salvy P. Russo‡ CSIRO Minerals, BayView AVenue, Clayton, Victoria 3169, Australia, and Applied Physics, School of Applied Sciences, RMIT UniVersity, Melbourne 3001, Australia ReceiVed: March 22, 2007

The stability of boron and carbon dopants (D) at vacant octahedral sites in rutile has been investigated using density functional theory quantum mechanical (QM) modeling. Three different types of sites were considered: vacant Ti lattice site (Dvac), vacant Ti site + Ti interstitial (DFrenkel), and interstitial site (Dint). The defect formation energies, Ed, at different temperatures and gas partial pressures were calculated from the QM total energies of the relaxed defect and host structures, combined with chemical potentials of the reservoir components that were corrected for both temperature and gas partial pressure variations. The contribution of vibrational free energy to Ed as a function of temperature was evaluated for the BFrenkel model using ab initio phonon density of states calculations. The calculated vibrational and configurational free energy terms were of opposite sign and partially cancelled, giving a relatively small (0.3 eV at 700 K) combined contribution to Ed. Under strongly reducing conditions at 1500 K, boron incorporation at interstitial and Frenkel sites is favored, with Ed values of 0.53 and 0.8 eV respectively, whereas the Ed values for carbon doping were high (>4 eV) for all three models under high-temperature reducing conditions. Under oxygen-rich conditions relevant to sol-gel processing, the Dvac model was favored for both boron and carbon. Further stabilization of the Dvac model for boron was obtained at a protonated vacancy site, giving Ed ) 1.16 eV for (B+H)vac at 700 K.

1. Introduction In the vast body of literature on the low-temperature photocatalytic properties and high-temperature electrical properties of TiO2, the role of titanium vacancies, VTi, has received little attention until recently. Generally, the properties are considered only in terms of oxygen vacancies and titanium interstitials.1 Mizushima et al.2 appear to be the first researchers to evoke VTi as a native defect in rutile to explain their photocurrent spectra results. Recently Nowotny et al.1,3 made high-temperature conductivity and thermoelectric power measurements on single-crystal rutile during prolonged oxidation experiments and were able to separate the influences of VTi defects from those due to oxygen vacancies and Ti interstitials. Phattalung et al.4 applied first-principles quantum mechanical (QM) calculations to study different native defects in anatase-form TiO2 and concluded that titanium vacancies were the lowest-energy acceptor defects under oxygen-rich growth conditions. In sol-gel prepared TiO2 for photocatalytic applications, direct evidence for VTi defects has been obtained from structure refinements.5-8 These studies show that, particularly for anataseform TiO2, the VTi concentration increases with decreasing crystallite size. Our recent studies,7 using the Debye function to model the diffraction properties of small TiO2 clusters containing VTi defects, are supportive of the VTi defects being distributed throughout the nanocystal rather than being concentrated at the surface. QM modeling7 showed that the VTi defects were stabilized by the incorporation of protons at the vacancy sites, analogous to the Ruetschi defects in MnO2.9 Given the stable character of VTi defects in TiO2, it might be expected that they would play a role in providing sites for dopant * To whom correspondence should be addressed. † CSIRO Minerals. ‡ RMIT University.

incorporation. The doping of TiO2 with nonmetals such as C, N, F, and B has provided a great impetus to photocatalysis research since the demonstration by Asahi et al.10 that such doping has the potential to improve the photocatalysis efficiency by modifying the band gap and extending the light absorption by TiO2 into the visible region. Asahi et al.10 complemented their experimental studies on N-doping with first-principles QM calculations of the changes to the band gap properties. Their work has been followed by numerous QM modeling studies of nonmetal doping of TiO2.11-18 Most of these studies,11,13,15-17 have focused on the incorporation of the dopants at oxygen vacancy and interstitial sites. Exceptions are the studies by Geng et al.14 on B incorporation in anatase and Valentin et al.12 on C doping of rutile and anatase, in which substitution at an empty Ti site was considered. In the latter two studies using density functional theory (DFT), the defect formation energies were calculated from the total energies of relaxed pure and defect structures in different sized supercells. An important aspect of such calculations is choosing appropriate chemical potential reservoirs for the supply of the dopant atoms and the acceptance of the O or Ti atoms that are replaced by the dopant atoms. Geng et al.14 used metallic titanium for the Ti reservoir so their results are relevant only to doping under strongly reducing conditions. Valentin et al.12 conducted a more comprehensive study in which the oxygen chemical potential was varied, allowing them to determine the relative stability of different modes of incorporation of C under both oxygen-rich and oxygen-poor conditions. They showed that incorporation of C at VTi sites was favored at the oxygen-rich conditions relevant to sol-gel synthesis conditions. Our interest in nonmetal doping of TiO2 goes back to an earlier study of rutile growth from borate fluxes in a reducing atmosphere at temperatures in the range 1373-1573 K.19 We

10.1021/jp072299e CCC: $37.00 © 2007 American Chemical Society Published on Web 06/29/2007

10916 J. Phys. Chem. C, Vol. 111, No. 29, 2007

Wilson et al.

found that significant levels of B (up to 2.4 atom%) were incorporated from the flux into the rutile crystals. A structural model for the boron incorporation at VTi sites was proposed, based on crystal chemical reasoning and chemical analyses. Direct confirmation of the model by a structure refinement was not possible because of the low scattering power of the B and the lack of long-range order. QM modeling provides an opportunity for testing the validity of the dopant incorporation model through the calculation of the defect formation energies to determine if the defect is stable. In this study, we present the results of DFT calculations of defect energies for B and C doping at empty octahedral sites in rutile-form TiO2. Three different cases are considered: (1) substitution for titanium at a VTi lattice site, (2) incorporation at a VTi site created by displacement of a Ti atom to an adjacent interstitial site (Frenkel defect), and (3) location within an interstitial site. The three different models for dopant incorporation are designated as Dvac, DFrenkel, and Dint, respectively, where D ) dopant (B or C). For the case of low-temperature (solgel) doping the incorporation of dopants at VTi sites stabilized by protons is considered. 2. Methodology 2.1. Computational Details. Quantum mechanical total energy calculations were based on the DFT GGA approximation and were implemented using the VASP code.20 The PW91 GGA exchange-correlation functional was used and the ultrasoft Vanderbilt pseudopotentials provided with VASP. A planewave cutoff energy of 396 eV was used, and the k-space integrations were performed using a Monkhorst-Pack sampling scheme,21 with a 7 × 7 × 7 k-point mesh. One dopant atom was included per 2 × 2 × 2 rutile supercell for the majority of the calculations. However, calculations were also performed in 2 × 2 × 3 and 3 × 3 × 3 rutile supercells for selected models to check the effect of dopant-dopant interactions on the total energy. The calculations involved relaxing the structures using a conjugategradients scheme to minimize the forces on all atoms and to obtain minimum energy 0 K configurations for both pure and doped rutiles. The relaxations were terminated when the forces on all atoms decreased below 0.008 eV/Å. Minimum energy relaxations were also performed on the different phases and components involved in the chemical potential reservoirs. The latter included crystalline TiO2, B2O3, Ti, B, and graphite and molecular O2, H2, H2O, and CO2. The molecular species were contained in periodic cells large enough to avoid moleculemolecule interactions. The dopant formation energies, Ed, at 0 K and 0 bar pressure were then calculated from Ed ) Edefect - Ehost + nTiµTi - nDµD

(1)

where Edefect and Ehost are the calculated total energies of the defect supercell and the host (TiO2) supercell, nTi and nD are the number of Ti atoms removed and the number of dopant atoms added, and µTi and µD are the chemical potentials of Ti and the dopant in their respective reservoirs. All results are reported in relation to molar quantities. Calculations of vibrational free energies at elevated temperature for pure rutile and the BFrenkel model were made using VASP for electronic relaxations and for the calculation of the Hellmann-Feynman forces, coupled with PHONON.22 The calculations were conducted using 2 × 2 × 2 supercells. For pure rutile, the full symmetry of the rutile structure was utilized, whereas for doped rutile, no symmetry was assumed (spacegroup P1). For the doped sample, the 2 × 2 × 2 cell contains 16 Ti,

32 O, and 1 B. Each atom was displaced in turn by 0.03 Å along each of the three independent directions, and the electronic structure was relaxed. This entailed 147 separate ab initio electronic relaxations and total energy calculations using VASP. The VASP conditions were the same as those described above for the total energy calculations. The forces were used in PHONON to calculate the total phonon density of states (DOS). Full details of these procedures for the case of pure rutile have recently been reported.23 Due to the high computational expense of these calculations, they were only performed for the one doping model. 2.2. High-Temperature Modeling. Expression 1 gives 0 K, 0 bar, and constant volume dopant formation energies. However, the experimentally observed incorporation of boron into rutile occurs at relatively high temperatures, >1300 K, and at low oxygen partial pressures.19 For the DFT results to be relevant to the experimental results, the effects of temperature and pressure on the different components of expression 1 need to be evaluated. The chemical potential terms in expression (1) come directly from QM total energy calculations at T ) 0 K, p ) 0 bar. These can be equated to µ(0, p0) chemical potentials at the standard pressure p0 of one bar (the effect of 1 bar of pressure on the interatomic distances and hence on E(0,0) is negligible). The µ(0,p0) values can then be corrected to the experimental temperature, µ(T, p0) by adding the increase in the Gibb’s free energy per mole from 0 K to temperature T, ∆G(∆T,1). This is obtained directly from tabulated thermodynamic data.24 In the case of gaseous components, the correction to the experimental pressure, p, is obtained using an additional term of the type RT ln(p/p0). For oxygen and carbon dioxide, the complete expressions, for po ) 1 bar are 1 1 µO(T,p) ) µO2 ) [EO2(0,1) + ∆GO2(∆T,1) + RT ln pO2] (2) 2 2 µCO2(T,p) ) ECO2(0,1) + ∆GCO2(∆T,1) + RT ln pCO2

(3)

For modeling of dopant formation energies at elevated temperatures, the first two terms on the RHS of (1) should strictly be free energies at constant pressure, which are related to the DFT energies, E(0,0), by expression 4 G(T,p) ) E(0,0) + Fvib + Fconf + pV

(4)

where Fconf and Fvib are the configurational and vibrational Helmholtz free energies at temperature T. A simple dimensional analysis analogous to that given by Reuter and Scheffler25 shows that the pV term will be of the order of tens of meV at the 1 bar pressure used in this study and so its contribution can be ignored. The vibrational term is expected to be significant at the high temperatures relevant to the B-doping experiments. In the harmonic approximation, the contribution from each vibrational mode of angular frequency ωi is given by 1 Fvib(T,ωi) ) pωi + RT ln(1 - epωi/RT) 2

(5)

This is convoluted with the vibrational DOS to obtain the Helmholtz vibrational free energy. We used the direct method22 to determine Fvib for pure rutile and for one of the B-doping models (the BFrenkel model). We found that Fvib(defect) - Fvib(host) is of a similar magnitude but of opposite sign to the corresponding configurational free energy term, and so the two contributions cancel to some extent, the difference being 0.32 eV at 700 K, rising to 0.57 eV at

Nonmetal Doping in Rutile

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10917

1500 K. We would expect a similar cancellation of Fvib and Fconf terms for the other models for both B and C doping as all of them involve short D-O bonds with similar high-frequency contributions to the phonon DOS and have similar configurational entropies, and so neglecting the configurational and vibrational free energy terms is not expected to change the relative stabilities of the different models. Thus in the presentation of the defect formation energies, Ed(T,p) will be calculated using expression 1, where the chemical potential terms are corrected for temperature and pressure but the Edefect - Ehost terms are from the 0 K total energy calculations. In the case of the BFrenkel model, the results will be given for the calculations of the vibrational and configurational energy terms. 2.3. Chemical Potential Reservoirs. Care was needed in setting chemical potential reservoirs and oxygen chemical potential ranges that would be appropriate to the different cases of high-temperature reducing conditions and low temperature sol-gel oxygen-rich conditions. The reservoirs used for Ti, H B, and C were TiO2, H2O, B2O3, and CO2. Within the stability fields of these bulk phases, the following chemical potential relations hold: µTi + 2µO ) µTiO2

(6)

2µH + µO ) µH2O

(7)

2µB + 3µO ) µB2O3

(8)

µC + 2µO ) µCO2

(9)

At the high temperatures relevant to the experimental Bdoping,19 equilibrium will be established between the different gaseous and solid species, and the chemical potentials in expressions 6-9 will be interdependent. This contrasts with a previous study on C-doping of TiO2,12 aimed at a lower applications temperature of 700 K, where CO2 was considered as an independent reservoir not in equilibrium with O2. As we are interested in evaluating dopant stabilities both under sol-gel conditions and under high-temperature reducing conditions, oxygen chemical potentials appropriate to the two situations need to be set. We define a reference state for the oxygen chemical potential, µ′O(T,p) by 1 µ′O(T,p) ) µO(T,p) - [EO2(0,1) + ∆GO2(∆T,1)] 2

(10)

Then, from expression 2, µ′O(T,p) ) 1/2 RT ln pO2 ) 0 when pO2 ) 1 bar. This is the oxygen-rich limit used to represent sol-gel preparative conditions. This reference poses a minor problem because the pure (p ) 1 bar) oxygen atmosphere precludes the presence of other gaseous species like CO2 and H2O as chemical potential reservoirs. The problem can be circumvented if the maximum oxygen partial pressure is set slightly below one. For C-doping, an appropriate oxygen-poor limit for µ′O(T,p) is the value corresponding to the onset of solid carbon deposition. This is determined from the two equilibria 2CO(g) ) C(s) + CO2(g)

(11)

2CO2(g) ) 2CO(g) + O2(g)

(12)

With pO2 + pCO + pCO2 ) 1, the partial pressures of the three gaseous components are uniquely determined at a fixed temperature, T, when the activity of C (s) ) 1. From the value of pO2 thus calculated, the oxygen-poor limit for µ′(T,p) can

Figure 1. Model showing the incorporation of a dopant atom, D, in rutile, at a vacant titanium site created by the displacement of Ti to an adjacent interstitial site (Frenkel defect). The rutile structure is viewed along . The c axis is horizontal.

be derived from (2) and (10). At 1500 K, the calculated oxygen partial pressure for C-deposition is 10-17 bars and the corresponding value of the oxygen-poor chemical potential limit is µ′O(1500, 10-17) ) -2.5 eV. For B-doping, the experimentally determined lower oxygen partial pressure limit for B-stabilized rutile at 1500 K occurs at a pO2 of 10-15 bars.19 The corresponding value for µ′O (1500,10-15) is -2.23 eV. We will use the same limit of -2.5 eV for µ′O as used for C-doping to plot the Ed vs µ′O results. 3. Structural Models for Dopant Incorporation The model originally proposed for B incorporation into rutile19 was based broadly on a model proposed by Bursill and Blanchin,26 to explain the incorporation of cations into vacant octahedral sites (interstitial sites) in rutile. According to this model, interstitial ions such as Ti3+ ions diffuse along the [001] channels of vacant octahedral sites in the rutile structure. At any one interstitial site along the channel, the presence of an extra ion gives rise to three face-shared octahedra along , Ti-Ti-Ti. The strong cation-cation repulsive forces are relieved by a terminal Ti atom in the linear clusters displacing from its lattice site to an adjacent interstitial site, giving two pairs of face-shared octahedra. The defect can be considered as a combination of a dopant atom incorporated at an interstitial site and a Frenkel defect. At high dopant concentrations, the defect clusters are ordered into planes, forming crystallographic shear structures.26 In the case of B incorporation, the small B3+ ion has a preference for three or four coordination. Rather than occupy the center of the vacant octahedral site, it displaces to one face, where it can adopt BO3 triangular coordination to three undersaturated oxygen atoms that are each coordinated to only two Ti atoms. The incorporation of B at a triangular site, and the displacement of Ti to an adjacent interstitial site (the Frenkel defect) is illustrated in Figure 1. To maintain charge balance, the addition of B3+ to TiO2 is accompanied by the reduction of three Ti4+ ions to Ti3+. The ratio of B:Ti3+ ) 1:3 required for charge balance was confirmed by chemical analyses on B-doped rutiles.19 The local structure of the defect cluster corresponds to a small element of the calcite, CaCO3 structure, which in fact is adopted by titanium borate TiBO3.27 A similar model can be evoked for C incorporation into rutile, locally forming triangular CO3 groupings. Small dopants such as B or C may also be incorporated into the vacant octahedral interstitial sites in the [001] channels in rutile, without the accompanying formation of a Frenkel defect. If B and C remain in their maximum oxidation state, then, formally, there needs to be a corresponding reduction of 3 and 4 Ti atoms respectively to the trivalent state for charge balance.

10918 J. Phys. Chem. C, Vol. 111, No. 29, 2007

Wilson et al.

TABLE 1: Chemical Potentials (eV) at Different Temperatures and at 1 Bar Pressure for Gaseous Species µ (0 K) ) E(0,0)

∆G (0-700 K)

µ (700 K)

∆G (0-1500 K)

µ (1500 K)

O2 (molecule) TiO2 (crystal) B2O3 (crystal) CO2 (molecule) H2O (molecule) Ti (metal) B(metal) C (graphite) H2 (molecule)

-9.80 -26.78 -40.48 -22.98 -13.99 -7.73 -6.71 -9.24 -6.75

-1.55 -0.42 -0.46 -1.62 -1.41 -0.28 -0.07 -0.07 -1.01

-11.35 -27.20 -40.94 -24.60 -15.40 -8.01 -6.78 -9.31 -7.76

-3.59 -1.56 -2.22 -3.89 -3.37 -0.85 -0.33 -0.28 -2.40

-13.39 -28.34 -42.70 -26.87 -17.36 -8.58 -7.04 -9.52 -9.15

Ti (in TiO2) B (in B2O3) C (in CO2) H (in H2O)

-16.98 -12.89 -13.18 -4.55

1.13 0.93 -0.07 -0.32

-15.67 -11.96 -13.25 -4.87

2.03 1.58 -0.30 -0.79

-14.95 -11.31 -13.48 -5.33

Finally, with respect to doping at vacant Ti sites, a third model can be proposed in which the dopant replaces Ti at a lattice site. For C-doping, there is no charge imbalance for C4+ replacing Ti4+, whereas for B3+ doping, some means of charge balance is required. We have previously shown that vacant Ti sites in anatase are stabilized by the incorporation of protons,7 and in this study, we explore the incorporation of protons at VTi sites in doped rutile. 4. Results and Discussion 4.1. Chemical Potential Results. Table 1 presents the DFT total energies per mole at 0 K for the various chemical potential reservoirs and the natural-state elements, as well as the ∆G(∆T,1) add-on terms and the resulting chemical potentials at temperatures of 700 and 1500 K. The former temperature is typical for the incorporation of dopants into sol-gel prepared TiO2,12 whereas the latter is representative of the hightemperature conditions used for B-doping of rutile.19 The E(0,0) values as shown in Table 1 are rarely presented in QM modeling publications because they depend on the QM method used. However, the use of the VASP code for titania studies is widespread enough4,28-30 to warrant presenting the raw data here for the benefit of other researchers. The DFT E(0,0) values in the upper part of Table 1 can be used to calculate the free energies of formation at 0 K for the different reservoir phases and compare them with published values. For example ∆Gf(TiO2, 0 K)bulk ) E(TiO2)bulk E(Ti)metal - E(O2)mol· in the case of TiO2. The calculated ∆Gf values for TiO2, H2O, CO2, and B2O3 of -9.25, -2.34, -3.94, and -12.35 eV, respectively, differ by typically ∼5% from the published experimental values24 of -9.80, -2.49, -4.08, and -13.26 eV, respectively. The results in the upper part of Table 1 were used to calculate the chemical potentials at 0, 700, and 1500 K for Ti and for the different dopants, D, in their respective reservoirs. The results at the oxygen-rich limit are shown in the lower part of Table 1. They show that the chemical potentials for Ti and for dopant elements B, C, and H have quite different temperature responses. For example µTi becomes more positive by 2.03 eV between 0 and 1500 K, whereas µH becomes more negative by 0.79 eV. The influence of these changes on the defect formation energy, Ed, for the different dopant models is readily obtained by substituting the values of µTi and µD at the different temperatures into expression (1). This leads to the following general observations on the variation of Ed with temperature for the different dopant models at the oxygen-rich limit. • The defect formation energy for the Dvac model becomes more positive (less stable) with increasing temperature. The magnitude of the destabilisation is much greater for C-doping

(+2.3 eV from 0 to 1500 K) compared with B-doping (+0.45 eV from 0 to 1500 K). • The defect formation energies for the DFrenkel and Dint models become more negative (more stable) with increasing temperature for B-doping but become slightly less stable for C-doping. • Incorporation of protons becomes less favorable with increasing temperature. It should be emphasized that the above trends concern only the influence of the temperature-dependent chemical potential terms on Ed. The overall effect of temperature on Ed will depend also on the contributions of the vibrational and configurational free energies. 4.2. Results for B-Doping. 4.2.1. Relaxed Structures. Polyhedral representations of the relaxed structures, viewed along [001] are shown for the Bint, BFrenkel, and Bvac models in Figures 2-4, respectively. In the case of the Bvac model, Figure 4 corresponds to the case where one H atom was included for charge balance (B3+ + H+ T Ti4+). Bond lengths and angles associated with the B-dopant environments in the three models are reported in Table 2. Relaxation of the Bint model, in which the B was located at the center of an octahedral interstitial site, resulted in movements of the surrounding oxygen atoms O(1) to O(4) to give 4-coordination of the B, with a mean B-O distance of 1.493 Å. This can be compared with a value of 1.478 Å reported as the mean tetrahedral B-O distance in borates, averaged over 120 observations.31 The O-B-O angles in Table 2 correspond to a flattened tetrahedron, with 4 angles of 102.9° and two

Figure 2. Relaxed structure for boron incorporated at an interstitial site in rutile (Bint). The structure is viewed along [001]. The atom labeling corresponds to that used in Table 2.

Nonmetal Doping in Rutile

Figure 3. Relaxed structure for boron incorporated at a Frenkel defect site in rutile (BFrenkel). The structure is viewed along [001]. The atom labeling corresponds to that used in Table 2.

Figure 4. Relaxed structure for boron incorporated at a vacant Ti lattice site in rutile (Bvac). The structure is viewed along [001]. The atom labeling corresponds to that used in Table 2.

opposing, much wider angles of 123.8°. The oxygen atoms O(1) to O(4) are each coordinated to three Ti atoms as well as the B atom. These oxygen atoms are thus formally over-saturated and Table 2 shows that their Ti-O bonds, to Ti(1) through Ti(8), are all considerably longer than the mean 〈Ti-O〉 distances. In addition, the mean 〈Ti-O〉 distances for Ti(1) to Ti(8) are in the range 1.99 to 2.00 Å, and are longer than the mean distance of 1.98 Å for the other 8 Ti atoms in the supercell. The observations are consistent with trivalent titanium being more concentrated at Ti sites adjacent to the interstitial B atom. The relaxed BFrenkel model, shown in Figure 3, has the B atom in triangular coordination, to oxygen atoms O(1) to O(3), with a mean 〈B-O〉 distance of 1.395 Å. This compares closely with the mean 〈B-O〉 distance of 1.394 Å in TiBO3,27 which has the same atomic arrangement as in the local environment around the dopant in the BFrenkel model.19 The full coordination octahedron is given in Table 2 for the Frenkel-displaced titanium atom Ti(1) and for Ti(2) with which the Ti(1) octahedron shares a common face. Ti(1) and Ti(2) are displaced in opposite directions by ∼0.3 Å normal to their common shared face, giving a Ti(1)-Ti(2) distance of 2.85 Å. The Bvac model is illustrated in Figure 4. In contrast to the other two models, which require trivalent titanium for charge balance, this model contains only tetravalent titanium. The

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10919 resulting shorter Ti-O distances, except for the strongly distorted Ti(1), are evident from an inspection of Table 2. The B has the same triangular coordination as in the BFrenkel model, with the same mean 〈B-O〉 distance of 1.395 Å. Figure 4 shows the location of the proton which was added for charge balance. The relaxed O-H distance of 0.99 Å is typical for such a bond. 4.2.2. 1500 K Defect Formation Energy Results. The relative stabilities at 1500 K of the three different B-doping models are shown in the plots of Ed vs µ′O in Figure 5. As discussed in section 2.2, the Ed values corresponding to the lines in Figure 5 are derived from the E(0,0) total energies for the host and doped rutile at 0 K, combined with the chemical potentials, µTi and µB, corrected to the temperature and oxygen partial pressure of the experiments. The values of the log of the oxygen partial pressure corresponding to the µ′O values are shown by the scale at the top of Figure 5. At the oxygen-rich limit at 1500 K, the Bvac model has a lower formation energy than the BFrenkel and Bint models. However, all three models have relatively large (>2 eV) positive formation energies, and the defects will thus be present at very low concentrations. With increasing reducing conditions (more negative µ′O) the BFrenkel and Bint models progressively increase their stability relative to the Bvac model. The plots of Ed vs µ′O for the Bint and BFrenkel models have the same slope, but the Bint model is more stable than the BFrenkel model by 0.6 eV, based on the 2 × 2 × 2 supercell calculations. For these two models, we tested the effect of increasing the supercell size on the calculated total energy. For the simple interstitial model, increasing the size to 3 × 3 × 3 supercell only changed the value of Ed by 0.1 eV. However, the more extended Frenkel defect model was sensitive to the supercell size. Increasing to a 2 × 2 × 3 supercell gave an Ed value that was lower (more stable) by 0.28 eV. A further increase to a 3 × 3 × 3 supercell gave a further small stabilization by 0.07 eV. The plot of Ed vs µ′O for the BFrenkel model in the 3 × 3 × 3 supercell is included in Figure 5. When the dopant-dopant interactions are minimized by using the 3 × 3 × 3 supercells, the difference in defect formation energies for the BFrenkel and Bint models is only 0.27 eV. At a µ′O value -2.23 eV, corresponding to the lower oxygen partial pressure limit for B-doped rutile established in the B-doping experiments,19 the defect formation energies for the Bint and BFrenkel models are 0.53 and 0.8 eV, respectively. These values can be compared with a value of 0.66 eV calculated using the experimentally determined maximum B content in rutile at 1500 K of 2.4 atom %, in the expression c ) Nsitese-Ed/RT

(13)

where c is the concentration of dopant and Nsites is the number of sites available to the dopant. In the case of the BFrenkel model, a more accurate estimate of the defect formation energy at elevated temperatures was obtained by calculating the Helmholtz vibrational and configurational free energy contributions. The results for Fvib and Fconf for pure rutile and for the BFrenkel model are given in Table 3 at temperatures of 700 and 1500 K, relevant to sol-gel doping and high-temperature doping, respectively. The individual vibrational free energies are very large (-21 to -22 eV at 1500 K), but there is considerable cancellation between the contributions from doped rutile and pure rutile in the term Fvib(defect) - Fvib(host). The B-doping involves the formation of short B-O bonds with associated higher vibrational frequencies than occur in pure rutile, and this gives a more positive Fvib value for B-doped rutile compared with that for pure rutile. The

10920 J. Phys. Chem. C, Vol. 111, No. 29, 2007

Wilson et al.

TABLE 2: QM-Relaxed Cell Parameters (Å) and Calculated Bond Lengths (Å) and Angles (°) Associated with Different B-Doped Rutile Modelsa Bint model

a

BFrenkel model

a ) 9.292 b ) 9.312 c ) 5.958 B-O(1) B-O(2) B-O(3) B-O(4)

1.49 1.50 1.49 1.49

O(1)-B-O(2) O(1)-B-O(3) O(1)-B-O(4) O(2)-B-O(3) O(2)-B-O(4) O(3)-B-O(4)

123.8 102.8 102.9 102.8 102.9 123.7

Ti(1)-O(1) Ti(1)-O(4)

2.12 2.19 2.00

Ti(2)-O(4)

2.09 1.99

Ti(3)-O(2) Ti(3)-O(4)

2.19 2.12 2.00

Ti(4)-O(2)

2.09 1.99

Ti(5)-O(2) Ti(5)-O(3)

2.12 2.19 2.00

Ti(6)-O(3)

2.08 2.00

Ti(7)-O(1) Ti(7)-O(3)

2.19 2.13 2.00

Ti(8)-O(1)

2.09 1.99

a ) 9.231 b ) 9.455 c ) 5.898 B-O(1) B-O(2) B-O(3)

Bvac model

1.39 1.40 1.40

a ) 9.310 b ) 9.238 c ) 5.913 B-O(1) B-O(2) B-O(3)

1.39 1.41 1.39

O(1)-B-O(2) O(1)-B-O(3) O(2)-B-O(3)

118.0 118.0 123.9

O(1)-B-O(2) O(1)-B-O(3) O(2)-B-O(3)

117.8 119.4 122.9

Ti(1)-O(4) Ti(1)-O(5) Ti(1)-O(6) Ti(1)-O(7) Ti(1)-O(8) Ti(1)-O(9)

1.84 1.88 1.88 1.96 2.23 2.23 2.00

H-O(4)

0.99

Ti(2)-O(7) Ti(2)-O(8) Ti(2)-O(9) Ti(2)-O(10) Ti(2)-O(11) Ti(2)-O(12)

2.13 2.02 2.02 1.90 1.89 1.89 1.975

Ti(1)-O(2) Ti(1)-O(3) Ti(1)-O(4) Ti(1-O(5) Ti(1)-O(6) Ti(1)-O(7)

2.09 2.35 2.00 1.75 1.98 1.97 2.02

Ti(2)-O(4)

2.00 1.96

Ti(3)-O(3)

2.00 1.97

Ti(1)-Ti(2)

2.85 (face)

Ti(4)-O(1)

2.07 1.97

Ti(3)-O(2)

2.04 1.985

Ti(5)-O(1)

2.07 1.975

Ti(4)-O(3)

2.04 1.99

Ti(6)-O(2)

2,06 1.97

Ti(5)-O(1)

2.05 1.985

Ti(6)-O(1)

2.05 1.985

Ti(7)-O(2) Ti(7)-O(3)

2.19 2.19 2.05

See Figures 2-4 for atom labeling.

resulting Fvib(defect) - Fvib(host) term is thus positive whereas the Fconf(defect) - Fconf(host) term is negative and of a similar magnitude, so a further cancellation of values occurs. As shown in Table 3, the resultant change to Ed is +0.32 eV at 700 K, rising to +0.57 eV at 1500 K.

TABLE 3: Calculated Configurational and Vibrational Free Energies (eV) for Pure Rutile and B-Doped Rutile (BFrenkel Model)a Fvib for BFrenkel model Fvib for pure rutile Fconf for BFrenkel model Fconf for pure rutile net contribution to Ed ) Fvib(defect) Fvib(host) + Fconf(defect) -Fconfhost) a

Figure 5. Plot of the defect formation energy, Ed, as a function of reduction potential, µ′O, for boron incorporation into vacant octahedral sites in rutile at 1500 K.

700 K

1500 K

-2.95 -3.52 -0.25 0 0.32

-21.04 -22.15 -0.54 0 0.57

Results apply to 2 × 2 × 2 supercells.

Detailed vibrational energy calculations were not conducted for the Bint model. However, based on the results obtained for the BFrenkel model, some qualitative comments can be made. In the Bint model, the B is 4-coordinated, with longer B-O bonds compared with those for 3-coordinated B in the BFrenkel model (see Table 2). The lower B-O frequencies in the phonon DOS will result in the Fvib(defect) being more negative (more stable) for the Bint model compared with the BFrenkel model. The configurational entropy terms are the same for both models, so greater cancellation of Fvib and Fconf terms with a resulting

Nonmetal Doping in Rutile

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10921

TABLE 4: Oxygen-Rich Defect Formation Energiesa defect model

Ed (eV) B-doping

doping at Frenkel defect 5.19 (2×2×2 cell) 4.91 (2×2×3 cell) 4.84 (3×3×3 cell) doping at interstitial site 5.67 (3-coord. B) 4.58 (4-coord. B) 4.57 (4-coord. B, 3×3×3) doping at VTi lattice site 2.22 (2×2×2 cell) 2.12 (2×2×3 cell) 1.16 addition of 1H at VTi 2.60 addition of 2H at VTi

Ed (eV) C-doping 7.94 (2×2×2 cell) 7.68 (2×2×3 cell) 9.61 (1-coord. C) 8.57 (4-coord. C) 2.65 (2×2×2 cell) 2.69 (2×2×3 cell) 3.81 5.40

a Sol-gel conditions (700 K) results for 2 × 2 × 2 supercell unless specified.

smaller contribution to the defect formation energy might be expected for the Bint model. A smaller destabilisation of Ed for the Bint model would ensure that the relative stabilities of the BFrenkel and Bint models would remain as shown in Figure 5, that is, with the Bint model slightly more stable. More generally, such cancellation of vibrational and configurational free energies is expected for small, light atom dopants, which form short D-O bonds. The high-frequency contributions to the phonon DOS associated with the dopant atom and the anions to which it coordinates, destabilise the doped phase relative to the pure host phase (Fvib(defect) - Fvib(host) is +ve), whereas the configurational entropy confers extra stability on the doped phase (Fconf(defect) - Fconf(host) is -ve). The corollary of this observation is that a structure modification that results in lengthened, weaker bonds will have a more negative vibrational free energy. In certain circumstances, this term may be dominant enough to stabilize the phase. Such vibrational energy stabilization has been proposed to explain the negative slope of phase boundaries in p-T diagrams for certain phase transitions where the higher-pressure phase has a higher coordination number (longer bonds) for the metal atoms.32 4.2.3. 700 K, Oxygen-Rich Results. Defect formation energies at 700 K at the oxygen-rich limit relevant to sol-gel conditions are given in Table 4. For the Dvac and Dint models, in which the defects are localized to single sites, increasing the size of the supercell above 2 × 2 × 2 had only a small effect on Ed (e0.1 eV). However, for the DFrenkel model, where the defect is extended over a lattice site and the adjacent interstitial site, increasing the supercell from 2 × 2 × 2 to 2 × 2 × 3 gave a decrease in the calculated Ed of 0.26 to 0.28 eV, indicating that defect-defect interactions are causing some destabilisation in the 2 × 2 × 2 supercell. A further increase in the supercell size to 3 × 3 × 3 gave only a small further decrease in Ed (by 0.07 eV). In the case of the Bint model, relaxations were conducted for two different boron coordinations. The 4-coordinated B case is more stable by 1.1 eV than for the model in which the B remained 3-coordinated during the relaxation. As previously noted for the 1500 K results, the Bvac model is considerably more stable than the other two models under oxygen-rich conditions. However, the calculated defect formation energies are still too high, at 2.1-2.2 eV, for B to be a significant dopant in sol-gel prepared rutile. Relaxations were conducted for the Bvac model with the addition of protons at the VTi site. The addition of one proton at the site, corresponding to the charge-balanced substitution of B3+ + H+ for Ti4+, conferred an extra stability of 1 eV, compared to pure B doping. Addition of a second proton, however, caused de-stabilization, with Ed being higher than for pure B doping.

Figure 6. Plot of the defect formation energy, Ed, as a function of reduction potential, µ′O, for carbon incorporation into vacant octahedral sites in rutile (2 × 2 × 2 supercell) at 1500 K.

4.3. Results for C-Doping. The variation in defect formation energy as a function of µ′O at 1500 K for the three different models for C-doping at vacant octahedral sites in rutile is shown in Figure 6. At the oxygen-rich limit, all three models have high Ed values, so C-dopant concentrations will be vanishingly small. The Ed values for the Cint and CFrenkel models decrease with increasing reducing conditions, but even at the oxygenpoor limit where solid C-deposition occurs, the defect formation energies for the three models all remain above +4 eV. A comparison of Figures 5 and 6 shows that the relative stabilities of the Dint and DFrenkel defects invert when D ) C compared with D ) B. Di Valentin et al.12 also applied DFT modeling to C-doping of rutile and reported the variation of Ed as a function of oxygen potential for the Cint and Cvac cases. However, there are some significant differences between the methodologies used by Di Valentin et al.12 and in our approach. In particular, Di Valentin et al.12 did not include the influence of temperature and pressure on the chemical potentials of Ti and C. As seen from Table 1, the change to µTi is 2 eV from 0 to 1500 K, and this directly affects the relative stabilities of the Cint and Cvac models. Second, Di Valentin et al.12 assumed that CO2 was an independent chemical potential reservoir for C, not in equilibrium with oxygen. This results in a different dependence of Ed on µ′O for the dopant models to those shown in Figure 6. A further difference is the definition of the oxygen-rich reference point for µ′O. Di Valentin et al.12 followed the procedure used by Reuter and Scheffler,25 and used µ′O ) µO - 1/2EO2. This results in a very large value of the oxygen pressure (∼1010 bars) at µ′O ) 0. In contrast, we equated µ′O to 1/2RT ln pO2, so µ′O ) 0 then corresponds to an oxygen pressure of 1 bar. When these differences are taken into account, reasonable agreement between the relative stabilities of the Cint and Cvac models at the oxygen-rich limit is found between the two studies, with the Cvac model being more stable by ∼5 eV. The Ed values obtained at the oxygen-rich limit at 700 K for the three C-doping models are reported in Table 4. The least stable defect involves location of the carbon at an interstitial site. Di Valentin et al.12 reported that their relaxed interstitial model had the carbon coordinated to only one oxygen atom, with C-O ) 1.287 Å. We were able to reproduce the singly coordinated defect, with a C-O distance of 1.253 Å in the relaxed structure. However, as seen from Table 4, we found that an alternative Cint model with 4-coordinated carbon (C-O distances of 1.413-1.414 Å) was 1.04 eV more stable than the singly coordinated carbon model. This is an example of a recurring problem with ab initio defect studies; the relaxed

10922 J. Phys. Chem. C, Vol. 111, No. 29, 2007 structure and the resulting total energy depend sensitively on the starting model and the way in which the relaxation is conducted. For the Cvac model, additional relaxations were performed for cases involving the incorporation of protons at the VTi site. In contrast to the case of B-doping, the addition of one proton at the C-doped VTi site resulted in a higher defect formation energy. The relatively high Ed values for all models reported in Table 4 suggests that carbon doping of rutile under sol-gel conditions is unlikely to occur via incorporation of the dopant at vacant octahedral sites in the structure, at an annealing temperature of 700 K. 5. Conclusions The present study confirms that high-temperature nonmetal doping of rutile can be reliably modeled using 0 K ab initio total energy calculations combined with corrections of the chemical potential reservoirs for free energy changes with temperature and pressure. The latter can be obtained from published tables of thermodynamic properties. Neglect of corrections of the total energies of the host and doped phases for vibrational and configurational free energy changes with temperature was found to introduce only modest errors, which did not change the relative stabilities of different dopant models. This was confirmed quantitatively by ab initio determinations of the phonon DOS for pure rutile and for one of the doped models and application of the results to calculate the vibrational free energies. It was found that the vibrational and configurational free energy differences between host and defect phases were of similar magnitude and of opposite sign, resulting in extensive cancellation, so that the net change to the defect formation energy due to these contributions was of the order of 0.3 eV at 700 K. The correction of the chemical potential reservoirs for free energy changes as a function of temperature has been shown to be important because such corrections can change the relative stabilities of alternative models, and of different dopants. For example, under oxygen-rich conditions, the Ed values for C and B doping at a VTi site (Dint model) at 0 K are 1.45 and 2.02 eV, respectively, and thus C-doping is more stable by 0.57 eV. However at 700 K, the corresponding Ed values for C and B doping are 2.65 and 2.22 eV. The stabilities are thus reversed at 700 K with B-doping more stable by 0.43 eV. The QM calculations confirm that incorporation of boron into rutile is energetically favorable at elevated temperatures under strongly reducing conditions, consistent with previous experimental boron-doping experiments.19 We had originally proposed that the boron is incorporated at Frenkel defects in rutile and has triangular coordination as in TiBO3.27 However, the QM modeling shows that incorporation of B at vacant octahedral (interstitial) sites within the [001] channels is more likely, as the defect formation energy is 0.27 eV more stable at 1500 K for the Bint model than for the BFrenkel model. In contrast to the results for B-doping, the QM calculations show that incorporation of C into rutile under reducing conditions at elevated temperatures will be negligible at any of the three different octahedral vacancy sites.

Wilson et al. Under oxygen-rich sol-gel conditions, the most favorable octahedral vacancy doping conditions are for B at a protonated VTi site. The defect formation energy is 1.16 eV at 700 K. At 0 K the (B+H)vac model is stabilized by a further 0.51 eV (Ed ) 0.65 eV). The 0 K defect formation energy for C-doping at a VTi site is 1.45 eV, rising to 2.65 eV at 700 K. These results suggest that annealing temperatures should be kept as low as possible to retain boron and carbon dopants in the rutile structure under oxygen-rich conditions. Acknowledgment. This work was supported by the Australian Partnership for Advanced Computing (APAC) and the Victorian Partnership for Advanced Computing (VPAC) supercomputer facilities. We thank the CSIRO Energy Transformed Flagship (ETF) for support of the project in which this study was conducted and for providing one of us (S.P.R.) with an ETF visiting scientist fellowship. Thanks to Graham Sparrow and Shoui Sun for internal refereeing of the manuscript. References and Notes (1) Nowotny, M. K.; Bak, T.; Nowotny, J. J. Phys. Chem. B 2006, 110, 16302. (2) Mizushima, K.; Tanaka, M.; Iida, S.; Goodenough, J. B. J. Phys. Chem. Solids 1979, 40, 1129. (3) Nowotny, M. K.; Bak, T.; Nowotny, J.; Sorrell, C. C. Physica Status Solidi B 2005, 242, R88. (4) Na-Phattalung, S.; Smith, M. F.; Kim, K.; Du, M. H.; Wei, S. H.; Zhang, S. B.; Limpijumnong, S. Phys. ReV. B 2006, 73. (5) Bokhimi, A.; Morales, A.; Novaro, O.; Lopez, T.; Sanchez, E.; Gomez, R. J. Mater. Res. 1995, 10, 2788. (6) Grey, I.; Madsen, I.; Bordet, P.; Wilson, N.; Li, C. AdV. Ecomater. 2005, 1, 35. (7) Grey, I. E.; Wilson, N. C. J. Solid State Chem. 2007, 180, 670. (8) Wang, J. A.; Limas-Ballesteros, R.; Lopez, T.; Moreno, A.; Gomez, R.; Novaro, O.; Bokhimi, X. J. Phys. Chem. B 2001, 105, 9692. (9) Ruetschi, P. J. Electrochem. Soc. 1984, 131, 2737. (10) Asahi, R.; Morikawa, T.; Ohwaki, T.; Aoki, K.; Taga, Y. Science 2001, 293, 269. (11) Di Valentin, C.; Pacchioni, G.; Selloni, A. Phys. ReV. B 2004, 70. (12) Di Valentin, C.; Pacchioni, G.; Selloni, A. Chem. Mater. 2005, 17, 6656. (13) Di Valentin, C.; Pacchioni, G.; Selloni, A.; Livraghi, S.; Giamello, E. J. Phys. Chem. B 2005, 109, 11414. (14) Geng, H.; Yin, S. W.; Yang, X.; Shuai, Z. G.; Liu, B. G. J. Phys.Condens. Matter 2005, 18, 87. (15) Lin, Z. S.; Orlov, A.; Lambert, R. M.; Payne, M. C. J. Phys. Chem. B 2005, 109, 20948. (16) Umebayashi, T.; Yamaki, T.; Itoh, H.; Asai, K. Appl. Phys. Lett. 2002, 81, 454. (17) Wang, H.; Lewis, J. P. J. Phys.-Condens. Matter 2005, 17, L209. (18) Zhao, W.; Ma, W. H.; Chen, C. C.; Zhao, J. C.; Shuai, Z. G. J. Am. Chem. Soc. 2004, 126, 4782. (19) Grey, I. E.; Li, C.; MacRae, C. M.; Bursill, L. A. J. Solid State Chem. 1996, 127, 240. (20) Kresse, G.; Furthmuller, J. Phys. ReV. B 1986, 54, 11169. (21) Pack, J. D.; Monkhorst, H. J. Phys. ReV. B. Condens. Matter 1977, 16, 1748. (22) Parlinski, K. PHONON 2002. (23) Sikora, R. J. Phys. Chem. Solids 2006, 66, 1069. (24) Barin, I. Thermochemical Data of Pure Substances; VCH: Weinheim, Germany, 1989. (25) Reuter, K.; Scheffler, M. Phys. ReV. B 2002, 65. (26) Bursill, L. A.; Blanchin, M. G. J. Phys. Lett. 1983, 44, L165. (27) Huber, M.; Deiseroth, H. J. Z. Kristallogr. 1995, 210, 685. (28) Harris, L. A.; Quong, A. A. Phys. ReV. Lett. 2004, 93. (29) Sullivan, J. M.; Erwin, S. C. Phys. ReV. B 2003, 67. (30) Thompson, S. J.; Lewis, S. P. Phys. ReV. B 2006, 73. (31) Baur, W. H. In Structure and Bonding in Crystals; O’Keeffe, M. a. N. A, Ed.; Academic Press: New York, 1981; Vol. II, p 31. (32) Navrotsky, A. J. Geophys. Res. Lett. 1980, 7, 709.