The problem of a metal impurity in an oxide: ab-initio study of electronic and structural properties of Cd in Rutile TiO2

In this work we undertake the problem of a transition metal impurity in an oxide. We present an ab-initio study of the relaxations introduced in TiO2 when a Cd impurity replaces substitutionally a Ti atom. Using the Full-Potential Linearized-Augmented-Plane-Wave method we obtain relaxed structures for different charge states of the impurity and computed the electric-field gradients (EFGs) at the Cd site. We find that EFGs, and also relaxations, are dependent on the charge state of the impurity. This dependence is very remarkable in the case of the EFG and is explained analyzing the electronic structure of the studied system. We predict fairly anisotropic relaxations for the nearest oxygen neighbors of the Cd impurity. The experimental confirmation of this prediction and a brief report of these calculations have recently been presented [P.R.L. 89, 55503 (2002)]. Our results for relaxations and EFGs are in clear contradiction with previous studies of this system that assumed isotropic relaxations and point out that no simple model is viable to describe relaxations and the EFG at Cd in TiO2 even approximately.


I. INTRODUCTION
The problem of metal impurities in oxides is a challenge in solid-state physics both from a fundamental and applied point of view. Description of impurity levels in oxide semiconductors has a great technological interest, but a complete theoretical description of the problem is very difficult. In effect, impurities introduce atomic relaxations in the host and modify the electronic structure of the system, being the interplay between these two effects one of the principal difficulties in the theoretical approach. Some experimental techniques used in material science introduce probe-atoms into the systems to be studied, 1,2 giving valuable physical information seen by the probe-atom, usually an impurity in the system. In particular, TDPAC spectroscopy 3 is a hyperfine technique that enables a high resolution determination of the EFG tensor at impurity sites. Due to the r −3 dependence of the EFG operator from the charge sources, the EFG is mostly originated in the electronic charge density close to the impurity nucleus, which in turn reflects the probe chemistry with its neighborhood. Thus, it would be important that the theoretical approach to the metal-impurity problem would be able to predict the EFG tensor with the required accuracy. Since the EFG is very sensitive to the anisotropic charge distribution close to the nucleus, for its accurate calculation the entire electronic configuration of the host, perturbed by the presence of the impurity, has to be determined. This kind of calculations should be performed at the ab initio level and can be done in the frame of the Density Functional Theory (DFT). A well proven method to solve the all-electron DFT equations in solids is the Full-Potential Linearized-Augmented-Plane-Wave (FLAPW) method. 4 This method was able to predict with high accuracy and precision the EFG at undisturbed lattice sites in a large variety of pure systems. 5,6 However, on systems with impurities, very few calculations have been performed and the state-of-the-art is far from being routinely in this field. 7,8,9,10,11,12,13,14,15 The first fully selfconsistent ab initio determination of the EFG tensor at an impurity site (Cd) in an oxide (TiO 2 ) has recently been reported 16 . In that work we performed a FLAPW calculation of the relaxations introduced by the impurity and studied their interplay with the electronic structure of the system, predicting highly anisotropic relaxations of the nearest neighbors of the impurity and a drastic change in the orientation of the principal component of the EFG tensor. This prediction was confirmed in the same work by a key PAC experiment. Due to the complexity of the physical situation we want to describe, even in the framework of ab initio calculations some variables have to be determined through comparison with experiments. For example, we performed the calculations for two charge states of the impurity and one of them could be discarded because it gave a value for the EFG asymmetry parameter, η, incompatible with experience.
The central purpose of this work is to present and discuss the electronic structure of the TiO 2 :Cd impurity system with more detail that we could do in Ref. 16. Here, we study the relationship between the calculated EFG at the impurity site -the key experimental quantity -and the interplaying physical quantities: the structural relaxations and the character and filling of the impurity state. We also discuss in detail the precision of our calculations studying its dependence on the different parameters that control the precision (impurity dilution, base dimension, average in k-space, etc.) and its dependence on the approximation used for the exchange-correlation potential.
Another important question to be addressed concerning all experimental techniques that introduce impuritytracers in solids is at what extent simple models usually used to predict the measured quantities are acceptable in view of the modifications introduced by the presence of the impurity. Concerning the EFG, simple approximations like the point-charge model 2 with antishielding factors 17,18 and isotropic relaxations of nearest-neighbor atoms 18,19 are commonly used. In this work we give arguments against the applicability of these simple models for the studied system. This paper is organized as follows. In section II we describe the method of calculation used in the present work and include some preliminary calculations for an unrelaxed cell in order to discuss methodologically the point of the charge state of the impurity. In section III.(A-C) we present and discuss the results of performing relaxations of the Cd nearest-neighbors (NN) in a 72-atoms super-cell for two different states of charge of the impurity, while in section III.D we study the accuracy of these results performing several additional calculations. Even if section III.D could be found rather technical by people not specially interested in first-principles calculations, we think it will be useful for the ab initio community interested in the problem of impurities in solids. In section III.E we compare our results with previous studies of this system. Finally, in section IV we present our conclusions.

