Optical absorption and luminescence energies of F centers in CaO from ab initio embedded cluster calculations

We calculated the optical absorption and luminescence energies of electrons trapped at oxygen vacancies in CaO using a consistent embedded cluster method which accounts for the long-range polarization effects and partial covalence of CaO. Optical absorption and luminescence energies of neutral F center and positively charged F+ center vacancies are calculated by means of time dependent density functional theory using the B3LYP exchange-correlation density functional. Our results demonstrate that using large basis sets to describe a diffuse nature of excited states, and properly accounting for long-range polarization induced by charged and excited defect states, is crucial for accurate predictions of optical excitation and luminescence energies of these defects. © 2006 American Institute of Physics. DOI: 10.1063/1.2337292


I. INTRODUCTION
Alkali earth metal oxides are used as supports for heterogeneous catalysis and photoelectrolysis, to grow high-T c superconductors, as anticorrosion coatings, and as ceramics and materials for microelectronics. [1][2][3][4][5] It is well established that point defects strongly affect the physical and chemical properties of these materials and in some cases are crucial for their applications. Oxygen vacancies are always present in metal oxides as isolated defects and in combination with metal vacancies and impurity ions. Formation of oxygen vacancies in wide gap oxides usually leads to the creation of defect states in the band gap of the material. Often these defects absorb light in the visible range, which may result in the sample coloration. Historically, such anion vacancy related color centers 6 are named F centers from the German word for color ͑farbe͒. By analogy, the notation of F centers is extended also to the vacancies, which do not absorb light in the visible range ͓e.g., MgO ͑Ref. 7͔͒. We note that this notation is common in physical community, whereas solid state chemists prefer a system of Kröger-Vink notations ͑see, for example, Ref. 8͒.
Several charge states of F centers are distinguished in alkali earth oxides, such as MgO and CaO. A vacancy of a neutral oxygen atom O 0 is referred to as an F center, a vacancy of O − ion is called an F + center, and the notation of F 2+ center is reserved for a vacancy of an O 2− ion. F centers in cubic ionic materials are characterized by having a maximum of electron density of electron͑s͒ localized inside a vacancy in the defect ground state. 9,10 However, in a range of other oxides electrons tend to rather localize on cations surrounding the oxygen vacancy ͑see, for example, review in Ref. 11͒.
In the case of CaO, the relative concentrations of F and F + centers generally depend on the method used in the defect production. For example, primarily F + centers are generated in electron-and neutron-irradiated samples, while predominantly neutral F centers are formed in thermochemically reduced CaO crystals.
F + centers are paramagnetic and have been studied using magnetic resonance techniques. 12 Optical absorption and emission spectra of both F and F + centers in CaO have also been extensively studied experimentally. [13][14][15][16][17][18][19] The optical absorption energies observed for F center ͑3.1 eV͒ and F + center ͑3.7 eV͒ correspond to 1 A 1g → 1 T 1u and 2 A 1g → 2 T 1u transitions, respectively. Optical excitation of these centers induces luminescence at 2.1 and 3.3 eV, respectively, for F and F + centers. The emission spectrum for F center can be obtained only by providing a sufficiently large excitation density at 3.1 eV. 14,16 It has been concluded that the 2.1 eV emission band is associated not with the allowed 1 T 1u → 1 A 1g transition, which has not been observed, 16 but with the spin-forbidden 3 T 1u → 1 A 1g transition. Such interpretation is consistent with a rather large Stokes shift ͑ϳ1 eV͒ observed for F center. In the case of F + center, the emission band at 3.3 eV has been associated with the spin allowed 2 T 1u → 2 A 1g transition, which exhibits a smaller Stokes shift of 0.4 eV.
While the spectroscopic characteristics of F centers in CaO are well established experimentally, there have been only a few theoretical studies of their electronic structure and properties. 20 cluster models. In this work we study the optical absorption and photoluminescence of F and F + centers in the bulk of CaO using a much more refined embedded cluster approach and density functional theory ͑DFT͒. The main challenges of this work concern building a consistent embedded cluster scheme for CaO, calculating optical excitation energies into delocalized excited defect states, and calculating luminescence energies. We demonstrate that using large basis sets to describe a diffuse nature of excited states, and properly accounting for long-range polarization induced by charged and excited defect states, is crucial for accurate predictions of optical excitation and luminescence energies of these defects. This paper is organized as follows. The embedded cluster scheme and the details of its parametrization are described in Sec. II. The electronic structure of the F centers and the calculations of the optical absorption and luminescence energies are presented in Sec. III. The conclusions are given in Sec. IV.

