Theory of the sp–d coupling of transition metal impurities with free carriers in ZnO

The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s,p{-}d$$\end{document}s,p-d exchange coupling between the spins of band carriers and of transition metal (TM) dopants ranging from Ti to Cu in ZnO is studied within the density functional theory. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+U$$\end{document}+U corrections are included to reproduce the experimental ZnO band gap and the dopant levels. The p–d coupling reveals unexpectedly complex features. In particular, (i) the p–d coupling constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_0\beta$$\end{document}N0β vary about 10 times when going from V to Ni, (ii) not only the value but also the sign of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_0\beta$$\end{document}N0β depends on the charge state of the dopant, (iii) the p–d coupling with the heavy holes and the light holes is not the same; in the case of Fe, Co and Ni, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_0\beta$$\end{document}N0βs for the two subbands can differ twice, and for Cu the opposite sign of the coupling is found for light and heavy holes. The main features of the p–d coupling are determined by the p–d hybridization between the d(TM) and p(O) orbitals. In contrast, the s–d coupling constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_0\alpha$$\end{document}N0α is almost the same for all TM ions, and does not depend on the charge state of the dopant. The TM-induced spin polarization of the p(O) orbitals contributes to the s–d coupling, enhancing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_0\alpha$$\end{document}N0α.

www.nature.com/scientificreports/ on its charge state can be reflected not only in the magnitude, but also in the sign of N 0 β . Moreover, N 0 β can be different for the light and heavy hole subbands depending on the detailed electronic structure of a TM ion.

TM impurity levels in ZnO.
Several magnetic properties of a TM dopant are determined by its energy levels relative to the VBM and the conduction band minimum (CBM). An exemplary band structure of TMdoped ZnO is discussed in Supplementary Information, see Figs. S2 and S3, while the relevant results, necessary to understand the mechanism of the s, p-d coupling, are given in Fig. 1. We first recall that the d(TM) shell of a substitutional TM ion in a zinc blende crystal is split into a e 2 doublet and a t 2 triplet higher in energy. Both states are spin split by the exchange coupling, stronger than the crystal field, and all TM ions in ZnO are in the high spin state, see Fig. 1e. Next, in ZnO the triplets are further split by the uniaxial wurtzite crystal field into singlets and doublets. This holds also for the p(O)-derived VBM, which is split by 64 meV into a light hole singlet and a heavy hole doublet, denoted by A 1 and E 2 in the following. When the U(TM) corrections are employed, the splitting of a t 2 state depends on its occupation. The U-induced contribution V U to the Kohn-Sham potential is 16 where φ m are the localized d orbitals occupied by σ m electrons, and ψ σ kν are the Kohn-Sham states for the wave vector k, band ν , and spin σ . The V U potential only acts on the contribution of the mth d(TM) orbital to the given (ν, k, σ ) state. As a result of the U(TM) correction, which favours fully occupied or completely empty orbitals over the partially occupied ones 16 , the splitting of the fully occupied t 2 is small, about 0.1 eV, which is close to the crystal field splitting of the VBM. Otherwise, the splitting is more pronounced, and can exceed the crystal field splitting. The calculated energy levels of neutral TM 2+ ions are shown in Fig. 1a. Considering the series Fe-Cu we see that t 2↑ is very close to the VBM. Interestingly, this trend reflects the d-shell energies of isolated TM atoms, see Supplementary Information, Fig. S1. Figure 1c,d show the charge state dependence of the energies of the Ni levels as an example. With the increasing occupation of the d shell the levels shift to higher energies, which is caused by the increasing intrashell �ε c = ε c↓ − ε c↑ , and the VBM, �ε vγ = ε vγ ↓ − ε vγ ↑ . Here, γ denotes the A 1 or the E 2 partner of the VBM. The exchange constants, N 0 α for the s-d coupling and N 0 β γ for the p-d coupling, are obtained directly from those splittings for supercells containing a single TM ion 12 . Since the splittings of the light hole A 1 and the heavy hole E 2 bands are different, we have where �S� = (N ↑ − N ↓ )/2 , the total spin of the supercell, is the difference in the number of spin-up and spindown electrons, N 0 is the density of the cation sites in ZnO, and x is the composition of Zn 1−x TM x O ( x = 0.028 for our supercells). Our definitions imply that both N 0 α and N 0 β are positive for the ferromagnetic (FM) and negative for the antiferromagnetic (AFM) coupling of conduction or valence electrons with the TM ion. Note that the FM (AFM) coupling for valence electrons implies the FM (AFM) coupling for holes as well. Figure 2a and Table 1 shows the central result of this paper, namely the calculated exchange constants N 0 α , N 0 β A , N 0 β E and their average N 0 β = (1/3)(N 0 β A + 2N 0 β E ) for the 3d TM in ZnO in both q = 0 and q = +1 charge states. One can observe the following features that characterize the results:

Discussion
Spin splittings and the exchange-correlation potential in DFT. Before a detailed discussion we present a few general remarks on the s, p-d coupling. When the electron gas is spin polarized, the exchange-correlation potential V xcσ depends on the direction of the electron spin σ . The spin splitting of the CBM is given by where H KS σ is the Kohn-Sham hamiltonian with the appropriate V xcσ . An analogous expression holds for �ε v . To simplify the discussion, we use two approximations. First, we assume that the orbital parts of the CBM wave functions for both spin directions are equal. (Validity of this approximation is discussed below.) Second, for illustrative purposes, we refer to the simple Xα approximation, according to which the spin-dependent exchange potential V xcσ is given by V xcσ (r) = A[n σ (r)] 1/3 , where n σ (r) is the density of electrons with the spin σ and A = −2e 2 (3/4π) 1/327 . With those assumptions, Eq. (4) simplifies to and the potential difference V xc determines the spin splitting of the CBM.
According to Eq. (5), the spin splitting �ε c is given by the product of V xc and the wave function squared. An insight into the mechanisms of the s, p-d coupling is obtained by inspecting those quantities. The first factor, V xc , is a functional of the densities n σ , and stems from the non-vanishing spin density n = n ↑ − n ↓ . It is mainly localized on the d-shell of the TM ion. The final charge and spin densities can be considered as a result achieved in two steps. First, in the absence of coupling between the ZnO host and, e.g., the Mn dopant, the ZnO:Mn system consists in the ZnO host with one vacancy and the Mn atom. The total spin of this system is 5/2. After switching on the coupling, i.e., the p-d hybridization, the Mn spin "spills" onto the ZnO host, mainly onto the first O neighbors of Mn. Thus, the Mn spin is somewhat reduced and the O neighbors become spin polarized, but the total spin is conserved and equal 5/2. The final spin polarization acts as a source of an additional attractive potential for the spin-up electrons, which determines the FM character of the coupling and causes the spin splitting of the band states.
n is shown in Fig. 3a for Co in ZnO. It is dominated by the Co orbitals, and also contains a contribution from the spin polarized O nearest neighbors. Figure 3b shows the difference ( n 1/3 ↑ − n 1/3 ↓ ), to which V xc is approximately proportional. The shape of the isosurface has a tetrahedral symmetry to a good approximation, and reflects that of n . Co 2+ has 7 d-electrons. From Fig. 1a it follows that the two e 2 orbitals are occupied with 4 electrons, two with spin-up and two with spin-down, and their total spin is zero. The spin density shown in Fig. 3 is thus dominated by the three t 2↑ orbitals. Because V xcσ is mainly given by the (1/3) power of n σ , V xc is smoother and somewhat more delocalized than the spin density, which enhances the role of the O neighbors.
The second factor which determines the spin splitting of a given state is its wave function. In a zinc blende crystal, the CBM of the Ŵ 1 symmetry contains the s orbitals only, while the Ŵ 15 VBM states are composed of both www.nature.com/scientificreports/ p and d orbitals. In the wurtzite ZnO these selection rules are relaxed by the hexagonal part of the crystal field, but this small effect can be neglected in the discussion. Those features are seen in Fig. 4a, which presents the wave functions of pure ZnO as a reference. In the case of ZnO:Co, Fig. 4b, hybridization between the TM and the host ZnO states results in the contribution from s(TM) to the CBM, and from d(TM) to the VBM, as it is discussed in detail below. Note that the effect of p-d hybridization on the t 2σ and e 2σ gap states is even more pronounced. It is reflected in their strong delocalization shown in Fig. S5 of Supplementary Information, especially when compared with the compact spin polarization n shown in Fig. 3.
Our results provide a good illustration of the fact that the Heisenberg form of the coupling, H ex = −J S · s , has an "effective" character. Indeed, the Heisenberg hamiltonian suggests that the exchange interaction acts on the spin component of the wave function. Actually, this is not the case, because V xcσ is spin-dependent but it acts on the orbital parts of the wave function. It leads to the differences between the spin-up and spin-down partners of CBM and VBM, and the main difference relies in the different contribution of the d(TM) to the wave functions shown in Fig. S4 of Supplementary Information. s-d exchange coupling. Theory of the exchange coupling between the s-like conduction electrons and the localized f-shell in rare earths (REs) was developed by Liu 3 . He considered a restrained set of states, namely the f-shell of the RE atom and the CBM, and applied the Hartree-Fock approach to obtain the exchange term. This direct (or potential) exchange stems from the overlap of the corresponding wave functions, and is FM. Because of the strong localization of the f orbitals, the overlap integrals extend over the volume of the magnetic atom only, and therefore the s-f coupling is driven by intra-atomic effects. The Liu's picture was then applied to the s-d coupling in DMSs in spite of differences between the two systems, which we discuss in the following.
The tight binding picture allows for an intuitive real space interpretation of the results at the atomic level. The CBM wave functions are represented as sums over the appropriate s R orbitals of atoms at sites R in the supercell: and are normalized to 1 in the supercell volume. The decomposition coefficients a Rσ are obtained by projection of the calculated ψ σ onto the atomic orbitals. Assuming again that the orbital part of the wave function is the same for both spin directions we have an approximate expression where the sum is limited to the first oxygen neighbors of the TM ion, in agreement with the strong localization of both n and V xc , see Fig. 3. In this expression the inter-atomic terms are neglected, which is justified by the smallness of the involved overlap integrals, and the spin splitting of the CBM is the sum of atomic-like contributions. The first one, �ε c (TM) , corresponds to the Liu picture of the s-d coupling, while the second one, �ε c (O NN ) , originates in the spin polarization of the O ions induced by the p-d hybridization with d(TM). where superscript indexes "0" mean the unperturbed d and VBM level energies, and V hop,σ is a spin dependent hopping integral between a TM ion and its neighbors. Average energies of the TM triplet levels are shown in Fig. 2b. Although the Figure presents the final self-consistent energies ε(t 2 ) rather than ε(t 0 2 ) , they can serve as a basis for discussion. Equation (8) together with Fig. 2b qualitatively explain the calculated characteristics of N 0 β in terms of the d-shell energy levels relative to the VBM. For the neutral TM dopants, the first term of Eq. (8) gives a positive (i.e., FM), while the second one gives a negative (AFM) contribution to N 0 β . The first term is dominant, and N 0 β is positive, since t 2↑ is closer to the VBM than t 2↓ . With the increasing atomic number, N 0 β increases due to the decreasing energy denominators, see Fig. 2a. In turn, the t 2↑ level for positively charged Fe, Co and Ni is below the VBM and thus both terms of Eq. (8) lead to AFM exchange coupling. The underlying mechanism is schematically shown in Fig. 5.
In the case of ZnO, the situation is somewhat more complex, since the VBM is a quasi-triplet formed by the A 1 and the E 2 hole subbands. The corresponding TM-induced spin splitting energies are given in Fig. 6a, b. According to our results, �ε A and �ε E can substantially differ, and the difference critically depends on both the Figure 5. Dependence of the p-d exchange coupling on the energies of the t 2σ (TM) levels. The hybridization occurs between the states with the same spin, and is spin-dependent. www.nature.com/scientificreports/ dopant and its charge state. In particular, the spin density of the Mn 2+ with the fully occupied spin-up d-shell is a fully symmetric Ŵ 1 object, it acts on A 1 and E 2 in a very similar fashion, and consequently �ε A and �ε E are almost equal. This is not the case of other TM ions. They are characterized by "non-spherical", i.e., non-Ŵ 1 , spin densities, which act differently on the A 1 and E 2 partners, and induce different �ε A and �ε E . For Cu, even the signs of the splittings are opposite. The impact of the p-d hybridization on the p-d coupling is well illustrated by Fig. 6. The Figure shows the decomposition coefficients a 2 TMσ of the VBM wave functions in the tight binding picture analogous to Eq. (6). By comparing Fig. 6c,d with g-h one observes that the hybridization is strongly spin-dependent, since the contribution of the spin-up and spin-down TM states to the VBM can differ by as much as one order of magnitude. This stems from the different energies of the TM gap states relative to the VBM, i.e., the different energy denominators in Eq. (8), which also control the mixing of wave functions. In most cases, the spin-down states are more distant from the VBM than the spin-up ones, and the contribution of d(TM) to the VBM is appreciably larger for the spin-up than for the spin-down channel. However, since the spin splitting is the energy difference, we have to consider the appropriate combinations of the a 2 TM↑ and a 2 TM↓ coefficients rather than their values separately. Depending on the actual level ordering, the combination is ±a 2 TM↑ − a 2 TM↓ , consistently with the Eq. (8) , where the +(−) sign holds when the spin-up TM level is above (below) the VBM. The results are shown in Fig. 6e,f. The very high level of correlation between the splittings �ε A and �ε E and the contribution of the d(TM) orbitals to the VBM is clear. We also note that because of the large differences between a 2 TM↑ and a 2 TM↓ , the approximate Eq. (7) cannot be applied to the VBM, and the appealing separation into the Liu intra-atomic contribution and the hybridization contribution does not hold. However, the intra-atomic contribution of d(TM) to �ε v does not vanish. Since the direct exchange leads to the FM coupling, it enhances the hybridization-induced values of the N 0 β s for all TM 2+ dopants, and reduces the negative values of N 0 β s for Fe 3+ , Co 3+ and Ni 3+ .
In a complementary approach, one can follow the impact of the p-d hybridization in the real space picture. The relevant wave functions for ZnO:Co are given in Fig. 4. In agreement with the attractive character of V xc and with Fig. 6, the spin-up VBM wave function contains a larger contribution of d(Co) than its spin-down partner. Indeed, from Fig. 4b it follows that the approximation of Eq. (5) is not justified for the VBM, because the contribution of Co to the spin-up (8.8%) and spin-down states (0.3%) differ over 20 times. The contributions of the d(TM) orbitals to the VBM for other TM dopants is shown in Fig. S4 of Supplementary Information. They depend on the dopant and its charge state.
Analyzing the consecutive TM ions we find that, as it follows from Fig. 6, in the case of the light V the contribution from the d(TM) shell to the VBM practically vanishes. Accordingly, the VBM spin splittings and the corresponding N 0 β s are small. Next, the average values of N 0 β increase to 0.4 -0.6 eV for Cr and Mn with the decreasing energies of t 2 . In the sequence Fe -Ni, t 2↑ is very close the valence band, which strongly enhances both the spin splitting of the VBM and N 0 β . Finally, Cu represents an interesting case, since the splittings �ε A and �ε E have opposite signs. This effect takes place because some of the energies of the d(Cu) states are below the VBM, thus changing the sign of the splitting in agreement with Fig. 5.
In the case of charged dopants with q = +1 , the TM-induced levels are lower in energy than for q = 0 , see Fig. 2b. This leads to higher values of N 0 β for Cr and Mn. More importantly, in the case of Fe, Co, and Ni, the t 2↑ levels are below the VBM, which changes the sign of the first term in Eq. (8), and drives the change of character of the p-d coupling to AFM, as displayed by the negative N 0 β in Fig. 2a. In particular, the AFM coupling can be expected for Fe, which typically assumes the Fe 3+ charge state 18 . Comparison with previous calculations and with experiment. The role of the p-d hybridization in the p-d exchange coupling was recognized early. Typically, it was taken into account by using the Anderson hamiltonian 4 . The calculations employing the Schreiffer-Wolf 5 transformation were performed and found a negative N 0 β for Mn in II-VI semiconductors like CdTe 6,11 . Indeed, it was properly recognized that the d(Mn) states in CdTe are placed about 2 eV below the VBM. As it is shown in Fig. 5, in this situation the p-d exchange www.nature.com/scientificreports/ coupling of holes with the Mn 2+ spins must be AFM independent of the details of calculations. This approach was used for other TM ions leading to AFM coupling for Mn 2+ , Fe 2+ , Co 2+7,9,28 and FM coupling for Sc 2+ and Ti 2+9, 28 . The configuration interaction and cluster-model calculations were also used to evaluate N 0 β for several II-VI hosts and TM dopants, and to interpret the experimental data; the obtained N 0 β 's were negative indicating the AFM coupling with holes [29][30][31][32] . In the papers above, simple expressions for N 0 β are given, in which the critical factor is the energy of the majority spin d(TM) level relative to the VBM. In the case of ZnO, the TM level energies were assumed to be similar to those in CdTe, which leads to AFM coupling with N 0 β of about −3 eV for Mn 2+ , Fe 2+ and Co 2+8,10,32 . However, the assumption that t 2↑ (Mn) is below the VBM in ZnO was invalidated by experiment. The measurements 33 proved that t 2↑ (Mn) is in the band gap, which was subsequently confirmed by both experiment and theory 17,25,34 . Consequently, in this case the p-d coupling is FM. This should be the case of Cr as well, since t 2↑ (Cr) is higher than that of Mn. Similarly, the energy of the gap levels of Co 2+ are well established 19,25 , and N 0 β > 0 is expected. This may imply that the interpretation of the X-ray data, leading to the negative N 0 β for Mn 31,32 , and possibly for Fe, Co, and Ni, was not correct. Finally, as observed in Ref. 33 , "location of the d levels below the VBM is an essential assumption behind the proposal of mid-gap Zhang-Rice-like states in ZnO:Mn 35 ". According to the results presented above this assumption is not correct.
The same issue, i.e., correctness of energies of the t 2 (TM) levels relative to the CBM and VBM, is present in the case of calculations based on the local density approximation. This approximation results in a severe underestimation of the ZnO band gap, and therefore wrong energies of the TM levels and the non-correct sign of the p-d coupling (e.g., N 0 β = −1.81 eV for Mn in ZnO 14 ). On the other hand, when the correct E gap of ZnO is used 25 , the energies of the TM gap states are close to the present results.
We now turn to the experimental results. In early magneto-optical measurements, a strong exchange interaction between band carriers and localized d(TM) electrons was reported for Mn, Fe, Co, Ni and Cu [36][37][38][39] . By contrast, the s, p-d coupling was not observed so far for Ti, V and Cr ions in ZnO 37 . This latter result is in a reasonable agreement with our findings. Indeed, the stable charge state of Ti is the nonmagnetic Ti 4+ ( q = 2 ) as in experiment [20][21][22] , and thus both exchange constants, N 0 α and N 0 β , vanish by definition. Also, V (stable in the V 3+ charge state according to both experiment 23,24 and our calculations), as well as Cr, are characterized by small N 0 β exchange constants.
In the case of Co, magnetooptical measurements showed that N 0 |β − α| = 0.8 eV 43 . The sign of N 0 β could not be determined experimentally due to the ambiguity in the valence bands ordering. By assuming N 0 α = 0.25 eV, it was concluded that the p-d coupling can be either FM with N 0 β = 1.0 eV, or AFM with N 0 β = −0.55 eV. Our calculated N 0 β depends on the charge state. For the neutral Co 2+ we find the FM coupling ( N 0 β γ = 4.2 and 2.6 eV for the A 1 and E 2 subbands, respectively), while for the positively charged Co 3+ the coupling is AFM ( N 0 β γ = −0.9 and −0.6 eV for the A 1 and E 2 subbands, respectively). Thus, we obtain that N 0 |β − α| = 2.6 eV for Co 2+ , and about 1.2 eV for Co 3+ . The latter value is reasonably close to magnetooptical data 43 .
N 0 β s determined for Fe and Ni are large and negative, namely − 2.7 eV 32 and −4.5 ± 0.6 eV 44 , respectively. These values are somewhat higher than our results for q = 0 ( N 0 β E = 2.0 and N 0 β A = 0.8 eV for Fe, and N 0 β E = 4.2 and N 0 β A = 2.3 eV for Ni), but importantly they are of opposite sign. However, we predict negative but smaller values for those dopants in the q = +1 charge states. Here, it should be mentioned that large and negative N 0 β s were proposed also for Mn (-3.0 eV 32 , -2.7 44 ) and for Co (-3.4 eV 32 , -2.3 44 ). However, as it is pointed out above, interpretation of these measurements is based on particular assumptions regarding the energies of the TM levels, which were subsequently questioned. We conclude that a reliable comparison with experiment for Fe, Co, and Ni requires the charge state of the TM ion to be established.

