Spin Polarization and Magnetic Properties of VGaON and VGaONInGa in GaN: GGA+U Approach

Electronic structure of a defect center containing the gallium vacancy and substitutional oxygen atom at nitrogen site (VGaON) in zinc blende and wurtzite GaN was analyzed within GGA+U approach. The +U term was applied to d(Ga), p(N), p(O), and d(In). Neutral VGaON is in the stable high spin state with spin S = 1. The defect structure is strongly dependent on geometry of the defect and the charge state. Two spin structures, which arise due to two different configurations in VGaON, with ON either along the c-axis or in one of three equivalent tetrahedral positions in wurtzite structure were analyzed. The weak ferromagnetic coupling between centers was found. The strength of magnetic coupling is increased when there is a complex containing VGaON with additional substitutional indium atom at the second neighbor to vacancy gallium site (VGaONInGa). Magnetic coupling between VGaONInGa is ferromagnetic due to strong spin polarization of p electrons of the nearest and distant nitrogen atoms.

III-nitride materials such as GaN and InGaN have found their applications in advanced solid-state lighting technologies [1][2][3] and optoelectronic devices including diodes or solar cells [4][5][6]. Moreover, the ferromagnetism (FM) in GaN or InN without doping by transition metal atoms was recently observed [7][8][9][10][11]. This FM was ascribed to the formation of native defects, such as cation vacancies or their complexes. For example in Ref. [9], analysis of the characteristics of hysteresis curves in irradiated GaN showed that the coercive field increases in line with the increase of concentration of gallium vacancy (V Ga ). The stable FM was discovered in n-type GaN, and it can be due to the presence of non-intentional donors, such as oxygen.
Measurements using positron annihilation spectroscopy of V Ga and its complexes containing the gallium vacancy and substitutional oxygen atom at nitrogen site V Ga O N have found them to be dominating defects in as-grown n-GaN [12][13][14][15][16][17][18][19][20][21]. Both defects were intensively studied in experiments for their optical properties [12,13,16,17,19,21] indicating defects as possible sources of the green (GL) and yellow (YL) luminescence in GaN. In Ref. [12], Son et al. detect spinpolarized V Ga O N by electron paramagnetic resonance and suggest that YL and GL bands can be explained by 0/− and −/−2 V Ga O N -optical transition levels. Moreover, in Ref. [13], two electronic structures were observed, which arise due to two different configurations of V Ga O N in wurtzite (w) GaN, one with O N either along the c-axis (axial configuration, referred below as a) and other one in one of three equivalent tetrahedral positions (basal configuration, referred below as b).
The increasing localization of the defect wave function has the opposite effect on the stability of the local magnetic moment and on the collective magnetization: in the former case, the increased localization stabilizes the high-spin (HS) state, while the coupling through the overlap of the wave functions of the neighboring defects is decreased [27,30,31,35]. For both the V Ga and the V Ga O N complexes, strong localization can lead to stable local spin moments [26,27,35] but it does not guarantee automatically a stable interaction between them [31,35]. Partial delocalization of defect-induced bands may reduce the stability of HS state defect but also be responsible for the long-range magnetic interactions. This stabilization can be due to p-d exchange interaction of impurity like Mn [36] or the spin polarization of p electrons in low In-content InGaN.
Both In concentration and microscopic In distribution strongly influence the electronic structure and physical properties of InGaN [37]. The localization of the valence band maximum (VBM) states and the domination of the light emission of InGaN with low In content were observed [37]. Similar to GaN, V Ga -complexes were also suggested to be the important non-radiative defects in InGaN quantum well [38,39]. Strong effects in electronic structure of V Ga O N -hydrogen complex were found in p-type InGaN for high In content [33].
To check the hypothesis that magnetic coupling between V Ga O N -complexes can be more stable in InGaN, in the present paper, we study V Ga O N -structure in zb (zinc blende)-and w-GaN and InGaN using GGA+U approach. After presenting the details of calculations in Sect. 1, the justification of the chosen approach is underlined. Next, we present results of calculations of the formation energies of defects (Sect.

Details of Calculations
Calculations based on the density-functional theory were performed using the ultrasoft pseudopotentials [40], the Perdew-Burke-Ernzerhof GGA exchange-correlation potential [41], including the +U term implemented in the QUANTUM-ESPRESSO code [42] along the theoretical framework developed in Ref. [43]. We employed ultrasoft atomic pseudopotentials and chose 3d, 4s, and 4p orbitals for Ga; 4d, 5s, and 5p for In; 2s and 2p for N; and O as valence orbitals. The plane wave basis with the kinetic energy cutoff (E cut ) of 40 Ry provided a convergent description of the analyzed properties. The Brillouin zone summations were performed using the Monkhorst-Pack scheme with a 2 × 2 × 2 k-point mesh [44]. Methfessel-Paxton smearing method with the smearing width of 0.136 eV was employed for obtaining partial occupancies. The zb 216-and 512-atoms and w 128-, 192-atoms supercells were considered, and ionic positions were optimized until the forces acting on ions were smaller than 0.02 eV/Å. The spin-orbit coupling was neglected. Formation energies were calculated according to Ref. [45] for N-rich conditions [30]. The +U corrections were imposed on d(Ga) and p(N) [30,31]   The band gaps E gap calculated within LDA/GGA are shrunk and amount only to about 1.4 and 1.6 eV for the zb-and wcrystals, respectively (see Fig. 1 and Refs. [30,31]). Moreover, both LDA and GGA give too high calculated  [46]. Inclusion of large U(Ga) = 10 eV term solves the problem with the correct position of d(Ga) [46][47][48], but E gap is still about 2.5 eV (see Fig. 1). Our previous works showed that increasing U(N) from 0 to 5 eVopens the band gap of w-GaN from 1.6 to 3.0 eV [30,31], but with d(Ga) centered about 13 eV below the VBM, similar to LDA/GGA results [30,46]. Here, we systematically analyze how the U(Ga) and U(N) terms affect the electronic structure in both the zb-and w-crystals. On-site +U parameter, varying from 0 to 10 eV, was applied separately on d(Ga) and p(N), and then together on both d(Ga) and p(N) orbitals, where U(N) = 5 eV. The energy E gap calculated as a function of U is shown in Fig. 1.
We found that U(Ga) = 3.0 eV along with U(N) = 5 eV reproduces the experimental E gap of 3.2 and 3.4 eV for both zband w-GaN [49] (Fig. 1, (3, 3′)) and the binding energy of Ga 3d level centered about 15.5 eV below the VBM-in agreement with Ref. [50]. These values are also in agreement with HyF results [51]. Such an underestimation of the band gap and band structure follows from the sublinear dependence of the LDA/GGA total energy on the occupation [43]. Moreover, the sensitivity of E gap in both U(Ga) and U(N) is explained by the orbital compositions of both the VBM and the minimum of the conduction band (CBM).
GGA+U calculations give the lattice constants a zb = 4.57 Å and (a w = 3.19 and c w = 5.2 Å) for zb-and w-GaN, respectively. These values are very close to the experimental data of a-

GGA+U vs HSE Results
It was noted above that HyF [16-18, 21, 22, 26-28, 32, 33] calculated V Ga -structure is in agreement with GGA+U results for U(N) = 5 eV [30,31]. Here, in order to verify the agreement further, we perform calculations for Heyd, Scusseria, and Ernzherof functional, based on the PBE functional where parameter α is a fraction of the exchange that is replaced by Hatree-Fook exchange [54]. Calculations were done for isolated V Ga O N and V Ga O N -V Ga O N (3NNs axial configuration in notation of Sect. 2.4) for w-GaN. Screening parameter α = 0.25 was set to reproduce band gap of~3.0 eV. Nevertheless, the results indicate that energies of spin polarization (ΔE PM −FM ) (defined in Sect. 2.1) agree to within 0.05 eVor less, and energies of magnetization (ΔE AFM−FM ) (defined in Sect. 2.2) agree to within 0.005 or less. That shows the good agreement for calculations of magnetic properties within these two approaches.
We note that the problem of choosing the α parameter for getting accurate defect levels is still an open issue [55,56], as well choosing the U parameter in GGA+U approach [30].

Formation Energy of Defects
Formation energy of charged V Ga O N was calculated. One geometry of a-V Ga O N was considered for cubic GaN and two aand b-V Ga O N configurations were analyzed for w-GaN. Because the aim was to understand the influence of In doping on spin-polarized properties, the number of configurations of complex V Ga O N was chosen: In Ga as second nearest neighbor to V Ga when forming O N -In Ga (referred below as o) or N-In Ga (referred later as n) chains where these O N and N are the nearest possible positions to V Ga neighbors. Hence, the two geometries as referred as a-o-and a-n-V Ga O N In Ga were considered for zb-crystal. And four geometry configurations referred as a-o-, a-n-, b-o-, and b-n-V Ga O N In Ga were analyzed for w-GaN. In this work, the formation enthalpy of InGaN was calculated and the binding energy, E bind , defined as a difference in the total energies of compounds that contain V G a O N In G a or (V G a O N and In G a ) were taken into consideration.
Formation energy E form of a defect was calculated by formula taken from Ref. [45] The first two terms on the right-hand side are the total energies of the supercell with and without the complex, respectively. n i is number with the +(−) sign corresponding to the removal (addition) of atoms. E VBM is the energy of the VBM of bulk GaN, and ε F is the Fermi energy referenced to this E VBM . The energy E VBM is determined from the total energy difference between the pure crystal with and without a hole at the VBM in the dilute limit by algorithm from Ref. [45]. μ i are the variable chemical potentials of atoms in the solid, which in general are different from the chemical potentials μ i (bulk) of the ground state of elements (Ga bulk, In bulk, and N 2 , O 2 ). Chemical potentials of the components in the standard phase are given by total energies per atom of the elemental solids: μ(Ga bulk) = E tot (Ga bulk), μ(In bulk) = E tot (In bulk), while μ(N bulk) = E tot (N 2 )/2 and μ(O bulk) = E tot (O 2 )/2 (ΔH f (NO) was neglected. In N-rich condition, μ(Ga) = E tot (Ga bulk) + ΔH f (GaN) and μ(In) = E tot (In bulk) + ΔH f (InN) are taken, where ΔH f is the enthalpy of formation per formula unit, and it is negative for stable compounds. ΔH f at T = 0 K is obtained by considering the reaction to form or decompose a crystalline GaN and InN from its components and dependent on an cohesive energy, E coh , of Ga, In, N, and O. The obtained results for E coh of Ga, N, and O were shown in Refs. [30,57]. Calculated E coh (In) and ΔH f (zb-GaN), ΔH f (InN) are 2.56 (2.5 [58]) and − 1.24 (−1.27), −0.36 (−0.32) eV [59]), (experimental values presented in brackets).
The last term, E correct , includes two corrections. The first one, ΔE PA , is the potential alignment correction of the VBM. The VBM in the ideal supercell and in the supercell with a (charged) defect differs by the electrostatic potential and is obtained by comparing the potential at two reference points far from the defect in the respective supercells with (P[D q ]) and without (P[0]) the defect, Second correction is an image charge correction as expressed by 2-order Makov-Payne form: where α M is the lattice-dependent Madelung constant, which for hexagonal structure is 3.5, W is the supercell volume, and ε is the static dielectric constant. E MP was calculated to be 0.2, 0.4 eV for charged defects (q = − 1, − 2). Results of calculations are presented in Sect. 2.2.

Results
This section summarizes the obtained results for the defect structure and formation energy and discusses magnetic interaction between defects.

Electronic Structure and Spin Polarization of V Ga O N
The defect states of V Ga O N stem from the result of the interaction between the vacancy orbitals and the O-impurity states. Local atomic configuration of this defect has the C 3v point symmetry. In both the zb-and the w-structures, the vacancy is tetrahedrally surrounded by three N and one O atoms, and the respective defect states are localized on the resultant broken sp 3 bonds (Fig. 2) that split into a singlet a g and a higher in energy quasitriplet "t 2 ". Energies of a g and "t 2 ", calculated below as relative to the VBM, depend on the crystal structure, geometry of defect and the charge state. a g is a resonance state with the valence bands."t 2 "-V Ga O N is located in the band gap and splits into a doublet e 2 and a singlet a 1 with a splitting energy of about 0.2, 0.5, and 0.7 eV for the zb-GaN, a-and bgeometries of V Ga O N in w-crystals, respectively (see Fig. 3, the green lines). The energy splitting contains the contributions of a weak perturbation by the w-symmetry in w-GaN and the U-induced so-called quasi-Jahn-Teller (JT) effect [30].
In the case of non-vanishing spin polarization, the exchange coupling splits e 2 into spin-up e 2↑ and spin-down e 2↓ states by the splitting exchange energy defined as Δε ex = ε(e 2↓ ) − ε(e 2↑ ), where ε is the energy of the defect level, and a 1 into a 1↑ and a 1↓ (Fig. 3, the blue lines). The Δε ex , in general, depends on the symmetry of defect, the charge state and U. The e 2↓ and a 1↓ of neutral V Ga O N in zb-GaN are localized in the band gap at 2.6 and 1.9 eV above the VBM, respectively. The e 2↑ and a 1↑ are resonances with the VBM (Fig.  3a). According to this point of view V Ga O N is a deep acceptor containing two holes.
As presented above in Section 2.3, in w-crystal, the defect can exists in two different geometries tagged as a-and b-V Ga O N . According to GGA+U calculations, single-electron level representations of electronic structures are different for a-and b-V Ga O N (Fig. 3 b, c, left panels). Structure of a-V Ga O N in w-GaN is similar to the one of zb-GaN, e 2↓ and a 1↓ are 2.45 and 1.65 eV with the respect to the VBM, and e 2↑ and a 1↑ are hybridized with the valence bands (Fig. 3b). Introducing O atom into basal plane of defect leads to strong symmetry perturbation, e 2↓ is split by 0.5 eV into two a 2↓(1) and a 2↓(2) singlet states. Thus, "t 2↓ " in this case is a composite band containing a 1↓ , a 2↓(1) , and a 2↓(2) levels located about 1.4, 2.3, and 2.8 eV above the VBM, respectively (Fig. 3c).
The energies of V Ga O N levels strongly depend on the charge state q. Single-electron energy levels of V Ga O N q for q = 0, − 1, and − 2 with their respective charge states are shown in Fig. 3. The physics behind the calculated electronic structure of charged defects is determined by the following counteracting effects [30]: (i) the intracenter Coulomb repulsion is dominant in the nonspin-polarized calculations. Without spin polarization, the energy of e 2 , a 1 increases by~0.5-0.6 eV with the q changing from 0 to − 2 (the levels are shown in green color in Fig. 3); (ii) The effect of the value exchange splitting, for example, in zb-GaN with q changing from 0 to − 2, Δε ex decreases from 2.6 to 0 eV; (iii) the U-induced potential which is attractive (repulsive) for occupied (unoccupied) orbitals [43]. This effect is clearly seen for e 2↓ that decreases with q changing from 0 to − 1 (Figs.3 a, b). These results are in agreement with HyF calculations from Ref. [17] where negative-U eff behavior of V Ga O N was observed.
"t 2 " defect state is occupied by 4, 5, and 6 electrons for q = 0, − 1, and − 2, respectively. As shown in Figs Table 1. Generally, ΔE PM−FM assumes the maximal value when "t 2 " is occupied with 4 electrons with the spin S = 1, and it vanishes when "t 2 " is fully occupied. Moreover, the energy of antiferromagnetic state (AFM) of single defect was considered also in the analysis [23]. Every considered geometry of the neutral V Ga O N is the magnetic centrum in HS state with the local magnetic moment  (Table 1). Finally, V Ga O N −2 is the nonspin-polarized centrum (in this case ΔE PM−FM = 0). According to our calculations, b-configuration stabilizes the The geometry of the defect affects also the localization of the V Ga O N states. The effect is displayed in the plots of the density of spin polarization, Fig. 2a,e-f. Figure 2 indicates that the V Ga O N states are dominated in all cases by the localized and spin-polarized contribution of the three sp 3 orbitals of the N nearest neighbors because O-atom is more electronegative than N-atom and two electrons with the opposite spins are located on sp 3 oxygen orbital. In contrast to GGA method in which electrons that occupy for example, "t 2↓ " are spread over four p orbitals of the nearest neighbors [31,60], GGA+U calculations showed that the partial occupancy is avoided [31,60]. Moreover, in the case of b-V Ga O N , one can observe the anisotropy for three sp 3 of N orbitals (Fig. 2d, f). The Uinduced symmetry breaking of e 2 level, i.e., the quasi-Jahn-Teller effect, was observed in b-V Ga O N . In a-V Ga O N contributions of the three N neighbors to the V Ga states of the vacancy, wave function are almost equal, whereas for b-V Ga O N the wave function is dominated by the two basal N ions located in the (x,y) plane, and the contribution of the remaining N ion is strongly reduced, see Fig. 2d-f.

Formation Energy and Transition Levels
The calculated E form of V Ga O N and assumed ε F at the VBM are given in Table 2. E form in all configurations is the same. The same trend was observed also in Ref. [32]. The obtained values 1.64, 3.54, and 5.73 eV in w-GaN are close to results of HSE (1.9, -, and 5.3 eV [17]) and HyF (B97-2-functional) (1.5, 4.0, and 6.2 eV) [32]. Our calculations demonstrate that E form is similar in all geometries.
The change in the defect charge state is determined by the transition level ε(q 1 /q 2 ), defined as the Fermi energy relative to the VBM at which formation energies of the q 1 and q 2 charge states are equal. We find ε(0/−) are 1.84 and 1.9, 2.1 eV for defects in zb-GaN, and a-, b-geometry in w-GaN, respectively, and ε(−/−2) are 2.1 and 2.2, 2.4 eV for defects in zb-GaN, and a-, b-geometry in w-GaN, respectively, which is consistent with V Ga O N energies shown in Fig. 3. Comparison of the calculated energies with the results for other exchange functionals is shown in Table 3. Table 3 contains dataset for only a-V Ga O N . Comparable values of ε(0/−) and ε(−/−2) were obtained with GGA+U and different HyF exchangecorrelation functionals. For example, our results are close to the obtained by HSE06 approach (with 20% exact exchange) [17,18].

Spin Polarization of V Ga O N In Ga
In Sect. 1.3 the geometrical configurations of V Ga O N In Ga were discussed. Although our calculations show that binding energy in such a complex is low,~0.1-0.3 eV, the formation energy is 1.9 eV which is a little higher than in the case of V Ga O N and considerably lower than in the case of V Ga [30,34]. Because the formation energy for different geometries (aor b-) is similar, both V Ga O N and V Ga O N In Ga can exist in nonequivalent atomic configurations.
In this section, the results for calculated energy of spin polarization, local magnetic moment, spin density, and density of charge of V Ga O N In Ga are presented. Moreover, the effect of addition of In impurity on the spin-polarized properties is analyzed by comparing these results with similar results for V Ga O N . Next, in order to get a clearer picture of the influence of In on electronic structure, the complex V Ga O N nIn Ga , where n = 2, 3, and 4, is investigated in zb-GaN.
The defect states of V Ga O N In Ga stem from the result of the interaction between the vacancy orbitals and the O-and Inimpurity states.    Table 2 The calculated E form of V Ga O N assuming ε F at the VBM (N-rich condition) Strong effect on the electron structure was observed for wurtzite crystal (see Fig. 3b, c, right panels). As with Oatom in basal plane, also in this case for both a-and b-V Ga O N In Ga , "t 2 " is split into three singlet states a 1↓ , a 2↓(1) , and a 2↓(2) due to strong tetragonal perturbation generated by In atom in the crystal structure. a 2↓(2) of b-n-V Ga O N In Ga is higher in energy, and it is a resonance state with the conduction bands (Fig. 3c).
Calculated ΔE PM−FM values for different configurations of V Ga O N In Ga are given in Table 1. According to our calculations, the inclusion of In Ga stabilizes HS state of the defect by 0.1-0.16 eV in comparison to V Ga O N (see Table 1). For example, in zb-GaN, the differences in the spin polarization energies of V Ga O N and V Ga O N In Ga are 0.08 and 0.14 eV, for q = 0 and − 1, respectively. The same energy gain takes place in w-structure, for example, the differences between ΔE PM −FM (b-n-V Ga O N In Ga ) and (b-V Ga O N ) are 0.16 and 0.08 eV, in cases q = 0 and − 1, respectively. Generally, the increase of ΔE PM−FM is in agreement with the changes of electron structure. Shifting up the single-electron levels (a 2↓ (2) in particular) to the CBM leads to the increase in the energy of exchange splitting, and therefore in the energy of ΔE PM−FM . The above results obtained for V Ga O N In Ga are typical: the stability of the HS state predicted for V Ga O N In Ga is 10-12% higher than that predicted for V Ga O N.
The increase in discrepancy of ΔE PM−FM was observed with the rising of In content. ΔE PM−FM is 1.7, 1.78, 1.88, 1.9, and 1.8 for n = 0, 1, 2, 3, and 4. For n = 4, ΔE PM−FM is smaller than for n = 2 and n = 3.
The spin (Fig. 2b, c, h-j) and electron (Fig. 2k, l) contour densities of V Ga O N In Ga are shown in Fig. 2. Like in the case of V Ga O N , V Ga O N In Ga states are localized mainly on the p orbitals of the three N atoms (Figs. 2 b, c, h-j). Spin density of V Ga O N In Ga is more delocalized, since they comprise longrange tails that involve p orbitals of distant N ions or even p orbital of oxygen (refer to V Ga O N 2In Ga , see Fig. 2c). These orbitals constitute the VBM (Fig. 4) of GaN. V Ga O N In Ga -spinup states are resonances degenerate with the valence bands, i.e., both a 1↑ and e 2↑ hybridize with the upper part of the VBM, which forces also a partial delocalization of their wave functions. It is evident from Fig. 2 that enhancing correlation effects by addition of In atom leads to the decrease of the localization of the spin density on the broken bonds and to the increase of its axial anisotropy.
The contour (Fig. 2k, l) is plotted in the (100) plane and shows that a large contribution to the electron density comes from the nitrogen atoms and strong ionic bonds resulting from sp 3 hybridization. The spherical symmetry around anions is observed, which indicates that the bonds in GaN:V Ga O N In Ga are dominated by the ionic component.
By calculating the contributions of individual atoms projected onto relevant atomic orbitals to the total DOS (Fig. 4a), one finds that the main contribution comes from the defects states p(N) of the N nearest neighbors of V Ga (Fig. 4b)-in agreement with Fig. 2. Both p(N) and p(O) orbitals also build the VBM of GaN (see Fig. 4b-d). The contribution of the d(In) orbitals to the spin density is non-negligible due to the substantial contribution of d(In) to the VBM and defect states (see Figs. 2 and 4f).

Magnetic Interaction of
In the present section, we study V Ga O N -V Ga O N and V Ga O N In Ga -V Ga O N In Ga defect pairs and analyze the impact of crystal distortion on their properties by comparing the results of magnetic interaction calculations. The magnetic coupling between vacancy complexes is discussed as the possible origin of the experimentally observed FM in GaN. Electronic structure of a defect pair is determined by three factors: (i) the distance between vacancy complex, (ii) the relative orientations of complexes with respect to each other and to the crystal axes, and (iii) the charge state. V Ga O N -V Ga O N and V Ga O N In Ga -V Ga O N In Ga configurations were considered, in which the defects are the third nearest neighbors (3NNs) and fourth nearest neighbors (4NNs) with respective spatial separation of about 5.2 and~6.5 Å (it is a distance between V Ga in the relaxed structure). In w-GaN, the defects can be located either in the same (x,y) basal plane perpendicular to the c-axis, which is referred to as the xy-case, or they can be oriented along the caxis, which is denoted here as the c-case. Finally, we mention that the defect pair in such configurations has eight nearest neighbor atoms (two O and six N atoms). Because the goal of work was to analyze the FM in GaN, we do not consider Ref. [15,20] 1.1 LDA Ref. [14] 0.65 1.2 LDA here 1NNs configurations with the seven nearest neighbor atoms. When the complex is spin polarized, the orientation of the two spins can be ferromagnetic (FM), ferrimagnetic (FiM, when spins of the two complexes are different because the defects are in different charge states), and antiferromagnetic (AFM). In the latter case, the total spin vanishes but the spin polarization energy is finite. Energy of magnetization ΔE AFM(FiM)-FM is defined as a difference between the total energies of AFM (FiM) and FM states and is positive when FM coupling is stable. In summary, two and four configura-  Table 4. In Table 4, we only briefly highlight the results. The obtained values are a little larger than those obtained by HSE06 approach in Ref. [26] due to shorter spatial separation between defects in our work.  Table 4).
The dominant mechanism of magnetic interaction between V Ga O N is determined by the interplay between the counteraction of the bonding-antibonding (BA) and the exchange of spin-up-spin-down (Δe ex ) splittings.  Table 4 Energy of magnetic coupling ΔE AFM(FiM)-FM together with total magnetic moment (in μ B ) of the complex in the charge state q (it is a sum q 1 and q 2, where) obtained for zb-and w-GaN calculations. For simplicity, both ΔE AFM−FM and ΔE FiM-FM are denoted by ΔE AFM−FM and the actual spin configurations and values of magnetic moment are given in the columns "μ" comprised of long-range tails which involve p orbitals of distant N ions (Fig. 5). Moreover, N neighbors of the complex are not equivalent, as displayed by the dissimilar contributions to the spin density. Indeed, the contribution of the axial N (or O) ions is almost vanishing (Fig. 5 a, c, e, g). Non-negligible spin density of In orbital demonstrates the important contribution of In atom into magnetic interaction between defects.