II. MATERIAL MODEL AND COMPUTATIONAL DETAILS
It has been demonstrated that the electronic structure of ionic crystals such as MgO can be studied successfully using cluster models where part of the crystal is modeled by a small number of atoms treated quantum mechanically while the rest of the solid is included in a more approximate way. The most common approaches embed the quantummechanical ͑QM͒ cluster into a large set of point charges ͑PCs͒ to account for the long-range electrostatic interactions. 23 An interface of total ion potentials ͑TIPs͒ between the quantum-mechanical cluster and the set of point charges is included to avoid artificial polarization of the electron density of the cluster towards the PCs and to account for the short-range interactions of the atoms of the cluster with the surroundings. 24,25 This approach does not include the long-range polarization of the lattice in response to the perturbation induced by a defect. This effect is particularly important for modeling charged defects and excited defect states. In the present work these effects are taken into account by embedding QM cluster into a polarizable lattice represented by a shell model ͑Refs. 26 and 27, and references therein͒ and treating the cluster electronic structure and the lattice polarization in a self-consistent way. [28][29][30] Another issue concerns choosing an accurate enough electronic structure method for achieving the best accuracy in predicting defect spectroscopic properties. In the recent years, a series of studies has been published aimed at describing the optical spectra of various oxygen vacancies in MgO bulk and surface F centers using configuration interaction ͑CI͒ based methods. [31][32][33][34] These calculations used an embedded cluster model that accounts for the Madelung potential and ion-size effects. They demonstrated that to accurately predict the optical spectra of the F, F + , and M centers in MgO one needs to employ extensive basis sets to account for electronic correlation effects through large CI expansions with a concomitant huge computational cost. On the other hand, extensive calculations of the bulk and surface defects in MgO, 34 SiO 2 , 35 and ␣-Al 2 O 3 ͑Ref. 36͒ have demonstrated that one can obtain quite accurate optical transition energies of these defects using DFT with so called hybrid functionals 37,38 and time dependent density functional theory 39,40 ͑TD-DFT͒ for calculating the excitation energies. In this work we extend this approach to studying defect properties in CaO.

A. Embedded cluster approach
Within the embedded cluster model used in this study, the bulk CaO crystal is represented atomistically using a cubic 24ϫ 24ϫ 24 ion cluster, referred to as a nanocluster. The nanocluster is divided into two regions: a spherical region I at the center of the nanocluster and the remaining partregion II. Ions in region II are represented using a classical rigid ion model, they are fixed in the ideal lattice positions and serve to reproduce correctly the potential energy surface and the electrostatic potential at the boundary of regions I and II. Region I has the radius of just over 11 Å and includes 437 ions. The innermost part of region I, containing up to several tens of ions, is considered quantum mechanically ͑QM cluster͒. The rest of region I includes up to several hundred ions represented using the shell model. 26 This simple mechanical approach allows one to simulate the dipole polarizability of the ions. The QM cluster is linked to the rest of the solid through an interface region. The interface ions have a dual nature: their interaction with the ions of the QM cluster is considered quantum mechanically, while their interaction with the ions of the shell model region is represented using classical interatomic potentials.
The total energy of the system is calculated as a sum of several contributions. 41,42 ͑i͒ The total energy of the QM cluster and the interface region in the electrostatic potential of all classical ions, is calculated using an ab initio method such as Hartree-Fock ͑HF͒ and DFT. ͑ii͒ The interaction energy of the classical atoms with each other, including contributions due to electrostatic and short-range interactions and ionic polarization. ͑iii͒ The short-range part of the interaction between the classical atoms and those of the QM cluster and of the interface.
The short-range interactions are represented using pairwise or many-body potentials, which depend only on the atomic coordinates. In this work we use the Buckinghamand Morse-type potentials: The parameters of interatomic potentials ͓A, , and C in Eq. ͑1͒ and A, r 0 , and ␣ in Eq. ͑2͔͒ are fitted ͑vide infra͒ to mimic accurately the non-Coulomb terms of the total energy arising from the Pauli repulsion and van der Waals attraction. The total energy of the system is minimized simultaneously with respect to positions of all species in region I: classical cores and shells and nuclei of quantum-mechanical ions. This scheme allows us to account self-consistently for both the defect-induced lattice polarization and the effect of the lattice polarization on the electronic structure of the defect itself. This is different from the usual embedding scheme based on the use of point charges 43,44 and total ion potentials [45][46][47] where the effect of lattice response is not included.
The scheme employed in this work is implemented in a computer code GUESS, 41,48 which is interfaced with a package for quantum-mechanical calculations used to compute the contributions to the total energy and forces. This particular realization of the embedded cluster approach has been previously applied to studying defects in NaCl, KBr, MgO, Mg 2 SiO 4 , and 12CaO · 7Al 2 O 3 . [49][50][51][52][53] In the following we give the details of the electronic structure calculations of the embedded QM cluster ͑Sec. II B͒ and outline the procedure used to provide a consistent description of the classical and quantum-mechanical parts of the system ͑Sec. II C͒.