A. A general scheme
Our aim is to calculate from first principles the structural relaxations produced in the TiO 2 lattice when a Cd impurity replaces a Ti atom and the electronic structure of the resulting system. In particular we are interested in the EFG tensor at the Cd site. The general scheme we adopted to deal with this problem was the following: we considered a super-cell (SC), containing a single Cd impurity, repeated periodically and performed firstprinciples calculations We studied the relaxation introduced by the impurity computing the forces on the Cd neighbors and moving them until the forces vanished. The calculations were performed with the wien97 implementation, developed by Blaha et al., 20 of the Full-Potential Linearized-Augmented-Plane-Wave (FLAPW) method and we worked in the LDA approximation using the Perdew and Wang parametrization for the exchangecorrelation potential. 21 In this method the unit cell is divided into non-overlapping spheres with radii R i and an interstitial region. The atomic spheres radii used for Cd, Ti and O were 1.11, 0.95 and 0.85Å, respectively. We took for the parameter RK MAX , which controls the size of the basis-set in these calculations, the value of 6 that gives basis-set consisting in 4500 LAPW functions for the SC described in section II B. We also introduced local orbitals (LO) to include Ti-3s and 3p, O-2s and Cd-4p orbitals. Integration in reciprocal space was performed using the tetrahedron method taking an adequate number of k-points in the first Brillouin zone. Once self-consistency of the potential was achieved, quantummechanical-derived forces were obtained according to Yu et al. 22 , the ions were displaced according to a Newton damped scheme 23 and new positions for Cd neighbors were obtained. The procedure is repeated until the forces on the ions are below a tolerance value taken as 0.025 eV/Å. At the relaxed structure, the V ij elements of the EFG tensor are obtained directly from the V 2M components of the lattice harmonic expansion of the selfconsistent potential (a more detailed description of the formalism used to compute the EFG tensor may be found in Ref. 6).
There is still an important point to take into account concerning the calculation of the electronic structure of an impurity system: the charge state of the impurity. We let this point to the end of this section as it will be easily handled after the analysis of some preliminary calculations. The rutile (TiO 2 ) structure is tetragonal (a = b = 4.5845 1Å , c = 2.9533 1Å ). The unit cell (shown in Fig. 1) contains 2 metal atoms (Ti) at positions 2a (0, 0, 0) and (1/2, 1/2, 1/2) and 4 anions (O) at positions 4f ±(u, u, 0; u+1/2, 1/2-u, 1/2), with u = 0.30493 7 . 24 The supercell considered in the present work consists of twelve unit cells of TiO 2 where one Ti atom has been replaced by a Cd atom. The resulting 72-atom SC (called 72A-SC in the future) has dimensions a ′ =2a, b ′ =2b, c ′ =3c and is also tetragonal with c ′ /a ′ =0.97 giving an almost cubic lattice. This SC keeps Cd atoms as far as possible from each other for the given cell volume. For checking purposes we have also considered SCs containing eight and sixteen unit cells. We assume that relaxations to be performed preserve point group symmetry of the cell in this initial configuration. Symmetry restricts O1 and O2 displacements to yz plane and x direction respectively. In order to check this assumption and the stability of the obtained solution, at the end of the relaxation process we performed new calculations with O1 and O2 atoms displaced from their symmetry positions and verified that these solutions are a minimum for each system studied.
The self-consistent calculations were performed taking 8 k-points in the irreducible Brillouin zone for the metallic system (situation (i)), and 2 k-points for the nonmetallic ones (situation (ii)). In order to plot the DOS we calculate eigenvalues at a denser mesh of 36 k-points.

