Donor and acceptor characteristics of native point defects in GaN

The semiconducting behaviour and optoelectronic response of gallium nitride is governed by point defect processes, which, despite many years of research, remain poorly understood. The key difficulty in the description of the dominant charged defects is determining a consistent position of the corresponding defect levels, which is difficult to derive using standard supercell calculations. In a complementary approach, we take advantage of the embedded cluster methodology that provides direct access to a common zero of the electrostatic potential for all point defects in all charge states. Charged defects polarise a host dielectric material with long-range forces that strongly affect the outcome of defect simulations; to account for the polarisation, we couple embedding with the hybrid quantum mechanical/molecular mechanical approach and investigate the structure, formation and ionisation energies, and equilibrium concentrations of native point defects in wurtzite GaN at a chemically accurate hybrid-density-functional-theory level. N vacancies are the most thermodynamically favourable native defects in GaN, which contribute to the n-type character of as-grown GaN but are not the main source, a result that is consistent with experiment. Our calculations show no native point defects can form thermodynamically stable acceptor states. GaN can be easily doped n-type, but, in equilibrium conditions at moderate temperatures acceptor dopants will be compensated by N vacancies and no significant hole concentrations will be observed, indicating non-equilibrium processes must dominate in p-type GaN. We identify spectroscopic signatures of native defects in the infrared, visible and ultraviolet luminescence ranges and complementary spectroscopies. Crucially, we calculate the effective-mass-like-state levels associated with electrons and holes bound in diffuse orbitals. These levels may be accessible in competition with more strongly-localised states in luminescence processes and allow the attribution of the observed 3.46 and 3.27 eV UV peaks in a broad range of GaN samples to the presence of N vacancies.