B. Computational details: QM cluster
Most of the results are obtained for a stoichiometric but not symmetrical QM cluster Ca 14 O 14 ͑see Fig. 1͒, which was used for the reasons described below in Sec. II C. A comparison to the results obtained using a nonstoichiometric but symmetrical cluster Ca 14 O 13 is also made. Hybrid density functional, which combines Becke's three parameter exchange functional 37 with a correlation functional by Lee, Yang, and Parr 38 ͑B3LYP͒, has been used for the electronic structure calculations. The QM contributions to the total energy and forces are calculated using the GAUSSIAN98 code, 54 which is interfaced with the GUESS code. The calculations of optical excitation energies were performed using the TD-DFT method as implemented in the GAUSSIAN98 package. Several Gaussian-type orbital ͑GTO͒ basis sets were used for the ions in the QM cluster. These were obtained as extensions of the standard 6-31G basis set as described in detail in Sec. II D. The interface Ca ions were represented using the large core effective core pseudopotential ͑ECP͒ LANL1 ͑Ref. 55͒ so as to reproduce the ion-size effect. These will be called total ion potentials ͑TIPs͒ below. In some calculations we also used the small core ECP LANL2 ͑Ref. 56͒ on Ca ions.

C. Implementation of consistent embedded cluster scheme
The simple scheme outlined in Sec. II A provides an adequate embedding potential and good description of the lattice for strongly ionic materials, such as NaCl and MgO. However, CaO is more covalent than MgO and whether one can use the same approach here is not clear a priori. To achieve consistency between the quantum and classical parts of the system and to account for larger degree of covalency in CaO as compared to MgO we used the following approach.
Firstly, reference calculations were carried out using the periodic CRYSTAL code, 57 implementing the same B3LYP functional and Gaussian-type atomic basis as used in our embedded cluster calculations. These calculations are required to obtain the lattice constant and the crystal band structure. Different basis sets reported by Habas et al. 58 ͑see Table I͒ have been used with up to 8s, 13sp, and 3d basis functions for Ca and up to 8s, 6p, and 1d functions for O atoms, respectively. The CaO lattice constant has been optimized using several basis sets ͑see Table I͒. We found that to achieve a good agreement with the experimental value ͓4.81 Å ͑Ref. 59͒ at room temperature͔, large all electron basis sets are needed with the best result being 4.818 Å. This value has been taken as a reference in further calculations and used to generate the nanocluster representing CaO bulk.
Secondly, we fitted the parameters of Buckingham interatomic potentials describing the interactions between classical CaO atoms to reproduce the lattice constant, obtained from periodic calculations, and other crystal properties. The electrostatic interaction of classical ions in ionic compounds is usually parametrized using formal ionic charges, e.g., ±2.0͉e͉ for alkali earth oxides. However, in spite of the largely ionic character of CaO, 60 it is expected that a covalent contribution to the bonding can be accounted for more accurately if the values of the charges are scaled from the formal ones in such a way that they are approximately equal inside and outside the QM cluster. We note that the use of scaled charges requires that a stoichiometric QM cluster is used in the calculations to ensure that a correct number of electrons is attributed to it.
A convenient procedure to derive scaled ionic charges is to use the results of a population analysis of a charge density obtained in quantum-mechanical calculations. For this pur-   61 ͑NPA͒ were averaged over anions and cations of the QM cluster, respectively, and then assigned to classical ions to calculate the electrostatic potential in the whole system. The procedure has been repeated several times until the value of the ionic charge has converged to ±1.74͉e͉, which was used in all further calculations unless specified otherwise. Note that this charge is smaller than ±1.85͉e͉ obtained for MgO in a similar manner.
Then the parameters of the short-range interactions between classical ions have been fitted to reproduce several properties of bulk CaO ͑see Tables II and III͒ including the bulk modulus ͑B͒, lattice constant ͑a͒, dielectric ͑ 0 and ϱ ͒, and elastic constants ͑C 11 , C 12 , and C 44 ͒. The target value for the lattice constant has been set to that obtained in the periodic QM calculations to minimize the mismatch between the classical and QM parts in the embedding scheme. For simplicity, only O ions were made polarizable, while the Ca ions were considered using the rigid ion model. The short-range part of the interatomic interaction was defined for O-O and Ca-O pairs and represented using Buckingham-type potentials. The van der Waals attraction for Ca-O interaction has been set to zero in analogy to other parametrization schemes, which is consistent with the rigid ion model used for Ca. As an initial guess we used the parameters of the interatomic potentials obtained earlier by Sangster and Stoneham. 62 The generalized utility lattice program 63 ͑GULP͒ was used for the fitting.
Finally, we introduced an additional corrective potential at the interface of the QM cluster. This is necessary to eliminate a small lattice distortion, which appears at the interface even though the interatomic distances in the QM and classical parts are the same. A Morse-type potential is introduced between the nearest neighbor interface Ca ions and O ions of a QM cluster, i.e., within the cutoff radius of 3.0 Å. We found that parameters A = 4.00 eV, ␣ = 1.00 Å −1 , and r 0 = 2.450 Å provide Ca-O distances in good agreement with those optimized using the CRYSTAL code ͑2.4091 Å͒ and minimize the local perturbations.