C. Preliminary study: Charge state of the impurity
In order to discuss the problem of the charge state of the impurity we have done some preliminary calculations. We have calculated the self-consistent electronic structure of the 72A-SC with all atoms in their initial unrelaxed positions. We want to analyze the changes in the electronic structure of the system caused by the presence of the impurity, neglecting for the moment the problem of structural relaxations. The density of states of this system is compared with the one of pure TiO 2 in Fig. 2.
is a wide band-gap semiconductor with the O-p band filled and the Ti-d band empty. Since Cd valence is 2+, when a Cd atom replaces a Ti atom in the SC the resulting system is metallic because of the lack of two electrons necessary to fill up the O-p band. The question that arises at this point is if the real system we want to describe provides the lacking two electrons (via an oxygen vacancy, for example) or not, and if this point is relevant for the present calculation. We will show first that this point is absolutely relevant for the present calculation. Comparison of Fig. 2(a) and (b) shows that the presence of the Cd impurity in the SC introduces the appearance of Cd-d levels and impurity states at the top and the bottom of the valence band in the corresponding DOS. The ones at the top of the valence band can be better seen looking at the band structure in this energy range shown in Fig. 3. The two bands that are immediately above and below the Fermi level when cross- ing A-point correspond to impurity anti-bonding states that are spatially located at Cd and at their O1 and O2 nearest-neighbor atoms. In particular, the wave function of the impurity state that remains almost completely unoccupied has character Cd-d yz , O1-p y , and O1-p z . Then, providing two electrons to the system implies a drastic change in the symmetry of the electronic charge distribution in the neighborhood of the impurity. Therefore, different nearest-neighbors relaxations and also different EFG could be expected for different charge states of the impurity.
In the present work, we present calculations for two charge states of the impurity system, corresponding to two different physical situations: (i) Cd 0 (neutral impurity state), corresponding, e.g., to an extremely pure crystal at low temperatures, and (ii) Cd − (charged acceptor state), the system provides two electrons via an oxygen vacancy, donor defects, etc. To study case (i) we used the 72A-SC already described, to study case (ii) we also used the 72A-SC but performed self-consistent calculations adding two electrons to the SC that we compensate with an homogeneous positive background in order to have a neutral cell to compute total energy and forces (this procedure is implemented since version wien97.9 of the FLAPW package). We have also simulated situation (ii) with an alternative procedure: we replaced in the 72A-SC the two most distant oxygen atoms by two fluorine atoms. In this way we provide two electrons to fill up the O-p band without introducing any artificial background. We expected that the difference between flourine and oxygen potentials should modify only slightly the re- sults. In summary, we performed self-consistent FLAPW calculations for the following systems:

A. Structural relaxations
Let us first consider the relaxation of only the 6 nearest oxygen neighbors of the Cd impurity (O1 and O2 in Fig 1.) In Table I we compared the results of the relaxation of these oxygen atoms for the different systems studied. We see that for both charge-states of the impurity the relaxations are quite anisotropic, with the Cd-O1 distance larger than Cd-O2 distance, opposite to the initial unrelaxed structure. This result is opposite to what other authors have assumed in previous studies of this system 18,19 and confirms the tendency predicted in our previous calculation with a much smaller SC. 25 As it can be seen in Table I the difference in the charge state of the impurity affects essentially the relaxation of O1 atoms that present a slightly larger relaxation for the charged impurity. Relaxation for the two ways of simulating the charged state of the impurity (case Cd − ) are very similar indicating that both approaches are well suited to deal  a considerable shortening in Ti-O2 bonds. However, stretching of Cd-O1 bond doesn't affect so much Ti-O1 bonds since the structure is more open in this direction. So, at the end of the relaxation, Cd-O1 bond stretches almost twice than Cd-O2 bond.