Summary and conclusions
Theoretical analysis of the s, p-d exchange coupling between free carriers and the 3d transition metal dopants in ZnO was conducted employing the GGA +U method. The present study reveals both the detailed characteristics for each ion, and general trends. A particular care was devoted to reproduce the correct band gap of ZnO and energies of the gap levels of the dopants. The calculated s-d coupling constant N 0 α is about 0.5 eV for all the TM ions, i.e., it does not depend on the dopant and its charge state. By contrast, the p-d exchange coupling reveals unexpectedly complex features. First, N 0 β s strongly depend on the chemical identity of the dopant, increasing about 10 times from V to Cu. Second, N 0 β is different for the two VBM subbands, the light hole A 1 and the heavy hole E 2 , and the corresponding values can differ by a factor 2, or even have opposite signs. Third, not only the magnitude but also the sign of N 0 β depends on the charge state of the TM ion. In particular, the coupling between holes and Fe, Co and Ni ions in the q = 0 neutral charge state is strong and ferromagnetic, while for q = +1 the coupling changes the character to antiferromagnetic. Finally, the stable charge state of Ti in ZnO is Ti 4+ , in which its spin vanishes.
Analysis of the wave functions reveals how the hybridization between the TM orbitals and the ZnO band states determines the s, p-d exchange coupling. The magnitude of the p-d coupling is determined by the energies of the d(TM) relative to the VBM. The most striking example is that of Cu, for which the t 2↑ (Cu) and VBM are almost degenerate, and the actual ordering of the d(Cu)-induced levels and the VBM explains different signs of N 0 β s of the light and the heavy holes, and their large magnitudes. Thus, the p-d(TM) hybridization leads to the Anderson-like picture of the p-d coupling, but its role is more complex. In particular, the p-d hybridization www.nature.com/scientificreports/ affects not only N 0 β but also the N 0 α constant. The main mechanism of the s-d coupling is grasped by the Liu's model, and it originates in the TM intra-atomic exchange interaction between the s and d electrons. However, the spin polarization of the oxygen neighbors of the TM ion induced by the p-d hybridization leads to the spin polarization of the s(O) orbitals, which contributes about 1/3 to the N 0 α constant. Comparison with experiment is satisfactory for Ti, V, Cr and Mn. In the case of Fe and Co, the definitive conclusions are not possible, because N 0 β depends on the dopant charge state, which was not assessed in experiment. An acceptable agreement for Co is obtained assuming the Co 3+ and not the Co 2+ charge state.