D. Basis sets used in the quantum-mechanical calculations
Our preliminary calculations have shown that, although the present cluster model predicted lattice parameter is close to that obtained from periodic calculations using the same DFT method and GTO basis sets of similar quality, the optical absorption ͑OA͒ energies of the F centers calculated with the cluster model exhibit a very strong dependence on the basis set. Therefore, we investigated the convergence of the OA energies on the size of the basis set as described below. The standard 6-31G basis set on all atoms in the QM cluster is considered as a reference ͑b1͒. Notice that the b1 set does not include basis functions centered at the vacancy. Basis set b2 adds an sp-type function at the vacancy site with the exponent factor of 0.27 bohr −2 and a single d function with an exponent factor of ␣ 0 ͑d͒ = 0.80 bohr −2 . The addition of extra N d functions at the vacancy site is indicated as b3.N; values of N from 1 to 5 have been used. The exponent factors of these d functions are calculated as ␣ 0 ͑d͒ ϫ ͑0.6͒ N , where 0.6 is a convenient scaling factor. For instance, the basis b3.2 contains the original b2 d function and two additional d functions with exponents ␣ 0 ͑d͒ ϫ 0.6 and ␣ 0 ͑d͒ ϫ ͑0.6͒ 2 . Basis set b4 is generated by extending b3.1 set so as to include polarization functions on the six nearest neighbors of the vacancy site, i.e., use 6-31G * basis for these atoms.
Two extensions of b4 have been considered. Set b5 is generated by using a triple zeta basis ͑6-311G͒ for the six cations next to the vacancy site. Alternatively, b6 defines the basis obtained by extension of b4 so as to use 6-31G * basis set for 12 anions next to the vacancy. Finally, extending set b5 by using a triple zeta plus polarization ͑6-311G * ͒ basis on the six cations next to the vacancy generates set b7. In order to verify that the results are converged with respect to the basis set, an even more extended basis has been considered for these cations, which leaves the three outermost primitive p functions of the b7 basis uncontracted; this is denoted as b8.