B. Electronic structure
In Fig. 5 we show the bands for the 72A-SC for the two states of charge considered for the impurity. Comparison of Fig. 5(a) with Fig. 3 shows that as a consequence of relaxation the outermost impurity level falls down in energy and goes into the O-p band becoming half-occupied. Relaxation increases the Cd-O1 distance and this produces a softening of the Cd-d−O1-p interaction and the fall down in energy of the anti-bonding impurity states. Figure 5 shows that the band structure of the neutral and charged relaxed structures are very similar, but the outermost impurity state is a slightly raised in energy when it is completely filled (case (ii)) as a consequence of the larger Coulomb repulsion. Comparison of Fig. 6(a) and To look at the orbital composition of the impurity states we plot in Fig. 7(a) and (b) the PDOS for Cd-d, O1p and O2-p for the two charge states of the impurity. The impurity state near the top of the valence band has Cd-d yz , O1-p y and O1-p z character and, as we mention looking at the bands, it is shifted upwards in energy when it is completely filled (case Cd − ). In Fig. 7(a)(Cd atom) and (b)(Cd atom) it can also be seen the presence of other impurity state with components Cd-d x 2 −y 2 and Cd-d 3z 2 −r 2 which is completely filled in both cases. This impurity state is located mainly at Cd-d x 2 −y 2 and O2-p x but involves also contributions from O1-p z , Cdd 3z 2 −r 2 and O3-p x (O3 is the NN of O2 atom in the x direction). From the present calculations the outermost impurity state has an occupation of around 1.3 e for the 72A-SC(Cd 0 ) and 2 e for the 72A-SC(Cd − ). We identify this 0.7 additional electrons in this state as the driving force that produce the slightly larger relaxation for O1 atoms in the charged cell with respect to the neutral cell. We want to mention that the fact that the impurity level falls at the Fermi energy in the Cd 0 case is not fortuitous. Due to Coulomb repulsion, the impurity level would fall below the Fermi energy if it was empty an above if it was filled being the only self-consistent solution to be half-occupied at exactly the Fermi energy. The occupation of the impurity level in the 72A-SC(Cd 0 ) is therefore a constant number and fairly independent of small fluctuations of charge in the cell.

