Investigation of pure and Co2+-doped ZnO quantum dot electronic structures using the density functional theory: choosing the right functional

The electronic structures of pure and Co2+-doped ZnO quantum dots (QDs) with sizes up to 300 atoms were investigated using three different density functional theory approximations: local spin density approximation (LSDA), gradient-corrected Perdew–Burke–Ernzerhof (PBE) and the hybrid PBE1 functionals with LANL2DZ pseudo-potential and associated basis set. Qualitative agreement among the three methods is found for the pure ZnO nanostructures, but only the hybrid functional reproduces the correct bandgap energies quantitatively. For Co2+-doped ZnO QDs, both LSDA and PBE incorrectly model interactions between Co2+ d levels and the valence band of ZnO, which will strongly impair predictions of dopant–carrier magnetic exchange interactions based on such computational results. Experimental observations are reproduced well in calculations at the hybrid PBE1 level of theory, making this the method of choice for exploring the magnetism of transition metal ions in ZnO QDs computationally. The qualitative features of the Co2+ 3d levels do not change appreciably with changes in cluster size over the range examined, leading to size-dependent dopant-band edge energy differences. The results presented here thus provide an experimentally calibrated framework for future ab initio descriptions of dopant–carrier and dopant–dopant magnetic exchange interactions in diluted magnetic semiconductors (DMS) nanocrystals.


Introduction
Diluted magnetic semiconductors (DMSs) have been the subject of intense theoretical and experimental investigation recently as potential materials for spintronics applications [1]- [4]. Dopant-carrier exchange interactions in DMSs give rise to 'giant' Zeeman splittings of the semiconductor band structure that may someday introduce new possibilities for information processing via control of carrier spins [5]. Of particular interest to emerging technologies are DMS nanostructures, including DMS quantum dots (DMS QDs) [6]. DMS QDs have been prepared by both molecular beam epitaxy [7] and colloidal synthesis techniques [8]- [11]. Experiments have shown DMS QDs to possess the attractive exchange-induced band structure Zeeman splittings of their bulk counterparts, with quantum confinement contributing additional electronic structure perturbations.
In parallel with experimental advances have been theoretical studies of DMS nanocrystals by ab initio computational methods. Sapra et al [12] investigated quantum confinement effects in Ga 1−x Mn x As nanocrystals using combined tight-binding and density functional theory (DFT) methods, and proposed a size-dependent hole-Mn 2+ binding energy that depends critically on the energy difference between the valence band (VB) edge and the Mn 2+ acceptor level. Huang et al [13] calculated the magnetic properties of Mn-doped Ge, GaAs and ZnSe nanocrystals using real space ab initio pseudo-potentials constructed within the local spindensity approximation (LSDA) and also described Mn-related impurity states becoming much deeper in energy with decreasing nanocrystal size. Using magnetic circular dichroism (MCD) spectroscopy, such size-dependent impurity binding energies were experimentally detected in Co 2+ : ZnSe DMS QDs by Norberg et al [14], who also compared the experimental energies to those calculated by LSDA ab initio methods. This comparison revealed that the calculations reproduced the trends correctly but did not reproduce the energies quantitatively, and therefore did not accurately predict at which sizes the size effects might become important. From a practical standpoint, DFT methods have also contributed in important ways to the ongoing discussion about how dopants are incorporated into colloidal semiconductor nanocrystals during synthesis [15,16].
The accuracy of theoretical descriptions of DMSs depends strongly on the applied formalism and computational strategy. Complications arise from the fact that the delocalized band structure of the semiconductor and the localized 'defect-like' magnetic ion electronic structures should both be described accurately. Historically, bulk semiconductors have often been modeled using a plane wave basis within the local density approximation (LDA) or generalized gradient approximation (GGA) of the DFT. Meanwhile, the energetics and physical properties of molecular transition metal complexes are better described with either meta-GGA [17] or hybrid DFT functionals [17,18], in which a certain amount of Fock exchange is added to part of the local or semilocal density functional exchange energy. Hybrid functionals might be anticipated to perform better than LDA or pure DFT functionals, but the accuracies of various DFT functionals for studies of DMS nanocrystals have not yet been systematically investigated. With few exceptions, past ab initio investigations of DMS nanostructures have lacked experimental calibration, and in most cases large systematic underestimation of bandgap energies has been tolerated. When modeling nanocrystalline DMSs, additional challenges arise from the sensitivity of the results to the structural design of the nanocrystal. In order for ab initio methodologies to contribute significantly to the practical understanding of magnetic exchange interactions and confinement effects in DMS nanostructures, reliable theoretical methodologies are needed that accurately describe both the delocalized and the localized electronic structure features of such DMSs.
In this paper, we examine different density functional theoretical approaches to the description of ZnO QDs with sizes up to 300 atoms. To execute these studies, a pseudohydrogen surface termination methodology [19,20] has been implemented within the Gaussian computational package that allows non-physical surface states to be shifted out of the bandgap energy range. These different DFT-based functionals were then applied to Co 2+ -doped ZnO QDs for comparison with experimental descriptions of this DMS. By introducing only one Co 2+ ion into the ZnO lattice, we have focused on description of the Co 2+ d levels and their interaction with the semiconductor host. The results presented here provide a framework for an experimentally calibrated ab initio description of DMS nanocrystals that can be used to address dopant-carrier and dopant-dopant magnetic exchange interactions in DMS nanocrystals in future investigations.