E. Tests of the model
To check the applicability of the model described above for studying the ground and excited states of oxygen vacancies in CaO we carried out a series of tests. First, the total  Fig. 1. In this case some distortion is introduced by the asymmetry of the QM cluster but the average Ca-O distances obtained for the asymmetric Ca 14 O 14 cluster are almost the same as those obtained with the symmetric cluster. We note that the parameters of the corrective Morse potential were fitted for 6-31G basis set only. When 6-311G * basis on the six cations next to the vacancy is used instead, the Ca-O distance at the interface becomes smaller than that obtained in a periodic model. At the same time, increasing the basis set ͑see Sec. II D͒ results in a better description of the band gap as represented by the energy of the lowest optical excitation calculated using TD-DFT, which becomes 6.95 eV, i.e., surprisingly close to the experimental value of 7.10 eV. 64 Thus the embedding model adequately describes the electronic structure of CaO. Excited states of oxygen vacancies can be delocalized. Therefore extensive tests have also been carried out to establish a cluster size and basis set necessary for describing these states as well as optical excitation and luminescence energies. The excitation energies and transition matrix elements calculated for the F center using TD-DFT in two QM clusters using different basis sets are summarized in Table IV ͑qualitatively very similar results were obtained for the F + center͒. In these calculations the ions in both QM cluster and shell model regions were fixed at their ideal lattice sites. We neglected the polarizability of O ions and used formal charges of ±2.0͉e͉ for the ions of the classical environment. This allows us to separate the cluster size and basis set effects from other effects, e.g., geometry relaxation. Two QM clusters were considered: the smallest possible cluster was obtained by removing the central oxygen atom in Ca 6 O QM cluster ͑Ca 6 V, where V defines the vacancy͒ and the Ca 14 O 12 V cluster was obtained by removing the central oxygen from the Ca 14 O 13 cluster.
The absorption energies show strong dependence on the basis set ͑see Table IV͒. For example, as the basis set is increased from b1 to b2, i.e., the O basis set is positioned at the vacancy site, the F center absorption energy decreases by ϳ0.2 eV, which is in line with previous studies for MgO. [31][32][33] The adsorption energy decreases further when diffuse basis functions are included. A considerable change in the OA energies is observed when the QM cluster is increased so as to include explicitly 14 Ca and 12 O atoms, surrounding the F center. The inclusion of polarization functions on Ca ions ͑b4͒ results in a further drop of the absorption energy to ϳ4.5 eV. However, when the number of sp functions is increased ͑b7͒ computed excitation energy becomes 3.67 eV. Further extension of the basis set to b8 does not introduce significant changes. Therefore, we suggest that basis set b7 provides the satisfactory variational freedom for calculating excitation energies.
In order to understand the effect of the basis set on the nature of the lowest excited state, we analyzed the expansion coefficients of the highest occupied molecular orbital ͑HOMO͒ and lowest unoccupied molecular orbital ͑LUMO͒ one-electron states over the basis functions for several basis sets including b3.1, b4, b5, b6, b7, and b8. In all cases the HOMO of the F center is formed mainly by s functions centered at the vacancy site with small contributions due to s and p functions of the six nearest Ca ions. In the case of b7 basis set, a small contribution due to d functions centered on these ions also appears. The LUMO state is dominated by p functions centered on the vacancy site and, similar to HOMO, contains small contributions due to s and p functions of the six nearest Ca neighbors for all basis sets except b7 where the d functions on the Ca near the vacancy have an important contribution. These contributions are clearly seen in the LUMO plot in Fig. 2.
These results suggest that a diffuse basis set centered at the vacancy site can be sufficient for describing accurately the ground state of F centers in CaO. However, to reproduce the excited states requires a much more extended basis set not only at the vacancy site but also on the neighboring atoms. These findings are in agreement with the previous discussion on the computation of the CaO band gap and with other studies of ͑CaO͒ n gas phase clusters. 65