Introduction
GaN is a key material for many technologically important applications including blue-light emitting and laser diodes [1], high-power microelectronics [2], solar cells [3,4] and catalysis [5]. Although the properties of native defects in GaN govern many fundamental characteristics of the material, such properties are not well understood, despite much research, and, indeed, which defects dominate is still a major topic of debate. GaN belongs to the family of III-V compound semiconductors, which has been widely studied experimentally and theoretically [6]. It stabilises in the wurtzite phase under ambient conditions. The most important native point defects to consider can be proposed by analogy with other III-V systems [7,8], which are much better characterised, such as GaAs, where vacancies and antisites are thought to dominate [9,10]. A comprehensive study of native point defects in GaN should necessarily take into consideration lattice vacancies, self-interstitials, and antisites.
A conclusion drawn from many studies is that the n-type character of as-grown GaN is due to nitrogen vacancies [11,12], while others claim it is due to unwanted impurities [13,14]. The various photoluminescence (PL) peaks observed in undoped GaN have also been attributed to many different native defects [15]. Point defects in GaN can exist as either donors, which under certain conditions donate electrons to the conduction band, or acceptors, which donate holes to the valence band. In applications, the criterion for the donor or acceptor nature is that the activation energy should be comparable with the thermal energy. Of course, many defects can both donate and accept electrons, which results in the change of their charge states, and may have donor and acceptor levels deep inside the band gap. Native defects can act as compensation, passivation, and recombination centres during doping, which makes understanding their properties a prerequisite to study impurities and extended defects. Due to the wide variety of defects, charge states, and ionisation levels possible, deciphering experimental results with respect to particular defect species is a difficult process requiring significant input from theory.
Two pioneering density functional theory (DFT) studies using the plane-wave supercell technique provided a comprehensive description of the structural and energetic properties of native point defects in wurtzite GaN, both performed over twenty years ago [16,17]. The compensation effect of native defects was stressed in the studies, and vacancies are identified as dominant in GaN, with N vacancies (V N ) dominating in p -type and Ga vacancies (V Ga ) in n-type GaN. Boguslawski et al found that the V N introduces a shallow donor level and may be responsible for the n-type character of as-grown GaN [16], in agreement with earlier studies using tight-binding methods [18,19], while Neugebauer and Van de Walle argued that the defect has high formation energy under n-type conditions and thus could be excluded as the source of the n-type conductivity in as-grown GaN [17]. Both groups reported that the split-interstitial configuration was most favourable for the N interstitial (N i ). Neugebauer and Van de Walle further argued that the large mismatch in the covalent radii of Ga and N renders the antisites (Ga N ) and (N Ga ) and Ga interstitials (Ga i ) energetically less favourable than vacancies [17]. Boguslawski et al, in contrast, found the concentration of Ga i ([Ga i ]) to be comparable to that of the vacancies (implying similar energies of formation) [16]. Finite-size effects resulting from employing the supercell method were not taken into consideration in either study due to restrictive computational loads and lack of a theoretical framework. Contemporary studies of native defects in the alternative, higher symmetry but metastable zinc blende (cubic) phase of GaN showed similar defect formation energies and electronic structures to those in wurtzite GaN [20][21][22].
The early results summarised above presented a baseline in our understanding of defect processes in GaN, which has been exploited in many experimental studies to interpret available or new data. With the development of new methods and increase in computational power, more sophisticated theoretical approaches have become more affordable in recent years and some of the early findings were challenged . While those earlier studies did report defect levels, their comparison with experiment remained unclear. Indeed, the resulting carrier concentrations derived from the defect formation energies were considered the most appropriate properties to compare with experiment. In later studies, both adiabatic and optical levels were calculated and compared with the widely-available PL and other spectroscopic data. The introduction of hybrid density functionals to the plane-wave DFT approach, where a fraction of exact Hartree-Fock exchange is included in the density functional, allowed more accurate electronic structures to be calculated, although the exact implementation of this method varied from study to study.
Despite many theoretical and experimental investigations on native defects in GaN, there remain many controversies and unresolved problems, which is partly due to the difficulty of directly detecting dilute point defects in experiments. Examples include the sources of yellow, green, red, blue and ultraviolet PL peaks observed in many samples; the source of intrinsic n-type conductivity in GaN; the stability of the N antisite; and the concentration of Ga interstitials. The commonly used DFT supercell approach for charged defects calculations in semiconductors is well known to suffer from finite-size effects, i.e. spurious electron confinement (in common with cluster techniques) within the supercell and periodic image interactions, which are particularly problematic for charged defects. Moreover, the periodic boundary conditions do not allow for a unique definition of the energy reference level for different unit cells and charge states [49][50][51][52]. Many efforts have been made to resolve these issues [42,[53][54][55]. The resulting variety in approaches is a plausible reason for the lack of consistency within different calcul ations and in comparison to experiments [42]. Another point is the choice of density functional, which has a significant effect on the calculated energy levels [56]; the level of exact exchange employed in hybrid DFT has varied over many different studies, which also can account for the range of different results.
Here, we report a systematic hybrid quantum mechanical/ molecular mechanical (QM/MM) embedded cluster study on the structural, electronic and optical properties of native defects in GaN. Our goal is to achieve a reliable, accurate and reproducible account of native point defects in wurtzite structured GaN using a method that avoids spurious interactions inherent in periodic models and thus provide a reference for comparison of computational and experimental studies of this material. By design, the embedded cluster method accounts for the short and long-range polarisation effects of charged defects on the host material; lacks periodic image interactions; and provides an unambiguous definition of the vacuum level, facilitating the calculation of ionisation energies with an absolute reference. Highly delocalised states, however, are poorly represented due to quantum confinement effects (within the QM region of the model), which, incidentally, can also enter the description of diffuse defect states by supercell methods. Efforts are underway to overcome this problem [57], but we note that, for sufficiently compact defect states, such as the majority of those studied here, confinement will present no significant problem, provided a suitable basis set is employed to describe the electronic levels. Moreover, the methodological framework provides access to a broader range of quantum chemical methods including hybrid density functionals than is commonly available in plane-wave DFT codes. We therefore employ two hybrid functionals, both of which have been designed to reproduce a comprehensive set of thermochemical data well, but one of which is also designed to reproduce kinetic data. We have indicated in a previous study that one of these functionals gives results in oxides that are very close to commonly applied hybrid functionals in planewave studies [68]; this functional then provides a benchmark for our calculations.
Our results show that vacancies will dominate in GaN, N antisites are highly unfavourable, N interstitials stabilise in split interstitial configurations, Ga interstitials preferentially occupy octahedral sites, and that thermodynamic concentrations of defects will be low apart from some special cases. We give a fully comprehensive list of all vertical transition energies, including both electron and hole ionisation, which may be used by experimentalists for comparison with various spectra. Our results broadly agree with earlier studies and can explain many experimental findings including PL, deep level transient spectroscopies and complementary techniques.
The remainder of the paper is structured as follows: in section 2 we provide details of our computational approach. In section 3 we present our results, beginning with the formation energies of each point defect studied, then the thermal transitions, defect concentrations, and the ionisation energies of all relevant charge states. Finally, we discuss implications of calculated defect energies for experimental observations and highlight remaining problems in defect assignment, and comment on the description of shallow states in hybrid QM/MM calculations.

Computational techniques
The hybrid QM/MM embedded cluster technique employed is designed to study localized states in ionic solids, where charged and strongly dipolar species predominantly interact via long-range electrostatic and short-range exchange forces [58,59]. This method and a plethora of closely related embedding approaches [60] split extended systems into the inner region, containing the central defect, described using molecular QM theories, and its surroundings, which are only slightly perturbed by the defect, modelled with MM approaches. Our choice of the QM methodology is DFT; the MM simulations employ the polarisable shell model interatomic potentials; and the interface between these two regions is based on cation centred semilocal pseudopotentials. The outer part of the MM region is held fixed at the pre-optimised geometry [61][62][63] and any polarisation here is treated a posteriori (see below). The entire cluster is surrounded by point charges, which are fitted in order to reproduce the Madelung potential of the infinite crystal within the inner QM region. A schematic of the model is shown in figure 1.
In our method, the inner cluster of 116 atoms of wurtzite GaN (pre-optimised using the MM model, see [61][62][63]) centred on the defect is treated using (i) the second-generation thermochemical hybrid exchange and correlation (XC) density functional B97-2 [64], which is similar to those commonly used in recent plane-wave supercell calculations (21% exact exchange compared with 25% for PBE0 [65]or HSE06 [66]), (ii) the SBKJC small-core pseudopotentials on Ga [67] within the cluster and large-core refitted pseudopotentials [61,62] in the interface that provide a short-range contribution to the embedding potential on the defect, and (iii) the atomic basis set of def2-TZVP form on N [68] and matching SBKJC basis on Ga [67]. For comparison, we use a second hybrid XC density functional employing 42% exact exchange (BB1k) [69], fitted to reproduce kinetic barriers and thermochemical data, which, given previous work on oxides [70][71][72][73], we expect gives a more accurate description of electron and hole localization than B97-2-see also our earlier work on GaN in [62,74]. This QM region of radius 6.8 Å is embedded in an outer cluster of radius 30 Å, which is treated with an MM level of theory using two-body interatomic potentials parameterised to reproduce the wurtzite GaN bulk structure and physical properties [62,63]. A possible mismatch between the description of the MM region and the QM region could introduce strain erroneously in the model; however, both the MM model and QM hybrid density functionals we employ reproduce the low temperature experimental lattice parameters and bond lengths very well, resulting in a small mismatch (see the supplementary material (stacks.iop.org/JPhysD/52/335104/mmedia) for a demonstration of the reproduction of experimental bond lengths within the centre of the QM region using either functional). The method has been implemented in the ChemShell package [58,59,75] that employs Gamess-UK [76] for the QM and GULP [77] for MM single point energy and gradient calculations. Further technical details are discussed elsewhere [58,78]. This method has been previously applied to gain insights into the defect properties of a range of ionic systems [50,73,[79][80][81] including earlier studies of GaN [57,62,82] and, moreover, the band alignment of polymorphs of TiO 2 [62,83,84].
Defect levels of various defect charge states we report are calculated using their ionization energy, which consists of vertical (optical) or adiabatic (thermal, or thermodynamic) transitions of electrons to vacuum or, alternatively, the conduction band of a semiconductor. Periodic calculations confirm that the VBM in GaN is composed of N 2p states, which, although extended throughout the infinite periodic system, remain well localised on anion sites. Using the embedded cluster approach with a suitably terminated QM cluster (i.e. with cations as terminating ions and an appropriate embedding potential), we can model such states well. Therefore, by calculating the ionization energy of an electron from the bulk system, we can probe the VBM and determine its position relative to vacuum.
Using the B97-2 functional (BB1k functional), we calculate the VBM to be 6.625 eV (7.340 eV) below the vacuum level, which is comparable with 6 eV as calculated using periodic slabs with the PBE-GGA [85] functional, 7.2 eV from surface GW calculations [86], and 6.8 eV from experiment [87]. The conduction band states, however, which are highly delocalised, are underbound in our embedded cluster calculations. To correct for this quantum confinement effect, we shift their position so that their minimum lies at an energy equal to the band gap above the VBM, relative to vacuum. For this procedure we use the experimental band gap value of 3.5 eV [88].
The formation energy of a defect X(E f [X]) is determined from the relevant defect reaction as: where ∆E(X) is the difference in energy between the embedded cluster with and without X, n i is the number of atoms of species i added (n i > 0) or subtracted (n i < 0) in forming X, µ i is the chemical potential of species i, q is the charge of X, and E f is the Fermi energy. µ i depends on the experimental growth condition, which can be either N or Ga rich [89]. The atomic chemical potentials we use in our calculations are discussed in the supplemental material.
For adiabatic processes involving charged defects, to account for the polarisation in the outer MM region that is kept fixed in our simulations (see figure 1), we use Jost's form ula [58]: where E pol is the polarisation energy, Q is the charge of the defect in the QM region, R is the radius of the active region (i.e. where relaxations are treated explicitly) and ε 0 is the static dielectric constant. For vertical processes, however, where the charge state changes from Q to Q + ΔQ but only electronic relaxation can occur, ionic polarisation due to the charge Q will remain 'frozen in' at the final state, Q + ΔQ. The generalised Jost's formula then takes the form: is the vertical polarisation energy, and ε ∞ is the high-frequency dielectric constant. This expression is an update of that which we have used in previous works [58,62,90], hence there are minor differences in energy for some optical processes that we report here in comparison with those presented previously. For the values of the dielectric constants (ε 0 = 10.09 and ε ∞ = 5.37, see the supplementary material) we use one third of the trace of the corresponding dielectric tensors, calculated using the interatomic potentials [61,63].
Polarisation energies in dielectric materials by charged centres in the QM/MM methodology include contributions from the QM, active MM and frozen MM regions, and Jost's formulae can provide a useful estimate of the two MM contrib utions using ratios of their respective radii. We thus calculate up to ca. 30% of the MM defect relaxation energy The inner cluster (red spheres) is treated at a QM level of theory, and the surrounding atoms are modelled using an MM level of theory, divided into an active region (blue spheres) where ionic and electronic interactions are accounted for explicitly, and a frozen region (green spheres) where polarisation is treated a posteriori. Embedding pseudopotentials are placed on cation sites within the interface (yellow spheres) between the QM and MM regions, and the entire cluster is terminated by point charges (orange spheres) chosen to reproduce the Madelung potential of the infinite crystal.
(more than 3 eV, for example, for triply charged defect states) to be due to the long-range polarisation effects. A similar conclusion can be drawn for periodic models, where long-range polarisation is usually not considered, but may prove to be vitally important in the total balance of defect energies.

Results
A summary of all native defects including their atomic and electronic structure along with vertical levels is given in the supplementary material while here we focus on the energetics of these defects and consequent observed physical properties.

Formation energies
We first consider the thermodynamic stability and adiabatic transitions of the point defects we consider in this work, viz. N and Ga vacancies, V N and V Ga , interstitials, N i and Ga i , and N antisite, N Ga . (Based on the Madelung potential and normal oxidation states of Ga, we do not consider the Ga antisite as a plausible point defect in GaN across the whole range of relevant chemical potentials.) Our results are presented in figure 2. Note that the relative stability of a particular defect or charge state does not preclude formation by kinetic means, which could be observed experimentally.
We find, using either the B97-2 or BB1k XC density functional, that V N is the most favourable native point defect for Fermi energies within the band gap, apart from a small range close to the conduction band under N-rich conditions with B97-2. It is a donor-like defect throughout the band gap, stabilising in the +3 charge state above the valence band with a transition to the +1 state in the middle part of the gap. The +2 state of V N has high formation energy and is not favourable at any value of Fermi energy, indicating a negative U nature of the defect [91], which is analogous to isoelectronic F centres in many oxide materials as discussed in section 4 in the supplementary material. The formation energy of V N in the +3 state near the valence band is negative, implying a spontaneous formation that would compensate any p-type carriers if the Fermi energy were close to the VBM. This behaviour fundamentally precludes any macroscopic samples from being p-type under thermodynamic equilibrium as any acceptors will be compensated by the vacancies; we discuss this point further in the next section. V N undergoes the (+3|+) transition at 1.29 eV (B97-2) and 1.83 eV (BB1k) above the VBM. Ganchenkova and Nieminen [32] reported that the −1 and −3 states of V N stabilised in the band gap below the CBM, in contrast to our calculations that show the stabilisation of the −1 charge state 0.32 eV (0.88 eV) above the CBM, determined using the B97-2 (BB1k) functional. Open-shell charged states 0 and −2, like the state +2, prove to be thermodynamically unstable under all physical conditions. Under N-rich conditions near the conduction band, per the B97-2 calculations, V 3− Ga , the only triple acceptor among the native defects in GaN, becomes more thermodynamically stable than V N . This defect is a donor in the +1 charge state near the VBM, and it stabilises as an acceptor in the band gap in the charge states 0, −1, −2 and −3. With both functionals, near the VBM, V Ga is found to be highly unstable.
Like the V N , the Ga i is most favourable in the +3 charge state above the VBM and in the +1 state below the CBM, thus it remains a donor across the band gap. Ga i undergoes a (+3|+) transition in the upper half of the band gap, at 2.65 eV (B97-2) and 3.32 eV (BB1k) above the valence band. The formation energy of Ga 3+ i near the VBM is negative (except under N-rich conditions with the B97-2 functional, where the value is slightly positive, ca. 0.25 eV), implying a spontaneous formation, compensating positive free charge carriers, similar to the case of V 3+ N . The formation of N interstitial defects can be considered as a competition between available octahedral and split interstitial sites. The octahedral site proves to be stable in the -1 state for BB1k across the band gap, while with B97-2 it undergoes a (−1|−2) charge transition at 0.11 eV below the CBM. This configuration, however, is unstable relative to the split interstitial for Fermi energies up to 1.37 eV (with B97-2) and 1.86 eV (with BB1k) above the CBM. On transition from the octahedral to the split geometry, the structure autoionizes, with an energy gain in excess of 2 eV. The B97-2 calculations show a (−1|0) thermodynamic transition at 0.2 eV below (with B97-2) and 0.28 eV above (with BB1k) the CBM. On lowering the Fermi energy into the band gap, the cationic charge states become more favourable-see figure 2.
We employ the climbing-image nudged-elastic-band method to investigate the transition from an octahedral geometry to a more favourable split interstitial configuration in the −2 charge state; the result is shown in figure 3. Note that this transition is spin allowed (both configurations are of the same spin ½), whereas in the thermodynamically stable charge state of −1, the octahedral site is in a quasi-atomic state of 3 P 2 , while the ground state of the split interstitial is 1 Σ g , i.e the transition is spin forbidden. We determine an energy barrier of 0.03 eV, which is comparable to the thermal energy at room temperature. Therefore, we expect a transition from the octahedral to the split interstitial configuration to occur spontaneously under normal conditions.
Antisites include nitrogen-on-gallium (N Ga ) and galliumon-nitrogen (Ga N ). As discussed above, Ga N will not be stable because of: (i) the large difference in the relative effective sizes of gallium and nitrogen ions; and (ii) the electron-poor valence shell of Ga (i.e. the electropositive nature of the metal Ga), which is highly unlikely to support Ga anionic states even in the presence of a favourable Madelung potential. Nitrogen, in contrast, can act both as an anion and cation in different chemical environ ments. N Ga is a donor-like defect near the VBM as it stabilises in the +3 charge state, but as the Fermi energy increases trans itions occur first to the +2 state at 1.15 eV (B97-2) and 1.71 eV (BB1k) above the VBM, then to the neutral state at 2.09 eV (B97-2) and 2.83 eV (BB1k) above the VBM. Increasing the electron concentration further, using B97-2 we find that N Ga stabilises in the −1 state below the CBM, while it remains neutral when using the BB1k functional.
When comparing our results with previous studies for charged defects (see table 1), we find lower (higher) formation energies for positive (negative) charges and hence transition levels closer to the CBM and, equivalently, the vacuum reference energy. Therefore, in some cases, we disagree with previous work on predicted lowest energy charge states. We note, however, that the calculated formation energies of neutral native defects are very close to those reported by other groups, as shown in tables 2 and 3. Gillen and Robertson [42] report that, for V Ga , the stronger exchange interactions in the self-interaction-corrected sX-LDA approach result in defect levels higher in the band gap relative to the VBM, compared with standard LDA calculations. Lyons et al calculated thermodynamic transition levels that are even higher above the VBM than those determined from sx-LDA calculations, and they link this behaviour to hole localization [44]. They emphasised the importance of full geometry optimisation with a hybrid functional and, as a demonstrative example, showed that hole localisation on nearest-neighbour N ions to a V Ga is accompanied by significant relaxation of those N ions [44], which we also observe in our calculations (see table S4). The trend for higher in-gap levels is confirmed by our hybrid DFT calculations, and, for the V Ga , is even more pronounced than that of Lyons et al who used the HSE functional. Indeed, we observe such a trend for all defects studied, with a stronger effect resulting from using the BB1k functional rather than the B97-2 functional. This observation is consistent with the more localised spin densities that result from BB1k functionals, as discussed in supplemental material.

Defect and carrier concentrations
From the computed formation energies presented in the previous section, we can calculate the equilibrium concentrations and self-consistent Fermi energy, E F at a given temperature T, under the constraints of thermodynamic equilibrium and charge neutrality, using our program SC-FERMI [94,95]. The program works by calculating the concentration of each defect where N X is the density of sites in which the defect may form, g q is the degeneracy of the state q (derived from the analysis in the proceeding section) and k is Boltzmann's constant, as well as determining the electron (n 0 ) and hole (p 0 ) concentrations  Under equilibrium conditions (see figure 4), we find that V N is the most dominant defect, but that the concentrations are low, as one would expect from the formation energies (figure 2). Using either XC functional results in very similar calculated concentrations and E F . For N-poor conditions, [V N ] > 10 13 cm −3 at T ca. 800 K, while n 0 = [V N ] for the temperature range studied (0 < T < 1500 K). E F remains in the band gap, from ca. 1.3 to 0.5 eV below the CBM, moving closer to the CBM as the temperature is increased. In N-rich conditions, due to the higher formation energies, n 0 and [V N ] Defect formation energies in the neutral state in N-poor conditions. Label ~ indicates values estimated from figures in references. Method used in: ab initio MD [16], LDA [17,23,28,92], PBE [25], HSE (28%) [38], LSDA [33], HSE (31%) [41], sX-LDA [42], GGA [93]    are several orders of magnitude below those of the N-poor case, and at higher T (above 1200 K), minority carrier concentrations (p 0 ) of one (for B97-2) and two (for BB1k) orders of magnitude below that of the majority carriers are obtained. E F remains closer to the middle of the gap, varying from ca. 1.7 to 1.2 eV below the CBM as T is increased up to 1500 K. These results demonstrate the intrinsic n-type nature of GaN, although the calculated carrier concentrations are low (below 10 17 cm −3 ). Our procedure also allows us to analyse compensation of ionised donors or acceptors, by introducing a fixed concentration of either and determining the resulting defect concentrations and Fermi energy given the charge neutrality constraint. Turning first to donors: in figure 5 we show the calculated E F and defect and carrier concentrations when a fixed ionised donor concentration [D + ] = 10 19 cm −3 is present (D here could represent for example a Si Ga impurity, assumed to be fully ionised). Because of the lack of low formation energy negatively charged defects (see figure 2), even close to the CBM, [D + ] is uncompensated and n 0 = 10 19 cm −3 over the full range of T studied. In N-poor conditions, we do see a small concentration of V N at high T, which is several orders of magnitude lower than [D + ]. E F decreases as T increases, following the curve expected for an ideal gas of fermions. In N-rich conditions, with the B97-2 functional, the formation energy of V Ga is below that of V N close to the CBM, but the resulting concentration of the order of 10 12  For acceptors, the situation is quite different. Shown in figure 6 are E F and the defect and carrier concentrations when a fixed ionised acceptor concentration of [A − ] = 10 19 cm −3 is present (here A could represent, for example, an ionised Mg Ga ). In this case, the acceptor is compensated by positive V N , rather than holes, while E F remains within the band gap, at least 0.9 eV above the VBM, and moves further from the VBM as T increases. Under N-rich conditions, using either XC functional results in E F that varies between 0.9 and 1.7 eV, with compensation by At T > 700 K, thermal excitation of holes results in p 0 > 10 13 cm −3 , which is still several orders of magnitude below that of [A − ]. Using the BB1k functional, we find that concentrations of other positively charged defects (Ga i and N i ) become greater than 10 13 cm −3 at T > 1300 K. Under N-poor conditions, the results are similar; but in this case, due to the lower V N formation energy, E F remains closer to the middle of the gap. As a consequence of this position of E F , at T > 400 K (600 K) for B97-2 (BB1k) the lowest energy state of V N changes from 3 + to 1 + , resulting in an increased compensation concentration (i.e. the compen- . At elevated temper atures, thermal excitation of carriers occurs, but as E F is closer to midgap both carrier types are observed. From these results, we can deduce how point defect compensation further restricts effective acceptor activation. We focus on the N-rich conditions, where V N are less favourable to form and, therefore, one would expect acceptors to be compensated less. We have assumed that the acceptors are fully ionised, but the best-case scenario for p -type activity is at T = 1500 K for B97-2, where p 0 = 10 −2 [A − ]. For Mg impurities, many studies claim that the ionisation energy is ~0.2 eV [15,96]. Our results show that, even at 1500 K, where one would expect p 0 of the order of 0.1 [Mg − Ga ] from the acceptor ionisation energy, the actual p 0 would be at most 10 −3 [Mg − Ga ]. At lower T, p 0 decreases to many orders of magnitude lower than that of the acceptor concentration. Therefore, in order to see any significant p -type activity, the compensating V N must somehow be controlled, most likely through forming complexes and/or more extended defects, which would explain why no p -type activity is observed unless either the Mg concentration is very high (at least 10 19 -10 20 cm −3 ), or particular structuring and design is employed (e.g. thin films [97], capping layers [78]). Under N-rich conditions, that help to suppress N vacancy formation, keeping the V N concentration one order of magnitude below the acceptor concentration will result in a 70% activation of the acceptors at room temperature (in our example, p 0 = 7 × 10 18 cm −3 ). Furthermore, we find that, at elevated temperatures, increasing [A − ] results in further compensation by other positive defects, particularly [N 2+ i ], but also [Ga 3+ i ] when using the BB1k functional. One of the puzzles in the defect chemistry of GaN is the apparent contradiction between the high formation energies for Ga vacancies calculated consistently by different groups [23,25,27,28,32,42,43,92] and several reports of substantial concentrations, up to 10 20 cm −3 , of this defect as seen in positron annihilation spectroscopic experiments [98][99][100][101].
Our calculations, as can be seen from figure 4, do not confirm the presence of such concentrations under any conditions in a pristine (formally undoped) material. A resolution for this contradiction has been presented by Saarinen et al who have shown that the Ga vacancies in GaN compensate for extrinsic donors such as O impurities and their concentration drops below 10 15 cm −3 when [O] is reduced [100]. Moreover, V Ga -O N defect complexes have been proposed as the main form in which the vacancy is accommodated by the material [102]. Modelling defect complexes is outside the scope of the cur rent study, but we do expect a significant stabilization of both the vacancy and impurity by the Coulomb interaction between them as nearest or next nearest neighbours, which would result in a substantial rise in the observed vacancy concentration; see [82].

Ionization energies as defect (transition) levels
The knowledge of defect formation energies in different oxidation (or charge) states-see [103]-alongside the vertical ionis ation potentials and electron affinities provides insight into many defect processes that play a key rôle in photoabsorption and luminescence, charge trapping, etc. In this analysis, we use a conventional configurational coordinate diagram (see [104]) shown in figure 7. For example, an electron trapped by an acceptor in the ground state (A − ) can be excited to the conduction band if it absorbs a photon with sufficient energy (E ab ), the resulting configuration being the neutral charge state plus a conduction electron, A 0 + e − . The defect in the ionised state will relax ( E * rel ) before emitting a photon and returning to its ground state via an electron-hole recombination. The relaxation lowers the energy of the ionised state as the surrounding atoms move to new stable positions. On emission, a conduction electron undergoes a vertical transition from the conduction band to the empty defect state (A − ), i.e. recombining with the hole bound to the acceptor, giving rise to a PL (E PL ). Again, the atoms around the defect relax back to the initial ground state (E rel ) after the optical emission. An analogous process will involve a hole ionisation from an acceptor, which, in other words, describes a capture by the acceptor of an electron from the valence band-see the upper part of the configuration-coordinate diagram shown in figure 7. The energy difference between the minima of the ionised and ground states can be identified in the PL spectra as a zero-phonon line (ZPL). In the adiabatic approximation following the Franck-Condon principle, the vertical electron transition is assumed to be much faster than the atomic relaxation, which is usually the case in defect processes [15], although we note that other more complex processes involving other photo-generated carriers may also be relevant [74]. Furthermore, the total energies of both the ground and excited electronic states are shifted up by the corresponding phonon contributions, which results in the next order corrections to the transition energies. The defect-to-band transitions are also often associated with activation of charge carriers trapped at defect sites, so the ZPL energies are considered as activation energies in the relevant electrical processes, which is the topic of the next section. The relationship with the defect level diagrams discussed in section 4 in the supplementary material is as follows: the position of the donor (occupied) level of a defect in charge state Q with respect to the conduction band corresponds to the photo-absorption in a D(Q|Q + 1)e process, while its separation from the valence band gives the PL energy in a D(Q|Q − 1)h process; the corresponding acceptor levels yield the energies of photo-absorption in a A(Q|Q − 1)h process and that of PL in a A(Q|Q + 1)e process. Further calculations of the potential energy surface along the configuration coordinate could be used to derive line shapes and intensities, but such calculations are beyond the scope of the current work.
We now can compare our results, presented in table 4, with pertinent available experimental data from spectroscopic measurements. We focus on the values computed using the BB1k functional, which is expected to provide a better account of electron localisation and, consequently, more accurate ionisation energies. The B97-2 derived results are provided in table 4 for comparison.
The observed PL spectra from unintentionally and intentionally doped n-type GaN almost always contain a broad yellow luminescence (YL) band with a maximum at about 2.2 eV [15]. Originally, the yellow band has been assigned to a radiative transition from a shallow donor or from the conduction band to a deep acceptor [105,106]. Previously, V Ga as a native defect has been suggested to be involved in YL both on computational [42,43] and experimental [98,101] evidence. In our calculations (see table 4), the negative and neutral states of V Ga are non-radiative electron emission centres, and its +1 state can give rise to an emission peak for radiative electron capture in the infrared region, V Ga (0| + 1)e, which we discuss in more detail below. More positively charged states of the defect prove to be unstable and we conclude that the original assignment is not supported by our data. It has been, however, reported that the observed increase in intensity of the 2.2 eV band results from redistribution of holes released on excitation to the valence band [107]. Indeed, we find an emission peak for radiative hole capture by the −1 and −2 states of V Ga at 2.17 and 2.21 eV; this process, however, will only be one of the sources of the widely observed YL. The respective calculated ZPL energies of 2.75 and 3.13 eV for both processes are for example significantly higher than the recently reported value of 2.57 eV for the YL1 centre [108]. We note that we calculate consistently high (>0.5 eV) defect relaxation energies, which is not consistent with such a low ZPL energy (by implication E rel ~ 0.37 eV). We find only one other defect state with the PL in the YL range: an octahedrally coordinated interstitial N (N O i ) in the −1 state would produce an emission peak at 2.17 eV on radiative hole capture. The octa hedral interstitial configuration with a triplet electronic ground state is, however, less stable than the singlet split interstitial structure. If formed it would be a relatively long-lived species, for example, in materials under ionising irradiation with energy above that of the N displacement threshold or grown under nitrogen rich conditions, and GaN thin films or powders exposed to nitrogen under pressure. We speculate that other sources of YL relate to defect complexes and/or impurities, which we have not considered in the present study.
Yan et al [41] have claimed that, in p-type GaN, V N in the +3 state can make a transition to the +2 state and give rise to an emission peak at 2.18 eV, which is in the yellow region of the spectrum. We do not observe such a peak for V N .
Another sharp peak at about 3.25-3.27 eV is observed in undoped GaN at low temperatures, which is known as the shallow DAP (donor to acceptor pair) band and is caused by transitions from the shallow donors to the shallow acceptors   figure 7, and all values are in eV. Results are given for two density functionals, but in general, we expect more accurate values when using the BB1k functional. Optical processes involving relaxed resonant states are not included here.
Defect states B97-2 BB1k    [15,[109][110][111]. Employing the present analysis, we do not recognize this peak in our calculations-see, however, further discussion in section 3.5. In GaN layers grown by metal-organic chemical vapour deposition (MOCVD) or hydride vapour phase epitaxy (HVPE), a blue luminescence (BL) band is often observed in the lowtemperature PL spectrum at about 2.9 eV with a 3.098 eV ZPL [107]. It has been attributed to transitions from the conduction band or a shallow donor to a relatively deep acceptor having an ionization energy of about 0.34-0.4 eV [15]. We calculate an emission peak for radiative electron capture by the −3 and −2 state of octahedral N i at 2.88 and 2.76 eV, respectively, but this configuration is significantly less stable than the split N i , whereas the ZPL energies are much higher than the observed value. Another peak at 2.85 eV for radiative electron capture by the −2 state of N Ga is found; however, this charge state remains metastable across the band gap.
In undoped GaN grown by HVPE or MOCVD methods, a red luminescence (RL) band is sometimes observed at about 1.8 eV, and is attributed to a deep acceptor level at 1.130 eV above the VBM [111,112]. We determine an emission peak for radiative hole capture by the neutral state of N Ga at 1.80 eV, while the +1 state of N Ga is calculated at 1.15 eV above the VBM (E abs for N Ga (+1| + 2)e is 2.35 eV, see figure S6). Of course, given the high computed formation energy of the N Ga , these peaks would only be observed in samples prepared by non-equilibrium processes. Katsikini et al find a RL1 appearing 1.7 eV below the absorption edge that is attributed to N interstitials-which we cannot confirm-and a RL2 appearing at about 1.0 eV above the absorption edge which is attributed to nitrogen dangling bonds [113]. We find, however, that the +2 and +1 state of Ga i would give rise to emission peaks for hole capture at 1.89 eV and 1.83 eV. As Ga + i would have the highest concentration of these defects in n-type Ga-rich conditions, we expect that this defect is a significant source of red PL in the corresponding samples of GaN. In N-rich conditions however, the concentration of Ga i approaches that of N Ga and, perhaps, both defects would compete for recombination with the photogenerated holes, although of course the two defects could form a complex with different spectroscopic signatures.
Next, a broad green band (GL) at about 2.5 eV and a broad red band at 1.9 eV have been detected in highly resistive GaN samples grown by MBE under extremely Ga-rich conditions, and have provisionally been attributed to internal transitions in some defects [15,114]. We find an emission peak for radiative hole capture N 0 Ga at 1.80 eV, and another emission peak for radiative hole capture by N − Ga at 2.54 eV, both with large relaxation energies that would account for the broad nature of the PL and which could provide a reasonable alternative attribution-see figure 8. Note that the N 0 Ga level lies in the middle of the band gap.
Furthermore, the GL bands will be originating from other intrinsic defects as, for example, a split N 0 i has a green emission peak for radiative hole capture at 2.49 eV, and a V −3 Ga has a green emission peak for radiative hole capture at 2.35 eV, which should then be expected to be seen in GaN prepared in N-rich conditions. As seen from table 4, we find that a V 0 N and V − N , on a hole capture from the VBM, are likely to give rise to the PL peaks at 3.46 eV and 3.41 eV that are practically ubiquitously observed in GaN. We have identified the V N as the major compensating centre for acceptor impurities in GaN (see section 3.2 above), where hole capture by electron-rich shallow donors explains the near band edge emission (see [62,74]). We discuss the 3.46 eV emission peak in further detail in section 3.5.
Finally, a near infrared (IR) PL signal at 0.95 eV (see [115][116][117][118]) appears in electron irradiated samples; this PL band mainly vanishes on annealing. Our calculations allow us to ascribe this band to electron capture by a V + Ga (see figure 9). The V Ga defects have previously been proposed as the obvious result of the irradiation and a convenient position of the 1.1 eV hole capture transition by V −3 Ga from early DFT calculations. The PL has been shown to be correlated with magnetically detected L5 and L6 defect centres in GaN, which in turn were associated with Ga i , situated at differing positions with respect to V Ga . This band is anticorrelated to the RL observed before irradiation. Above, we have attributed the RL to either Ga i or N Ga . Electron irradiation should remove N Ga and result in a significant concentration of V Ga , which will bind electrons on nitrogen ions decorating the vacant site (reminiscent of the RL2 centres as proposed by Katsikini et al). On photoexcitation, significant hole concentrations will form at the vacancies (as the photoexciting frequency is sub-band gap), rather than in the delocalised band states, which will suppress the hole capture process involved in the Ga i RL, and promote the V Ga electron capture process described above. Therefore, we would expect the RL to reduce in intensity and the IR PL to increase on irradation, as is seen experimentally. Subsequent annealing will remove some of the vacancies, while those  remaining will adopt the lowest energy charge state (−2 or −3 under n-type conditions, see figure 2), which cannot luminesce in the infrared, resulting in a suppression of IR PL, again as seen experimentally.
In this section, we have compared our results to a selection of experimental PL studies. We do not attempt to be exhaustive in our comparison with the available literature, or to discuss the spectroscopic processes in intricate detail. Such a study on the luminescence properties of undoped GaN is beyond the scope of the present work. Our results given in table 4 should, however, provide a useful comparison tool for the interested reader.

Defect as traps
Defect transition levels and configurational diagrams could be usefully employed in considering diverse phenomena that are not limited to the optical processes. A large number of deep electron and hole traps (see e.g. [119,120]) have been reported using a host of modern electrical, dielectric and related techniques. In our analysis, we distinguish fast optical (vertical) and slow thermodynamic (adiabatic) processes. Interpretation of experiment, however, can become quite involved and, in most cases, will require a separate study. Depending on the relative thermodynamic stability of the terms corresponding to the initial and final states, we can recognise thermally or optically activated ionisations, followed by recombination accompanied by luminescence or radiationless processes. Below, we consider examples of attribution of experimentally reported trap levels based on our defect calcul ations employing the BB1K functional.
We calculate the energy of V 0 N at E c − 0.04 eV (see figure S1) in excellent agreement with a recent optical admittance spectroscopy studies of a 0.051 eV state [121] and its singly negative state at E c − 0.09 eV that falls in the range of ET1 [120]. These levels indicate the shallow donor nature of the vacancy, and were discussed in a previous publication [62]. Both states, however, are metastable (as seen from the negative ZPL energies), and the majority of defects are expected to be in the singly or triply positively charged states depending on the position of the E F -see our discussion of the corresponding defect concentrations in section 3.2 above. In this case the vertical ionisation potential is our estimate of the activation energy. A hole trapping by V + N at 1.69 eV is also close to the 1.76 eV energy of the H5 trap [119], for which we, however, find more candidates, e.g. a V N (+3| + 1)2h transition that lies at E V + 1.83 eV.
A dominant electron trap in variously prepared GaN is the E3 [119] (or ET10 [120]) centre, with its level at E c − 0.56 eV [122,123] and high hole capture cross sections; thus acting as an efficient recombination centre. The concentration of the centres with the optical ionization energy of ca. 1.3 eV was found to correlate directly with the intensity of the YL band. These measurements are consistent with our calculated optical ionisations levels of 1.15 and 1.29 eV for V −3 Ga and V −2 Ga , respectively, with the thermal transition between the two lying about 0.02 eV below E c . Under the thermal activation, the latter defect is, therefore, expected to dominate, and the hole recombination on the defect centre, we predict, will result in an optical transition in the YL range (emission at 2.35 eV). As the PL energy in the electron ionisation case is calculated to be negative, the activation energy is ill defined, which will qualitatively change the kinetic behaviour of the system. The E3 level could be attributed to V −2 Ga or, indeed, any of the three negative states of V Ga if the activation energy lies somewhere between vertical and adiabatic levels. Provisionally, the activation energy for V Ga (−1|0)e would thus be expected to be between 0.75 and 0.90 eV, which would place this transition in the energy range of ET12 donor centres whose presence, however, is usually linked to dislocations [120].
The presence of a split interstitial nitrogen has been convincingly demonstrated by single and double electron spin resonance techniques [31], which is magnetically active in the neutral charge state, N S0 i . We calculate its thermal ionisation level at E c − 0.12 eV, close to the activation energy of ET2 centre [120]. Its further thermal ionisation energy of 1.76 eV is not accessible by usual DLTS methods but we calculate a counterpart hole trapping in the N S i (+2| + 1)h process as 1.75 eV, which is again in a marked agreement (or coincidence) with the E V + 1.76 eV position of the H5 trap [119].
For the interstitial Ga, the experimental evidence is again from optically detected magnetic resonance [35], which is observed within a broad IR PL band around 0.95 eV. The IR PL band emerges on electron irradiation of GaN, and its magnetic structure has been attributed to the Ga +2 i . The negative values of E PL we calculate for the neutral and positively charged states of Ga i are indicative of the radiationless character of the involved processes and the corresponding bounds on the activation energies are rather wide. The neutral state is predicted to be autoionized, and the thermal ionisation energy of Ga +1 i to be 0.38 eV. Perhaps, more usefully, the energy of the Ga i (+1| + 3)2e transition is found to be 0.18 eV, which places this defect midway between ET3 and ET4 [120].
In n-type GaN, the N Ga defect adopts a neutral charge state with the thermal ionisation potential of 0.93 eV, which makes it a primary candidate for the ET12 centre or, possibly, a slightly lower energy ET11 [120,124], and which could also be related to the E8 traps as seen in [119]. A similar attribution has been proposed to a E c − 1.04 eV trap with the thermal activation of 0.90 eV from I-T characterisations by Polenta et al [125]. The activation energy of further ionisation N Ga (+1| + 2)e can be estimated to lie between 0.43 and 0.56 eV. It is unlikely that this process is responsible for the ubiquitous E c − 0.56 eV trap discussed above-we expect the antisite to be only a small minority defect, but it is consistent with the behaviour of rare ET8 traps at E c − 0.49 eV [126], which have previously been associated with N Ga, and show a barrier of about 0.15 eV for electron capture (see discussion in [120]). Finally, the hole capture by N +3 Ga at 1.71 eV is again consistent with the position of H5 traps [119], for which, however, we have already seen more candidates.
In conclusion we reiterate that identification of particular defect processes is not possible on the closeness of the calculated energies with the reported trap levels or activation energies, but will require further careful analysis.

Effects of delocalization
Analysis of the defect formation and ionisation energies presented above assumes a compact nature of the defect electronic states. Such an assumption, however, is only an approximation as in real dielectric materials any charged centre will necessarily trap one or more charge carriers in shallow hydrogenic states (see [57,127]), which are typically probed by DLTS and related techniques. Due to the highly delocalised character of such states the exact nature of the defect becomes less important and different defects of the same charge would have very similar energies. When comparing our results to PL and other experiments, we have followed the standard procedure of considering ionisation to and from bands. A more complete picture would involve alternative levels associated with the diffuse states, which, however, lies outside of the scope of this study. As a guide to the interested reader we provide the characteristic energies and effective Bohr radii of hydrogenic defect traps both for electrons (m * = 0.22 m e [from [57]]) and holes (m * = 1.50 m e [128]) in table 5. As one can see, the electron states in GaN are actually diffuse, whereas holes remain quite compact.
We find that the diffuse hole level is 0.2 eV above the VBM, indicating that a shallow acceptor would bind a hole with this ionis ation energy. The only means of achieving p-type GaN is through doping with high concentrations of Mg, and an associated level at ~0.2 eV is observed experimentally. An adiabatic ionisation energy of ~0.2 eV is commonly attributed to the (0/− 1) transition of Mg Ga from previous DFT studies. The diffuse hole level, however, indicates that a much shallower acceptor that binds effective mass like holes must predominate in p-type GaN, given the agreement with the experimentally observed level. Moreover, as discussed above, we attribute the ubiquitous 3.46 eV UV PL peak seen in undoped and p-doped GaN to the recombination of an electron bound to a neutral V N with a hole at the VBM (see figure 10). Another peak observed experimentally in a broad range of samples, in particular p-type GaN, is DAP peak at 3.25-3.27 eV. If, instead of recombining with a hole at the VBM, the electron bound to the V N recombines with a diffuse hole bound to a shallow acceptor, the emission energy will shift by 0.2 eV, resulting in a PL peak at 3.26 eV. The balance between these two competing processes can explain the observed peaks, and demonstrates the power of including the diffuse states in analysis of spectroscopic processes.

Summary and conclusions
We have reported a detailed systematic study of point defect processes in GaN from embedded cluster calculations. Quantitative energy levels of each native defect in GaN were provided. We find a shallow donor level of V N which accounts for the n-type conductivity of as-grown GaN. V Ga has a deep acceptor level. Ga i is a donor with its donor level in the conduction band. Split N i has a deep acceptor level. Neutral N Ga has an energy level near the middle and thus is stable in the middle of the band gap. Losing or gaining electrons introduces shallow acceptor or donor levels in the band gap. This character of N Ga makes it a plausible candidate for the exper imentally observed broad red and green luminescence that detected only in high-resistivity GaN samples. We also give formation energies of each defect, among which vacancies are the most favourable ones in different conditions. Our calculated ioniz ation energies for the different charge states of the various defects rationalize the observed PL.
Acknowledgments Z X thanks the China Scholarship Council (CSC) for support. We are grateful to Alex Ganose and Stephen Shevlin for technical help, and Chris Van de Walle and Audrius Alkauskas for useful discussions. The EPSRC is acknowledged for funding (EP/K038419; EP/I03014X; EP/K016288). A W was supported by the Royal Society. Computational resources were provided through the Materials Chemistry Consortium on EPSRC Grant No. EP/L000202.