Method of calculations
The calculations are performed within the density functional theory 45,46 in the generalized gradient approximation (GGA) of the exchange-correlation potential V xc 47 , supplemented by the +U corrections 16 . We use the pseudopotential method implemented in the quantum espresso code 48 , and employ ultrasoft pseudopotentials, which include nonlinear core correction in the case of Co and Ni. The valence atomic configuration is 3d 10 4s 2 for Zn, 2s 2 p 4 for O, and 3s 2 p 6 4s 2 p 0 3d n or 4s 2 p 0 3d n for TM ions with n electrons on the d shell. For V, Ti, Cr, Mn, Fe and Co, the plane-waves kinetic energy cutoffs of 30 Ry for wave functions and 180 Ry for charge density are employed. Convergence was assessed by test calculations with cutoffs of 40 Ry. Following the recommendation of quantum espresso, for Ni and Cu the cutoff is increased to 45 Ry.
Spin-orbit interaction is neglected. We justify this approximation by the results of experiments regarding TM dopants in ZnO. The interaction manifests itself in optical measurements, where the lines of intracenter d − d transitions reveal rich structures, and are split in particular by the spin-orbit coupling. According to the results for Co 49,50 , the spin-orbit splittings of the initial and final states of the T 4 2 (F) -A 2 (F) emission can be estimated as 19 and 6 cm −1 , respectively, which corresponds to 1-3 meV. Similar values, lower than 10 meV, were reported for other TM dopants 50,51 . In the case of Mn 2+ and Fe 3+ , due to the absence of orbital momentum in the d 5 configuration, the second order spin-orbit interaction leads to splitting energies below 1 meV 51 . Moreover, the spin-orbit splitting of the VBM in ZnO is about 10 meV, i.e., it is very small compared to the spin splittings of the order of 1 eV characterizing heavier atoms and semiconductors such as InSb, CdTe or PbTe. The smallness of the spin-orbit coupling in ZnO:TM justifies its neglect both in our [17][18][19] and in the previous ab initio calculations 25,26 .
The electronic structure of the wurtzite ZnO is examined with an 8 × 8 × 8 k-point grid. Analysis of a single TM impurity in ZnO is performed using 3 × 3 × 2 supercells with 72 atoms, while k-space summations are performed with a 3 × 3 × 3 k-point grid. For too small supercells, the spurious defect-defect coupling can distort final results. Convergence of the results with respect to the supercell size was checked for Mn and Co with 6 × 6 × 4 supercells with 576 atoms. We obtained that the N 0 α s are the same for both supercells, while N 0 β s are lower by 10 [17][18][19]60 . For other TM dopants, the U corrections are taken to be 2-3 eV, as suggested in the literature 25,26 . We checked that in most cases a variation of the U(TM) value by 1 eV alters the impurity levels by about 0.1 eV, and the N 0 α or N 0 β values by less than 0.05 eV. We also mention that in spite of the energetic proximity and strong hybridization between the VBM and the TM-induced levels, a non-ambiguous identification of hole states was always possible based on the analysis of wave functions. However, since the p-d coupling depends on the inverse energy distance between the TM-induced levels and the VBM, the results are less accurate for Co, Ni and Cu than for Cr and Mn.
Various charge states of the TM dopants are considered. In general, in the absence of additional dopants, a TM ion occurs in the neutral q = 0 charge state, denoted as TM 2+ , and other charge states q can also be assumed when defects are present. Generally, the stable charge state of a defect in a semiconductor depends on the Fermi level. Transition level ε(q/q ′ ) of a defect is defined as the Fermi energy at which the stable charge state changes from q to q ′ , or in other words, as the Fermi energy at which formation energies of q and q ′ are equal: where E(q) is the total energy of the doped supercell and ε v is the VBM energy of pure ZnO. The finite size effects are taken into account by including the image charge corrections and potential alignment for charged defects 61,62 . The energies ε(q/q ′ ) in the gap determine the possible charge states of TM ions.
Finally, the spin splitting energies of the VBM and CBM are taken directly from the Kohn-Sham levels. Alternatively, they can be obtained from appropriate excitation energies, as discussed in Supplementary Information.

Data availibility
The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.