A. Electronic structure and optical absorption of F centers
The results reported in this section were obtained using stoichiometric, but nonsymmetric, Ca 14 O 14 QM cluster and basis set b7. To find the equilibrium structure of the F and F + centers their total energies have been minimized with respect to coordinates of all atoms in region I, as described in Sec. II. The calculated defect-induced displacements of nearby atoms are summarized in Table V. In the case of the F center, QM cluster is used instead, the lattice relaxation is symmetric and the displacements of the Ca ions from the vacancy site are ϳ1% for the F center and ϳ5% for the F + center. The displacements of the ions located more than 4 Å away from the vacancy are almost negligible, in particular, for the F center with displacements of less than 0.001 Å. For the charged F + center the relaxation of those distant ions is small but still significant ͑ϳ0.010 Å͒. We would like to mention that for the F center, the atomic displacements computed here are in a good agreement with those obtained in previous first principles calculations based on a periodic approach. 66 The electronic structure of the fully relaxed F and F + centers is shown in Fig. 2. The occupied 1 A 1g state of the F center is at 2.40 eV above the top of the valence band ͑VB͒, while unoccupied 1 T 1u states are in the range of 0.93-1.03 eV below the conduction band ͑CB͒. Notice that, because of the use of an asymmetric cluster model, the threefold degeneracy of the t 1u level in the ground state is slightly broken with a splitting between the resulting levels of 0.03 and 0.10 eV for F + and F centers, respectively. A qualitatively similar electronic structure was found for the F + center with the occupied 2 A 1g state being 1.92 eV above the VB and unoccupied 2 T 1u states at ϳ1.15 eV below the CB.
The excitation energies calculated using the TD-DFT for the most intense transitions of these defects are summarized in Table VI and also indicated in Fig. 2. The results clearly demonstrate that a good agreement with the experimental data can be achieved when the defect-induced lattice relaxation is taken into account and a sufficiently large basis set is used to accommodate the delocalized excited states. In the case of the F center, the excitation energy of the 1 A 1g → 1 T 1u transition ͑3.52 eV͒ is within 0.4 eV or 12% of the experimental value. Better accuracy has been achieved for the excitation energies for F center in MgO bulk using also TD-DFT and relatively small basis sets. 34 However, very large basis sets are required when accuracy of the same order is sought by making use of sophisticated explicitly correlated wave function calculations. 33 In this respect, it is worth pointing out that basis set requirements for the TD-DFT accurate prediction of the optical excitation of the F center in CaO are more stringent than for MgO. Nevertheless, even in FIG. 2. Schematic representation of the location of the energy levels and absorption transitions ͑in eV͒ associated with the F and F + centers ͑left panel͒ and contour plot of the molecular orbitals of the 1 A 1g and 1 T 1u symmetries located around the oxygen vacancy ͑right panels͒. Blue and red spheres stand for Ca and O atoms in the QM region, respectively; white spheres represent part of the interface Ca atoms. The relative positions of the different energy levels are given with respect to the levels representing the valence band ͑VB͒ and conduction band ͑CB͒ of the perfect crystal. TABLE V. Amount of relaxation with respect to the bulk structure ͑⌬r and %͒ of the first six Ca 2+ neighbors of the F / F + center represented by an embedded Ca 14 O 13 V QM cluster ͑Fig. 1͒. The meaning of A, B, and C labels is explained in Fig. 1. The last row includes the results for the symmetric Ca 14 O 12 V cluster. ⌬r is defined as d͑Ca-V͒ def − d͑Ca-O͒ reg , where d͑Ca-V͒ def and d͑Ca-O͒ reg are the Ca-V and Ca-O distances for the defective and regular structures. The percentage of relaxation is calculated with respect to d͑Ca-O͒ reg : % =͑⌬r / d͑Ca-O͒ reg ͒ ϫ 100.

Model
F center F + center Ca 14  F center Expt.: 3.3 eV this case there is a noticeable difference with respect to experiment which could be reduced by using a more accurate treatment of electronic correlation effects in the geometry optimization within the present embedding scheme.
In the case of the F + center, the excitation energy of the 2 A 1g → 2 T 1u transition ͑3.85 eV͒ is within 4% of the experimental value. We note that the difference between the excitation energies calculated for the unrelaxed and relaxed F + center is 0.46 eV, which emphasizes the importance of the lattice relaxation effect on the electronic structure of this charged defect.
Our results suggest that both the TD-DFT calculations and these based on wave function methods require sufficient variational freedom to describe excited states for this kind of defects. Indeed, present calculations for the embedded Ca 14 O 12 V QM cluster based on a complete active space selfconsistent field ͑CASSCF͒ wave function and further including electronic correlation by perturbation theory up to second order 67,68 ͑CASPT2͒ show a significant reduction in the excitation energies when a large basis set for the six Ca atoms nearest the oxygen vacancy is included. For the unrelaxed system, the CASPT2 values obtained using basis set b7 are 4.18 and 4.84 eV for the F and F + centers, respectively, which are ϳ0.4-0.5 eV larger than the corresponding values computed with TD-DFT B3LYP. When the geometry relaxation around the vacancy, as predicted by the shell model, is included, the CASPT2 values decrease by ϳ0.1 eV for the F center and 0.4 eV for the F + center. Therefore, the structural relaxation effect is almost the same for TD-DFT and CASPT2 calculations. The fact that CASPT2 overestimates the excitation energies is due to a limited active space used in these calculations.