Crystal Structure zb-and w-GaN and Relaxation
All above results are demonstrated for atomic relaxed structures. The increase of FM stability when In content grows can be explained by the distortion effects in crystal structure due to the atomic displacement in the process of relaxation. As it can be seen in Figs. 2 and 5, significant and complex perturbations in the crystal structure of GaN are experienced when vacancy is introduced and it can be attributed to the fact that the radius of In is larger than the one for Ga and the radius for N is larger than the one for O. The interatomic distances between the neighbors of a complex defect are mainly determined by the lattice constants of the host crystal (but here, for low In content, it was neglected), and also by the atomic relaxations, i.e., displacements of the neighbors. In the studied compounds, the outward relaxation of the nearest neighbors elongates bonds by about 5-12%, i.e., 0.   Fig. 6, it follows that distortion (change in bond length) for the V Ga O N In Ga is larger than that for the V Ga O N . But we note that structural distortions are more complex. Although, the displacements of the second and third neighbors are an order of magnitude smaller, the effect of atomic relaxations around defects, involving not only the nearest but also more distant neighbors, cannot be neglected. The states of the defect complexes are determined by the overlap of the N and O dangling bonds given by the N-N (N-O) distances. In ideal structure (after relaxation without defect), N-N is equal to 3.18 Å. In V Ga O N , the N-N and N-O dangling bonds are 3.55 and 3.48 Å, respectively. But in V Ga O N In Ga , these values are shorter; three N-N are 3. 58, 3.52, and 3.48 Å, respectively, and N-O is 3.42 Å. It implies that the defect states are more localized, and the energy of spin polarization of such defects is higher.