Methodology
All calculations were carried out using the development version of the Gaussian computational packages [21] with our implementation of pseudo-hydrogen atoms as described below. We compared selected density functionals from three major DFT categories: the LSDA [22]- [25], pure generalized gradient (Perdew-Burke-Ernzerhof (PBE)) [26,27] and its hybrid version (PBE1) in which 25% of Hartree-Fock exchange is included in the exchange functional [26,27]. The Los Alamos double-ζ pseudo-core potential (LANL2DZ) [28]- [31] and associated basis set was used for all the atoms with the Zn (4s, 3d), Co (3s, 3p, 4s, 3d), and O (1s, 2s, 2p) electrons treated explicitly. Although generated from the Hartree-Fock atomic calculations, this pseudopotential was shown to perform well for the DFT-based calculations of the first row transition metal complexes, actinides and lanthanides [32]- [35].
We constructed five different-sized ZnO QDs with the wurtzite structure and nearly spherical shapes as shown in figure 1. The diameters of the QDs range from ∼4 Å (Zn 6 O 6 ) to ∼18 Å (Zn 141 O 141 ). The lattice parameters were obtained from experimental data [36]: a (Å) = 3.249, c (Å) = 5.204 and u = 0.382. Our initial QDs were built to have surface atoms with no more than two dangling bonds. Calculation of a ZnO QD with unsaturated bonds results in unphysical surface states inside or near the bandgap of the semiconductor. In addition, convergence of the self-consistent field (SCF) solution for a bare QD is more difficult to obtain due to excessive flexible electronic degrees of freedom. To solve the problem of surface states, we implemented a capping scheme with pseudo-hydrogens. As proposed by Wang and Li [19] and Chelikowsky et al [20], a surface atom with the formal valence charge m is bound to a hydrogen with nuclear charge q = (8 − m)/4. In the current study, hydrogen atoms with nuclear charges of 0.5 and 1.5 are used to terminate surface O 2− and Zn 2+ ions, respectively. The O-H (1.057 Å) and Zn-H (1.731 Å) bond distances were taken from optimized ZnH 4 and OH 4 tetrahedra. The surface capping results in states from pseudo-hydrogen atoms that are well separated from the bandgap: ∼5 eV below the topmost filled orbital and ∼3-4 eV above the lowest empty one. Although we chose to use the experimental lattice structure in this paper, a full geometry optimization of a pseudo-hydrogen-capped Zn 33 O 33 cluster leads to a change of ZnO bond length by 3% or less, depending on the functional. This indicates that the capping scheme used in these calculations has little effect on the electronic and geometric properties of the ZnO QDs while stabilizing the SCF process.