B. Luminescence energies of F centers
Accurate calculations of luminescence energies require the knowledge of the potential energy surface of the corresponding excited state. This is relatively straightforward to achieve in the case of the F center. First, the energy of the lowest 3 T 1u excited state, which is responsible for the luminescence, was minimized with respect to the coordinates of all centers in region I using the unrestricted B3LYP density functional. Then the luminescence energy of the 3 T 1u → 1 A 1g transition was calculated using TD-DFT as the 1 A 1g → 3 T 1u excitation energy for the geometry corresponding to the 3 T 1u energy minimum. The obtained luminescence energy of 1.81 eV ͑see Table VI͒ is 0.3 eV lower than the experimental value. The luminescence energy can also be estimated from the total vertical energy difference between the singlet and triplet states at the triplet minimum energy geometry ͑⌬SCF method͒. The result thus obtained is 1.70 eV which as expected is close to the TD-DFT value for the singlet-triplet excitation but somewhat smaller due to the inclusion of the density relaxation in the triplet state whereas this is not the case in the TD-DFT calculations. Notice also that TD-DFT introduces some configuration mixing, which is absent in the ⌬SCF scheme. It is worth mentioning that in the optimized geometry of the excited state, the displacement of the axial Ca is different from the displacement of the equatorial Ca because of the Jahn-Teller effect produced by the partial occupancy of the t 1u ͑p-like͒ orbitals. Notice that strictly speaking the Jahn-Teller effect in the excited state reduces the symmetry of the system and therefore the notation 3 T 1u should not be used. Nevertheless, we find it useful to continue using this term and distinguish the x, y, and z components as illustrated in Figs. 2 and 3.
In the case of the F + center the situation is more complex because the two states involved in the luminescence transition are of the same spin symmetry ͑ 2 T 1u and 2 A 1g ͒. The GAUSSIAN98 code used in these calculations does not provide for calculating forces in excited states using TD-DFT. Therefore, to obtain the luminescence energy, we constructed an approximate potential energy surface of the 2 T 1u state using the following approach. First, we defined a set of local geometries and total energies of the F + center by displacing six Ca ions near the vacancy along the crystal axes equivalent to FIG. 3. Contour plots corresponding to the potential surfaces involved in the emission process of an F + center in CaO. These involve the 2 A 1g ͑a͒, degenerate 2 T 1u ͑x͒ and 2 T 1u ͑x͒ ͑b͒, and nondegenerate 2 T 1u ͑y͒ ͑c͒ states arising from the symmetry breaking originated from the Jahn-Teller distortion ͑see text in Sec. III B͒. ␣ and ␤ correspond to scaling factors which modify proportionally the distance between the six nearest cations and the defective center, taking as reference the geometry relaxed for the 2 A 1g state. ͑d͒ Potential curves for the ground and excited states fixing ␣ to 1.00.