Summary and Conclusions
In summary, spin states of V Ga O N and V Ga O N In Ga complexes in both zb-and w-GaN, and magnetic coupling between them, were studied within GGA+U calculations. The U(Ga) =3 eV and U(N) = 5 eV terms were imposed on d(Ga) and p(N) leading to the correct band gap of GaN.
Charge states q from 0 to − 2 of V Ga O N were considered. In both crystal structures for neutral V Ga O N with S = 1, high-spin configuration is stable. Wave functions of V Ga O N have a multi-orbital character, being composed of three p(N) and one p(O) orbitals of vacancy neighbors. But the main contribution to the spin density comes from sp 3 N orbitals.
Two different electronic structures, which arise due to two different geometry configurations of V Ga O N , with O N either along the c-axis and in one of three equivalent tetrahedral positions in w-GaN, were analyzed. The latter geometry configuration assumes stronger stability of HS state and more delocalization of defect state.
Introducing In Ga as second neighbor to V Ga on the one hand imposes changes to the electronic structure, on the other gives rise to the delocalized wave function of the defect as the crystal structure is perturbed, and finally contributes to long-tail spin density distribution. Magnetic moments originate mainly from sp 3 N orbitals, and the contribution of p orbitals of distant N, d(Ga), and d(In) states is about 20%.
Various relative orientations of the defects and several charge states (q = 0 and q = − 2) were considered, and consequences regarding the observed FM in GaN were pointed out. Using a relation predicted from mean field model, T c = 2zS(S + 1)J/3k B , the room temperature of FM implies that (assuming z = 6 neighbors and S = 1) the coupling constant J(V Ga O N In Ga ) = 0.01 and J(V Ga O N ) = 0.0045 eV, i.e., 920 K and 415 K, respectively. These values have a limited reliability as the distribution of defects is random and the coupling depends on the distance between defects.
Comparing the obtained results with experiments, we note that, according to the results, the observed collective ferromagnetism in GaN systems [7][8][9][10][11] can originate from magnetic interaction between V Ga O N defects. And in Ga-rich InGaN alloys, we predict even stronger FM.