Pure ZnO QDs
The experimentally measured low-temperature bandgap of bulk ZnO is about 3.4 eV [37,38]. Theoretical studies of bulk ZnO yield a broad range of estimated bandgap energies, however, ranging from less than 1.30 eV by LDA and pure generalized gradient density functionals [39]- [42] to 3.2-3.7 eV by self-interaction-corrected LDA [39,42] and hybrid functionals [40,43]. To determine the accuracy of different density functionals for the QDs of figure 1, the results of our calculations were compared with recent experimental data [44], and extrapolated to the bulk limit. Figure 2 plots calculated bandgap energies (E g ) as a function of the diameter of the ZnO QD. Experimental data for ZnO clusters over the same range of sizes are also included in figure 2 [44]. The ground state energy gap ( E) between the highest occupied oneelectron molecular orbital (HOMO) and the lowest unoccupied one-electron molecular orbital (LUMO) is often used as an approximation of the bandgap of a material. This can be a good approximation in a bulk structure when a single electron promotion does not greatly affect the overall band structure. However, since the bandgap energy is usually measured from UV-vis absorption experiments, and the first excited states in the QDs have excitonic character, the Exp. [44] Exp. [ [44] were collected at room temperature, whereas the calculations relate to low-temperature energies. In ZnO QDs, E g increases by ∼0.1 eV upon cooling from room temperature to cryogenic temperatures [45].
electron-hole interaction can play an important role in the experimental energies of these excited states. This electron-hole interaction is not included in a calculated HOMO-LUMO energy gap.
In figure 2, we compare the commonly used HOMO-LUMO gaps ( E) with the experimental bandgaps. For all methods considered here, the HOMO-LUMO gap ( E) decreases with increasing size of the cluster. This behavior is understood as the quantum confinement effect. Although quantum confinement behavior is observed in the results calculated by all methods, the LSDA and pure PBE methods underestimate the bulk limit bandgap severely compared to experimental values. On the other hand, the hybrid PBE1 method predicts a decreasing E approaching the experimental bandgap of bulk ZnO. The electron-hole interaction plays a more important role in smaller QDs, as anticipated. For small clusters, the HOMO-LUMO gap is not a good approximation because it does not account for the electron-hole interaction. As a result, the deviation between experimental and theoretical results is more pronounced for smaller clusters, but experimental uncertainties in cluster size are also the greatest for these data. As the size of the cluster increases, electron-hole interaction becomes smaller compared to the total energy of the system. The HOMO-LUMO gap ( E) and the experimental value (E g ) therefore converge at large sizes. Figure 3 illustrates the comparison of the density of states (DOS) for a range of ZnO nanoclusters obtained with the three DFT functionals. The DOS is plotted as the number of one-electron orbitals within a bin size of E = 0.2 eV. Vertically through figure 3, the plots PBE and (c) PBE1 with the LANL2DZ pseudo-core potential and associated basis set (energy relative to the VB maximum). Experimental photoemission spectrum (d) for bulk ZnO is adapted from [40,46]. describe QDs ranging from 12 to 282 atoms (not including pseudo-hydrogen atom capping) for a given functional. The VB, spanning from 0 to −5 eV, is mostly made up of O 2p orbitals. The Zn 3d band is deeper in energy, ranging from ∼ −7 eV to -10 eV. The conduction band originates mainly from Zn 4s orbitals with about 10% mixture of O 2p character. As seen from figure 3, discrete energy levels are broadened into continuum bands and the position of the Zn 3d band shifts towards lower energy with increasing QD size. The results for the VB DOS of larger ZnO QDs agree well with photoemission results for bulk ZnO [46]- [50] as illustrated in figure 3(d). When the DOS plots for different density functionals are compared (horizontally through figure 3), the widths and positions of the VB relative to the HOMO are found to be quite similar for the three density functionals used here. All bands are very similar in shape and in characteristic atomic orbital component, regardless of the absolute position. The hybrid PBE1 method affords a larger bandgap and a higher energy 3d band.  The LSDA, PBE and PBE1 methods used in combination with the LANL2DZ pseudopotential and corresponding basis set thus all predict qualitatively similar electronic structures for undoped ZnO QDs. The hybrid method results in a substantially larger (and closer to the experimental) bandgap compared to the LSDA and GGA. Because the bandgap energy of a semiconductor is one of the main factors defining the spectroscopic and other physical properties of DMS nanocrystals based on that host semiconductor, a correct description of this value is important.

Co 2+ -doped ZnO QDs
To investigate the properties of transition metal-doped ZnO clusters, a Co 2+ dopant is introduced by replacing a Zn 2+ ion. The Co 2+ ion sits in the core of the QD at the largest possible distance from the surfaces. The ground state of Co 2+ is high spin with five spin up and two spin down d-electrons (total S = 3/2). In the pseudo-tetrahedral crystal field cation environment of wurtzite ZnO, electrons with the same spin are split by symmetry into two lower energy e (d 2 x−y and d 2 z ) and three higher energy t 2 levels (d x y , d yz and d x z ) [51]- [53], while exchange interactions lead to the splitting between electrons with different spins. Figure 4 shows the DOS for Co 2+ -doped ZnO QDs, with the alpha (spin-up) and beta (spin-down) spins separated. To clarify the location of the Co 2+ d-electrons, the projected DOS 8 for these molecular orbitals are shaded and magnified. The valence and conduction bands of ZnO in the Co 2+ -doped clusters are very similar to those of the pure ZnO clusters for the three functionals studied. The treatment of the Co 2+ 3d-electrons, however, is very different for the pure and hybrid DFT functionals.
Both LSDA and PBE functionals predict that most of the Co 2+ 3d-electrons are located inside the ZnO bandgap, and are localized around the Co 2+ ion. This agrees with earlier theoretical studies of bulk Co 2+ -doped ZnO [52], [54]- [56]. The crystal field splitting between the occupied Co 2+ e and t 2 d-orbitals is about 0.5 eV for spin-up electrons in the LSDA and pure PBE theories. These are not resolved as separate spin-up peaks in figure 4. The non-hybrid DFT functionals have a large (∼2 eV) exchange splitting, leading to well-separated spin-up and spin-down Co 2+ d levels.
In sharp contrast to the results computed by LSDA and GGA functionals, the hybrid method predicts a substantially different picture of the Co 2+ d-electron distribution. The overall energies of the occupied Co 2+ d-orbitals are much lower than those computed with LSDA and PBE. The exchange splitting between spin-up and spin-down Co 2+ d-electron energies is less than 0.5 eV, reducing the energy gap between the spin-up and spin-down peaks. These results agree qualitatively with the recent self-interaction-corrected LDA calculations of Toyoda et al [55] and the LDA+U computations by Hu et al [56]. As seen in figure 4, the hybrid PBE1 yields more covalent delocalization of the five Co 2+ spin-up d-electrons into the VB of the ZnO cluster than found with the other functionals, while the spin-down d-electrons are localized around the Co 2+ ions, with energies just above the ZnO VB edges. Figure 5 illustrates PBE1 spin-up molecular orbitals near the top of the VB (a) and around −4 eV ((b) and (c)), of a doped Zn 83 CoO 84 QD. From these plots, the peak near −4 eV is a bonding molecular orbital involving Co 2+ 3d and O 2p levels, while the peak at the top of the VB is the anti-bonding-type molecular orbital. Interestingly, some of the deep orbitals with high Co 2+ d character are more localized than others (figure 5(b) versus figure 5(c)). Figure 5(d) shows the spin-down molecular orbital at the edge of the VB with much more localized Co 2+ 3d character.
Recently reported photoemission studies on bulk Co 2+ -doped ZnO [49,57,58] describe two key results: (1) pure 3d Co 2+ states are located at the edge of the ZnO VB, well separated from the conduction band; (2) a broad satellite peak is observed at about 4 eV into the VB. While the pure DFT methods shift the occupied Co 2+ 3d levels and put them too close to the conduction band to be consistent with experiment, the hybrid PBE1 correctly places the occupied Co 2+ 3d levels above the ZnO VB edge, cleanly separated from the conduction band. Because the spinup and spin-down levels occur at similar energies for the hybrid calculations, these levels could be observed as a single broad Co 2+ 3d band extending above the VB edge in the photoemission studies. The experimentally observed satellite peak is described as due to highly hybridized Co 2+ 3d and O 2p states. This deeper peak is observed in the calculated results only when using the hybrid PBE1 functional.
Importantly, the calculated results in figure 4 show that while the ZnO band edges shift due to quantum confinement, the Co 2+ d-electron levels remain largely invariant over the same range of quantum confinement, resulting in size-dependent energy differences between Co 2+ and ZnO orbital energies. This description is consistent with experimental results showing no measurable change in the Co 2+ ligand field transition energies between nanocrystalline and bulk Co 2+ -doped ZnO [10,14,59,60], and recent experimental MCD observations and theoretical LSDA results showing pinned Co 2+ d levels but size-dependent band edge energies in Co 2+ -doped ZnSe QDs [14]. An interesting observation can be made regarding the separation between the occupied and unoccupied spin-down levels, an energy that is often related to the tetrahedral crystal field splitting parameter, Dq [52,53]. The energy between the occupied Co 2+ d spin-down e and unoccupied Co 2+ t 2 levels is less than 1 eV in LSDA calculations, and grows to ∼6 eV in the calculations performed with the hybrid functional. The mixing of HF exchange into DFT changes the DFT virtual orbitals from having energies associated with single electron excitations to rather having energies associated with electron affinities [61]. Therefore, the separation between occupied Co 2+ d e and unoccupied Co 2+ t 2 levels should not be interpreted as directly yielding the crystal field splitting parameter when results obtained from a hybrid functional are discussed. In order to interpret the calculated Dq parameter correctly, a different method should be used, for example TDDFT, which takes into account the electron-hole interactions and provides results comparable to the spectral observables. Indeed, the first excited state calculated with linear response TDDFT/PBE1 hybrid functional associated with the Co 3d-3d electronic transition is only ∼0.9 eV above the ground state (to be presented in a forthcoming paper), in reasonable agreement with the experimental 4 A 2 (F) → 4 T 2 (F) transition energy of ∼0.5 eV [59,60]. From the ligand field theory, the energy of this transition is simply 10Dq, meaning these calculated and experimental energies both yield values of Dq below 1.0 eV for Co 2+ -doped ZnO QDs and bulk crystals. A more complete description of optical transition energies in various TM 2+ -doped ZnO QDs calculated using the methods developed here will be the subject of a forthcoming publication.

Conclusion
In summary, the electronic structures of pure and Co 2+ -doped ZnO QDs have been investigated using three different DFT functionals: pure and gradient-corrected (LSDA and PBE), and hybrid PBE1 with 25% Hartree-Fock exchange. To achieve SCF convergence and eliminate surface states in the bandgap, a procedure to create nearly spherical clusters capped with pseudo-hydrogen atoms having modified nuclear charges was implemented within the Gaussian computational package. Our results show that smaller QDs retain most of the electronic structure features of larger QDs and that, except for the bandgap, the computed electronic structure features of pure ZnO clusters are similar for the three density functionals considered. LSDA and PBE computations of the HOMO-LUMO gaps are also similar. However, both energies are substantially lower than experimental bandgap energies at the ZnO bulk limit. In contrast, the hybrid PBE1 method yields bandgap energies that approach the correct asymptotic bandgap energy of ZnO in the bulk limit. While the PBE1 HOMO-LUMO energy difference approaches experimental values for large clusters, differences between the two increase as the cluster size decreases, demonstrating that the electron-hole interaction needs to be taken into account for smaller clusters. Thus, the approximation of the HOMO-LUMO energy difference being equal to the bandgap is generally inaccurate for ZnO nanostructures in the quantum confinement regime.
For Co 2+ -doped ZnO QDs, the qualitative features of the Co 2+ 3d levels do not change appreciably with changes in cluster size over the range examined. Both LSDA and PBE incorrectly model interactions between Co 2+ d levels and the VB of ZnO, which ultimately limits the usefulness of these approaches for predicting or interpreting magnetic exchange interactions in DMS nanostructures. The experimental observations are reproduced well in calculations at the hybrid PBE1 level of theory, making this the method of choice for exploring the magnetism of DMS nanocrystals computationally.