C. Electric-Field Gradients
In Table II we show the results for the V ii principal components of the EFG tensor for the three systems studied. The resulting EFGs for the two approaches used to simulate the charged impurity are very similar, the difference of 0.4 x 10 21 V/m 2 in components V XX and V Y Y is within one could expect for the small difference found in the oxygen positions, since EFG is very sensitive to small structural changes. These results agree very well (in magnitude, symmetry and orientation, see Table II) with the experimental results obtained for the EFG at Cd impurities substitutionally located at cationic sites in rutile TiO 2 . 16 The difference between the EFGs obtained for the charged and neutral cells is very remarkable: the sign, direction and absolute value of the largest V ii component (V 33 ) are different in both situations, and also the value of the asymmetry parameter η. The high η value obtained for 72A-SC(0) shows that the electron availability present in the sample leads the impurity to be in a charged state.
In order to investigate the origin of the difference in the EFG for the two charge states of the impurity we concentrate in the valence contribution to the EFG which originates in the asymmetry of the valence charge distribution inside the muffin-tin sphere. The valence contribution is usually dominant in FLAPW calculations and can be split in the different orbital symmetries. 6 In Table III we show the total valence contribution to V ii and its components arising from p and d orbital symmetries. We see that the largest differences correspond to    Table III shows that the changes δV ii are quite well described by this estimation giving for I d a value around 83 × 10 21 V/m 2 . Another interesting point is the presence of d contributions to V ii in the 72A-SC[Cd − (2e)] although the Cd-d band is completely filled (see Table III). This curious point may be explained because radial d wave functions are energy dependent, so, electrons in anti-bonding impurity states (at the Fermi level) contribute to V ii with a different 1/r 3 factor than the ones of the same symmetry in Cd-d levels (between 6 and 7 eV below the Fermi level). We have found an increment of 30% in 1/r 3 when going from an energy of -7 eV to the Fermi level. This dependence of 1/r 3 with energy was also mentioned and verified by Blaha et al. The main sources of error of the present study in order to compare with experiment are the size of the SC considered, the size of the basis, the LDA approximation used for the exchange-correlation potential and the relaxation process that has been restricted to the six nearest oxygen neighbors of the Cd atom. To check the accuracy of the present study we have performed several additional calculations.    b. Size effects In order to check how appropriate are the dimensions of the 72A-SC (2a × 2b × 3c) used in the present work we have performed self-consistent electronic structure calculations for a 48A-SC (2a × 2b × 2c) and a 96A-SC (2a × 2b × 4c). In these calculations we put the nearest oxygen neighbors of Cd at relative positions from Cd that correspond to equilibrium in the 72A-SC. The size of the forces that oxygen atoms O1 and O2 experiment in the 48A-SC and 96A-SC is a measure of the convergence of the present calculation. We obtained that in the 96A-SC the forces on O1 and O2 atoms point outwards (respect to Cd) and are around 0.15 eV/Å. Forces of this size produce, during the relaxation process in the 72A-SC, changes in distances smaller than 0.01Å and changes in EFG of about 0.2 × 10 21 V/m 2 , so, similar changes would be expected if relaxation in the 96A-SC was performed. In the case of the 48A-SC forces point inwards and are of the same magnitude for O2 atoms but they are about 0.65 eV/Å in O1 atoms. We see that there is a size effect that makes relaxations to be larger for larger SCs, but for the 72A-SC the effect is quite small and no significant variations should be expected if relaxations in bigger SCs were considered.
c. Basis size We compute the self-consistent electronic structure of the 72A-SC[Cd − (2e)] for the relaxed positions of Table I increasing the size of the basis to 7100 LAPW functions (RK MAX =7). Forces on O1 and O2 atoms are below the tolerance value of 0.025 eV/Å indicating that the result of the relaxation performed is unaltered for a substantial increase of the basis size. In fact, the forces on all the other atoms in the cell have the same values than using RK MAX =6 (within the tolerance) except for the Ti neighbors of Cd in the z direction where the forces differ in 0.12 eV/Å. Changes obtained in V ii components of the EFG tensor are smaller than 0.05 × 10 21 V/m 2 . We have also considered the inclusion of a Cd-4d LO to improve linearization. When a Cd-4d LO is introduced with energy at the Fermi level in order to improve the description of the impurity states no influence in the forces is detected. The change in V ii components of the EFG is smaller than 0.1 × 10 21 V/m 2 .
d. Exchange-correlation potential We performed electronic self-consistent calculations and relaxation of Cd nearest oxygen neighbors for the system 72A-SC[Cd − (2e)] using the generalized gradient approximation (GGA) 30 instead of LDA. With the use of this parametrization for the exchange-correlation potential we obtained 2.18Å and 2.10Å for Cd-O1 and Cd-O2 distances, respectively (i.e. only a small change of 0.01Å in Cd-O2 distance). For the EFG we obtained V 33 =V Y Y =+4.94 × 10 21 V/m 2 and η=0.37.
e. Further relaxations In order to study the effect of relaxing the coordinates of atoms beyond the nearest neighbors of the Cd atom we have performed the following two additional relaxations: (a) we allow to relax the coordinates of all atoms within a cutoff radius R C =4Å centered at Cd (this involves 24 atoms) and (b) idem for R C =4.6Å (this involves 42 atoms). For radius larger than 4.6Å, atoms to be relaxed would be nearer to the images of Cd atom from neighboring cells than to the Cd itself, so we consider this radius as a limit for the present SC. In Table IV we compare the results obtained for both relaxations with the ones corresponding to the NN relaxation. We observe that there is not any qualitative change in the results already discuss for the NN relaxation. There are, however, some small variations in the values predicted for the EFG. The differences between results from relaxations (a) and (b) (that are as larger as the differences between relaxations (a) and NN) are in part caused because in (b) the relaxation of O3 atom is allowed. Atom O3 is the NN of O2 atom in the X direction and its relaxation of about 0.04Å allows a further relaxation of about 0.012Å of the O2 atom. The fact that atom O3 is the one that experiments the largest relaxation (not considering O1 and O2) shows that directional bonding plays an important role in this structure.
In summary, we observed that none of the factors considered influence qualitatively the results, but their effect is neither negligible. We therefore confirm our prediction about sign (positive) and direction (Y ) of V 33 , but it is not possible to perform an exact prediction about its magnitude. Ours checks shows that d(Cd-O1) and d(Cd-O2) are converged within 0.01Å and V 33 and η within 0.5 × 10 21 V/m 2 and 0.1, respectively. Taking the values from the last row of Table IV as our predicted values for V ii we obtain discrepancies with experiment: 16,26 of 0.37 × 10 21 V/m 2 for |V 33 | and 0.12 for η that could be attributed to precision errors of the present calculations.