͑4͒
Here ␣ and ␤ are parameters used to scale the distances between the vacancy site and its Ca neighbors in the regular structure ͑d A and d B,C ͒ and d A Ј and d B,C Ј are the corresponding distances in the generated structures. The need for two parameters arises from the nature of the excited state, which shows a Jahn-Teller distortion as commented above for the excited state of the F center. In this way we constructed a grid of 80 points defined by ␣ and ␤ ranging from 0.96 to 1.12. The rest of the atoms present in the system have been frozen at the equilibrium geometry of the 2 A 1g ground state. The potential energy surface of the three electronic states derived from the original 2 T 1u one for the ideal geometry was constructed by adding the calculated excitation energies of 2 A 1g → 2 T 1u transitions for each point on the grid to the ground state energies. The calculated potential energy surfaces for the 2 A 1g electronic ground state, the degenerate 2 T 1u ͑x͒ and 2 T 1u ͑z͒, and nondegenerate 2 T 1u ͑y͒ states are shown in Figs. 3͑a͒-3͑c͒, respectively. The lowest excited state energy minimum was found for the state 2 T 1u ͑y͒ at ␣ = 1.000 and ␤ = 1.024. The corresponding 2 T 1u → 2 A 1g emission energy of 3.79 eV is equal to the excitation energy at the grid point corresponding to the minimum on the excited state potential energy surface. Thus a calculated Stokes shift of 0.1 eV is noticeably lower than the experimental one ͑0.4 eV͒. This discrepancy could be caused by an approximate account for the lattice relaxation in the F + excited state, which does not include the full lattice response to the delocalized excited state.
To roughly estimate a possible error made by this approximation one can consider an extreme case that the excited electron is fully delocalized and does not screen the vacancy. Then the lattice relaxation in the excited state of the F + center could be approximated by that corresponding to the F ++ center. The luminescence energy can be calculated as a 2 A 1g → 2 T 1u vertical optical transition in the F + center at the relaxed F 2+ center geometry. The luminescence energy ͑3.69 eV͒ and the Stokes shift ͑0.16 eV͒ calculated using this approach agree marginally better with the experiment than those described above. However, the origin of such small value for the calculated Stokes shift remains unclear although it points towards an insufficient relaxation of the excited state. Using the same method one can also estimate the effect of weak vacancy screening in the excited state of the F center. In this case one should use the equilibrium geometry found for the ground state of the F + center as the equilibrium geometry for the excited state of the F center. Then the luminescence energy can be calculated as the excitation energy of the 1 A 1g → 3 T 1u transition for this geometry. The obtained value ͑2.09 eV͒ is almost identical to the experimental one ͑2.10 eV͒. When comparing the final geometry of the triplet state of the F center and of the doublet state of the F + center one finds that the geometry relaxation of the six nearest Ca to the vacancy differs by 3%-8% being, of course, larger for the charged defect.

IV. CONCLUSIONS
The optical absorption and luminescence energies of the F and F + centers in the bulk of CaO have been studied using the embedded cluster model. The electronic transitions corresponding to both absorption and emission have been computed by TD-DFT method employing the hybrid B3LYP density functional and compared with the CASPT2 explicitly correlated method. The results suggest that Ca 14 O 13 is the minimum quantum-mechanical cluster required for describing the properties of bulk CaO and predicting reasonably accurate excitation energies for both F and F + centers. Using extensive basis sets to represent the first nearest neighbors to the O vacancies is crucial for obtaining accurate enough results for the excitation of these point defects in CaO. The requirement for very large basis sets is in agreement with previous results on the same centers of MgO carried out using the explicitly correlated CASPT2 wave function method. 33 However, for CaO the requirement is still more stringent because of the more delocalized character of the excited state.
Calculated absorption energies for both F and F + centers are in good agreement with experimental data, especially after taking into account the geometry relaxation around the point defect. This effect is more important for the F + center, as expected from the charged nature of this defect. The calculated emission energies are also in reasonable agreement with the experiment. However, for the F + center the discrepancy with respect to the experimental value is larger due to the approximations made in determining the geometrical relaxation for the excited state involved in the emission process.
To summarize, embedded cluster model calculations that take into account the long-range polarization effects due to the creation of the point defect coupled with hybrid TD-DFT methods are able to reproduce the main spectroscopic features associated with the presence of trapped electrons at oxygen vacancy sites. Quantitative agreement with experiment requires employing very large basis sets on the cations surrounding the vacancy and accounting for the geometry relaxation effects around the point defect.

ACKNOWLEDGMENTS
One of the authors ͑J.C.͒ thanks the Spanish Ministerio de Educación y Ciencia for a predoctoral grant. Financial support has been provided by the Spanish Ministry of Education and Science ͑Project Nos. CTQ2005-08459-CO2-01 and UNBA05-33-001͒ and, in part, by Generalitat de Catalunya ͓Project Nos. 2005SGR-00697 and 2005 PEIR 0051/69 and Distinció per a la Promoció de la Recerca Universitària de la Generalitat de Catalunya granted to another author ͑F.I.͔͒. Another author ͑P.V.S.͒ is supported by Grant in Aid for Creative Scientific Research ͑Grant No. 16GS0205͒ from the Japanese Ministry of Education, Culture, Sports, Science, and Technology. Part of computer time was provided by the Centre de Supercomputació de Catalunya, CESCA, and Barcelona Supercomputing Center, BSC, through generous grants from Universitat de Barcelona, Fundació Catalana per a la Recerca and BSC.