Neutral impurity: Cd 0
Let us briefly discuss the case of the 72A-SC(Cd 0 ). Due to the factor four that exists because of k-sampling, calculations are much more time consuming in this case. Hence, we have checked the size of the SC performing calculations only on a 48A-SC, and also checked the size of the basis repeating the calculation for the 48A-SC with RK MAX =7. Comparison of the values obtained for EFG and forces let us conclude that errors in the 72A-SC(Cd 0 ) are expected to be of the same magnitude than the ones obtained for the 72A-SC[Cd − (2e)]. Only the effect of adding a Cd-4d LO and the use of GGA instead LDA produce a larger variation of the EFG in this case than in the charged cell, and is understandable because of the larger d contribution to V ii in the case of the neutral state of the impurity. Calculations in the 72A-SC(Cd 0 ) with this parametrization for the exchange-correlation potential and introducing Cd-4d LO predict values of V 33 =V XX =-7.90× 10 21 V/m 2 and η=0.80. In addition, we found a force of 0.04 eV/Å inwards at O1 atoms, so a small refinement in O1 position would be expected. About further relaxations, we have only performed relaxation (a) (R C =4Å) for the 72A-SC(Cd 0 ) and also obtain variations of the same order than in 72A-SC[Cd − (2e)], in particular we obtained V 33 =V XX =-6.8× 10 21 V/m 2 and η=0.97. Then, the same conclusions about accuracy than in the 72A-SC[Cd − (2e)] hold for the 72A-SC(Cd 0 ) but note that, due to the high value of η in this case, the sign and direction of eq could change because of precision.
Step of relaxation

E. Comparison with other calculations
The simplest and most widely used approximation for the calculation of the EFG at a probe-atom is the pointcharge model (PCM). 2 In this approximation the EFG tensor at the probe site is ( ij is the EFG tensor produced by valence nominal charges located at the ion positions in the lattice, and γ ∞ is the Sternheimer antishielding factor 17 that depends only on the probe-atom. In this way, this model assumes that the symmetry and orientation of the EFG tensor at impurity sites are unaltered by its presence. PCM gives for TiO 2 (Cd): V 33 =V XX =-2.27 × 10 21 V/m 2 and η=0.40 when a value of -29.27 is used for γ ∞ . 31 In Ref. 18 the authors assume that relaxation is responsible for the disagreement of PCM predictions with experiment and speculate that an isotropic relaxation of 0.04Å outwards of all the Cd nearest oxygen neighbors would produce the desired result for |V 33 | and η. Our results indicate that relaxations are not isotropic and they are so much larger than this, but even if the relaxed coordinates from our calculation (see Table I) were used, the PCM would fail in the description of the EFG giving V 33 =V XX = -8.20 × 10 21 V/m 2 in clear contradiction with V 33 =V Y Y = +4.55 × 10 21 V/m 2 that we obtain from the self-consistent FLAPW calculation. The disagreement of PCM with the FLAPW prediction could not either be attributed in this case to a change in the value of γ ∞ since the sign and directions of both predictions are different. Thus, it is clear that the problem of EFG at cationic sites in TiO 2 is too complicated to be described even approximately by simple PCM calculations.
In our previous work 25 we performed self-consistent electronic structure calculations of this system with a twelve-atom super-cell (12A-SC: a × b × 2c). Relaxations were performed only for the neutral charge state of the impurity and at the end of the relaxation process two electrons were added to compute EFG without computing the self-consistent potential of the charged cell. The structural relaxation obtained by this proce-dure were smaller than in the present work but account for the inversion of Cd-O1 and Cd-O2 distances [d(Cd-O1)=2.12Å, d(Cd-O2)=2.07Å] with respect to the unrelaxed structure (see Table I). The description obtained for the EFG was very similar to the one of the 72A-SC[Cd − (2e)], but we understand that this agreement is somehow fortuitous because usage of 12A-SC relaxed coordinates to compute the EFG in the 72A-SC[Cd − (2e)] gives very different results.
In Ref. 19, the authors have used a very similar approach as that we used in Ref. 25 using a 12A-SC, but they assumed that the relaxations of the nearest oxygen neighbors of the Cd atom were isotropic. As a consequence, they obtained d(Cd-O1)=2.04Å, d(Cd-O2)=2.08Å, and also a very different result for EFG (see Table I), V 33 =V ZZ =-5.09 × 10 21 V/m 2 and η=0.39 (for a carefully comparison of these two calculations see Ref. 25).
A question that arises at this point is if self-consistent electronic FLAPW calculations performed with a converged SC give EFG compatible with experience or not in the case that relaxations of the oxygen NN of the Cd atom are constrained to be isotropic. This is an interesting point in order to know if available experimental data is enough to refute the assumption that the oxygen NN relax isotropically. We have performed self-consistent calculations for the system 72A-SC[Cd − (2e)] for different positions of O1 and O2 atoms but moving them outwards, keeping the relation d(Cd-O1)/d(Cd-O2) constant. We determined the equilibrium position of oxygen atoms as the one that produce a minimum in the energy [see Fig. 8(a)]. We obtain that Cd-O distances are 2.12Å and 2.16Å for O1 and O2 atoms, respectively [a relaxation of 9% of the unrelaxed distances, see Fig. 8(a)]. If we compare Figs. 8(a) and (b), we can note that due to the assumption of isotropic relaxations, there is no inversion in the Cd-O distances and, in consequence, the strong change in the EFG components does not take place. In particular, there is not change in sign and orientation of V 33 as in the case of our free relaxation. At the equilibrium position we obtained for V 33 a value of -4.46 × 10 21 V/m 2 (pointing in [001]-direction) and a high η value of 0.91, confirming that an isotropic relaxation is not consistent with the experimental data.

IV. CONCLUSIONS
In this work we have studied through a series of first-principles calculations the problem of a Cd impurity substitutionally-located at the cationic site in rutile TiO 2 . The main result of our work, i.e. that Cd introduces in the host fairly anisotropic relaxations of its nearest oxygen neighbors and that this produces a change of orientation of V 33 from the [001] to the [110] direction when a Ti atom is replaced by a Cd atom in pure TiO 2 , was briefly presented in a recent work with the experimental confirmation of the last prediction 16 .
In this work we have presented details about the electronic structure of the different impurity systems considered. We have considered atomic relaxations and electronic structure self-consistently and have obtained that both aspects of the problem interact with each other. We obtained that atomic relaxations are different for the charged and neutral state of the impurity and that, on the other side, the relaxation process produces a drastic variation in the asymmetry of the charge distribution near the probe-atom for a given charge state of the impurity, which is detected in the strong variation of the EFG tensor. We have shown that the huge difference in the values of the asymmetry parameter η between the charged and neutral state of the impurity arises because of the filling of the impurity level at the Fermi energy. This difference in the η value determines, through comparison with experiment, that Cd is in a charged state when it is introduced as impurity in TiO 2 at room temperature. From these results we have confirmed that the EFG tensor is a very useful magnitude because it is sensible to subtle details of the electronic structure and it can be determined experimentally with high resolution. We have performed a series of checks of the accuracy of the present calculations in order to show that all the predictions of this work are the same if an increment in the basis size, k-mesh or size of the SC are considered or if a different exchange-correlation potential is used. We have also shown that considering relaxations beyond nearest neighbors does not produce any qualitative change in our results. Finally, we have checked that the hypothesis of isotropic relaxations and the use of PCM approximations give results incompatible with experiments and with our calculations. From our results it is clear that the problem of the EFG at Cd impurities in TiO 2 is too complicated to be described (even approximately) by simple models like PCM, antishielding factors and isotropic relaxations. We can conclude that a proper theoretical description of electronic properties of metal impurities in oxide semiconductors should consider self-consistently the chargestate of the impurity and the impurity-induced distortions in the host, specially in the first shell of neighbors of the impurity.