Formation of intrinsic point defects in AlN: a study of donor and acceptor characteristics using hybrid QM/MM techniques

Abstract

The wide-gap material aluminium nitride (AlN) is gaining increasing attention for its applications in optoelectronics, energy, and quantum computing, making the investigation of its defect properties crucial for effective use in these fields. This study employs a hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster method and uses three different hybrid-level density functional theory (DFT) functionals (B97-2, PBE0, and BB1K) to investigate systematically the thermodynamic stability, electronic properties, and donor/acceptor nature of intrinsic charged point defects in AlN. Our approach allows for the examination of defects within a significantly larger simulation cell, enhancing the reliability of the calculations. Our findings highlight the stable structures of defects including the symmetry-breaking reconstruction of octahedral-centred nitrogen interstitial and nitrogen split-interstitial defects, as well as the potential of nitrogen vacancies and aluminium interstitial defects as sources of shallow donor species. Additionally, we compute equilibrium defect concentrations at different annealing temperatures up to 2800 K, elucidating the roles of nitrogen and aluminium vacancies and interstitials in the conductivity of undoped, n-type doped, and p-type doped AlN. This comprehensive study advances the understanding of defect stability and electronic properties in AlN, paving the way for its enhanced application in advanced technologies.

---

1 Introduction

Aluminium nitride (AlN) is an important material for energy,¹ optoelectronic,^2–4 and potentially quantum computing devices.^5,6 The ultra-wide band gap of AlN enables efficient emission of ultraviolet light in UV light-emitting diodes (LEDs)² for sterilisation and water purification, as well as excellent thermal and electrical insulating properties for high-power electronic applications.⁷ The investigation of defect properties of the material is crucial for understanding many of its fundamental characteristics and enhancing its performance in these applications.⁸ To control the electrical functionalities of the material, doping with various elements, including silicon for n-type⁹ and magnesium for p-type conductivity,^10,11 has been applied, although the process requires careful control to achieve the desired defect concentrations and minimise compensating defects. Characteristic defect signatures of specific luminescence emissions have been reported,⁸ aiding in the identification and understanding of these defects.

The many properties of AlN that relate to the defect structure remain unclear. It is widely considered that AlN tends to be n-type due to native defects such as nitrogen vacancies and oxygen impurities. In experiment, while the presence of N and Al vacancy defects has been identified by various techniques;^12–19 other intrinsic point defects including interstitial and antisite defects have not been directly observed. These native vacancies,^12,20–28 as well as related defect complexes with oxygen impurities,^{19,21,22,25,29–44} are proposed as possible sources of deep-state luminescence/absorption within the band gap. Other unintentionally contaminating or intentionally doped impurities include silicon,^23,45,46 carbon,^46–51 magnesium,^10,52–54 and other species^8,55 which are also found to be stable in the material. All these defects can act as either donors or acceptors, depending on their charge states, leading to different energy states within the band gap region. However, the imprecise control over eliminating unwanted impurities and introducing intended impurities leads to discrepancies and gaps in our knowledge, which persist regarding the nature and properties of these defects; and it is often difficult to identify whether a particular type of defect or a group of different defects is associated with a certain wavelength of photo-response. Consequently, a comprehensive understanding of these defects and their complex donor/acceptor behaviour is essential for the development and optimisation of AlN-based applications.

The stability and donor/acceptor properties of these intrinsic defects can be explored by computational techniques, such as the density functional theory (DFT) methods. Earlier DFT studies, using the local density approximation (LDA) and generalized gradient approximation (GGA) functionals, on defects in zinc blende and wurtzite phase AlN provided insights into structural and electronic configurations of intrinsic point defects, as well as many impurities.^56–63 They were able to explore several defect types that had not been identified in experiment, including nitrogen split-interstitials and antisites. However, these LDA and GGA simulations suffer from the self-interaction error, which causes excessively delocalised charge densities and an underestimated band gap value, resulting in inaccurate representation of charged defect species for the Fermi level within the gap region. In recent years, the rapid growth in computer power has enabled the use of more chemically accurate methods with larger simulation cells. The Heyd–Scuseria–Ernzerhof (HSE) hybrid GGA functional has been widely used in defect studies of AlN,^64–76 as the fraction of Hartree–Fock exchange is tuneable for better agreement with experimental physical properties including lattice constants, enthalpy of formation, and band gap.⁷⁵ Experimental photo-absorption and photoluminescence energy bands have often been associated with the optical transitions involving the deep and shallow states of defects, as calculated in those studies using hybrid DFT functional methods.

Current DFT methods are mostly implemented for a system constructed with periodic boundary conditions (PBC) using plane-wave based software. For charged defects, such an implementation can give an ambiguous description of the reference potential level, and therefore an unreliable description of electron ionisation processes. Additionally, the supercell system, which is often constructed for these DFT calculations, suffers from the finite-size effect.⁷⁷ For defect calculations, the defect is periodically repeated causing spurious interactions between defects in adjacent supercells, which requires the system either to be corrected⁷⁸ for these interactions or to be scaled up so that the interactions are minimised due to distance. Because PBC-based methods become more expensive as the size of simulation cell increases, there is a substantial trade-off between choosing more accurate exchange–correlation functionals against scaling up the system.

In this work, we use a hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster method for a systematic investigation of the electronic properties and thermodynamic stability of intrinsic point defects in AlN. This method allows us to study point defect properties at the hybrid DFT level in a simulation cell, which when including the MM region contains about 10 000 atoms, 100 times larger than the 96-atom supercell used in the most recent AlN defect study using periodic boundary conditions.⁷⁵ The defect is embedded in a central region which is treated quantum mechanically (QM) and an outer region where short-range and long-range interactions and polarization effects are modelled using an interatomic potential, or molecular mechanics (MM) approach. In this way, for defects at different charge states, the ionisation process has an unambiguous reference, or vacuum level, which improves the reliability of the calculations. There have been previous studies all using PBC-based techniques, but they either used a lower-level DFT theory or covered an incomplete set of defect types,^{61,64,73,75,79} making it difficult to draw conclusions about the stability of each defect species. Here, by accessing a wider range of hybrid functionals of different HF exchange fractions, we present a comprehensive survey based on QM/MM methods with a critical examination of the effect of different DFT functionals.

While systematically investigating the geometries and electronic structures of all types of intrinsic defects in all their stable charge states, we are able to explore the novel structures and stability of the symmetry-breaking reconstruction of the octahedral-centred N interstitial and the N split-interstitial, both of which are found to be stable; the low formation energies of the Al interstitial, which were found to be comparable to the N vacancy; and the diffuse states of compact defect species, which indicate that the N vacancy and Al interstitial are a possible source of intrinsic shallow donor species. Our results using the B97-2, PBE0, and BB1K functionals show that, even at the hybrid level, different DFT theories can lead to different conclusions for the defect stability. In the final section, equilibrium defect concentrations are computed at different annealing temperatures up to 2800 K to study defect and charge species compensation in scenarios of undoped, n-type doped, and p-type doped materials at all relevant preparative conditions, highlighting the significant roles of the N vacancy, Al vacancy, and Al interstitial for the material's conductivity. Our previous study using the classical force field and the Mott–Littleton embedded cluster approach has already shown some important properties of defect stability, including the identification of the new lower-symmetric, lower-energy N interstitial configuration, and in-lattice migration mechanisms in their formal charges.⁸⁰ In this study, we build on this earlier work by using the QM/MM approach. We focus on key properties of charged intrinsic point defect types in wide-gap AlN.

2 Computational approach

2.1 Hybrid QM/MM embedded cluster technique

The hybrid QM/MM embedded cluster method is applied to model point defect formation in AlN using ChemShell (Tcl version^81,82 and Python version⁸³). There are five concentric regions/layers in our QM/MM system (see the schematic illustration of regions 1–5 in Fig. S1, ESI†). In the inner-most central region (region 1), where the localised defect state is positioned, DFT calculations are performed using NWChem code.⁸⁴ An 86-atom cluster (∼5 Å radius) is chosen to achieve computational efficiency and is sufficient to obtain accurate results (similar size QM clusters have been chosen in previous reports^85–87). The QM cluster is modelled with three different DFT hybrid functionals, B97-2,⁸⁸ PBE0,^89,90 and BB1K⁹¹ (21%, 25%, and 42% exact Hartree–Fock (HF) exchange fractions, respectively). Table S1 (ESI)† shows all three functionals can reproduce structural parameters and optimized bond lengths from previous experiment measurements and calculations. Later in the study, we compare how the three different hybrid functionals could affect the local environment of defect species. The Def2-TZVP⁹² basis set is used for both Al and N. The outermost diffuse, as well as high angular momentum basis functions, are removed for both species (f for N, and d & f for Al), as deleting them results in a very small change in defect formation energies and saves a significant amount of computational time.

Surrounding the central region, there is an “active” MM region (region 3), where all the atoms and shells are explicitly controlled by the classical force field model and the shell model for the polarization of the nitrogen ions.^80,93 In the subsequent “frozen” MM region (region 4), the positions of atoms are fixed. The combined thickness of both active and frozen MM regions is ∼15 Å (containing approximately 10 000 atoms in total). The force field model used in both MM regions is the two-body interatomic potential fitted empirically to a broad range of structural and physical properties of bulk wurtzite AlN. We have successfully implemented the same potential model to predict ionisation potentials, defect formation energies, and defect migration mechanisms of AlN previously.⁸⁰ All the MM calculations in our QM/MM system are performed with the GULP code.^94,95 Lastly, at the surface of the whole QM/MM system, point charges are added to reproduce the Madelung potential of the infinite crystal around the defect site (region 5).

At the interface (region 2) between the QM cluster and MM region, a thin layer of cations is associated with a local electron localising pseudopotential, which also compensates for force mismatches between the QM region and the MM region. The pseudopotential takes the form of a linear combination of three Gaussian functions with respect to the distance r between electrons and respective nuclei:

r²U_p(r) = A₁r exp(−Z₁r²) + A₂r² exp(−Z₂r²) + A₃r² exp(−Z₃r²) (1)

The A and Z parameters are fitted separately using the fit_my_ecp⁹⁶ software, where a global search is done by fulfilling two objectives: (i) minimising the gradient difference between the atoms in the QM, interface and active MM regions; (ii) minimising energy scatter of innermost localised states on anions. In general, the pseudopotentials need to be adjusted to work with the force field implemented in the technique. The resultant pseudopotential for the cations at the interface takes the form with coefficients given in atomic units:

r²U_p(r) = −36r exp(−25r²) + 42.6r² exp(−3.4r²) + 0.45r² exp(−0.55r²) (2)

During a QM/MM calculation with a polarisable shell model potential, ChemShell alternates between QM and MM calculations iteratively to reach ultimate convergence of the whole system. A number of reports have been published using the same method for GaN^85,97,98 and for other ionic materials.^87,99–103 A more detailed discussion of our QM/MM implementation can be found in a previous publication.¹⁰⁴

2.2 Calculations of defect formation energies

The processes of forming point defects in bulk AlN can be represented by the following reactions (using a slightly revised form of Kröger–Vink convention¹⁰⁵) depending on two extremes of the chemical conditions, namely, the “N-rich/Al-poor” condition and “N-poor/Al-rich” condition (e⁻ is an electron and h⁺ is a hole):

Al vacancy (V_Al)

N-poor/Al-rich: Al⁰_Al → V^q_Al − qh⁺ + Al_(s) (4)

N vacancy (V_N)

N-poor/Al-rich: N⁰_N + Al_(s) → V^q_N + qe⁻ + AlN_(s) (6)

Al interstitial (Al_i)

N-poor/Al-rich: Al⁰_Al + Al_(s) → Al^q_i + qe⁻ (8)

N interstitial (N_i)

N-poor/Al-rich: Al⁰_Al + AlN_(s) → N^q_i − qh⁺ + Al_(s) (10)

N antisite (N_Al)

N-rich/Al-poor: Al⁰_Al + N_2(g) → N^q_Al − qh⁺ + AlN_(s) (11)

N-poor/Al-rich: Al⁰_Al + AlN_(s) → N^q_Al − qh⁺ + 2Al_(s) (12)

The formation energy (E_f) of a defect (X) at charge state q is defined as a function of Fermi energy level (E_F) above the valence band maximum (VBM):

where E_d is the energy difference between the total QM/MM energies of the perfect and defective systems. Here, the formation of electrons/holes in the system is determined by the electron chemical potential, whose value typically lies within the gap region (VBM ≤ E_F ≤ CBM). n_i is the number of defect species added to (n_i > 0) or removed from (n_i < 0) the system. The chemical potential of the species, μ_i, is dependent on the thermochemical formation enthalpy of AlN (ΔH_f(AlN) = −3.296 eV):¹⁰⁶

μ_Al + μ_N = ΔH_f(AlN) (14)

In the N-rich/Al-poor condition, μ_Al = ΔH_f(AlN) and μ_N = 0. In the N-poor/Al-rich condition, μ_N = ΔH_f(AlN) and μ_Al = 0. Our formation energies are calculated via the Born–Haber cycles for these defect reactions (see detailed explanation from our previous study⁸⁰). For the intrinsic defects in chemical equilibrium, the source (or the sink) for forming defects can be N₂ gas, solid metallic Al, and solid AlN. Here we calculate the total energies of N₂ gas-phase molecule for the three different functionals using NWChem (PBE0: −109.44 Ha; B97-2: −109.529 Ha; BB1K: −109.521 Ha). The total energies of solid Al are calculated by subtracting the total energies of the Al atom in its gas phase (PBE0: −242.248 Ha; B97-2: −242.367 Ha; BB1K: −242.383 Ha) from the experimental enthalpy of sublimation of solid Al (0.1257 Ha).¹⁰⁷

Since in the “frozen MM region” (region 4) the atomic positions of all the atoms (cores) do not respond to the presence of the defect centre and only shells react to the potential change from the defect, to account for the missing long-range polarisation effect in the frozen MM region and beyond to infinity, a correction term (E_Jost) is applied on any non-neutral charged defect states, where we use the correction formulated by Jost,¹⁰⁸

where R is the radius of the active region (from the cluster centre to the edge of the active MM region), and ε_r = ε₀, the static dielectric constant of the material in the case of an adiabatic process (e.g., the fully relaxed states). Table S2 (ESI)† summarises the Jost correction energies for our defect species from charge −3 to charge +3.

In the case of a vertical process (e.g., electron ionisation and optical excitation) where relaxation involves only electrons, an appropriate form of Jost's correction is applied involving the high-frequency dielectric constant ε_∞. The description of Jost's correction in the case of vertical transition from q to q + Δq thus is:

According to eqn (15) and eqn (16), the accuracy of the Jost correction can depend on the dielectric constants, as well as the size of the simulation cluster, which may vary according to the specific material system. Here we take ε₀ = 10.98 and ε_∞ = 4.68, obtained by averaging over dielectric tensor components calculated in GULP.⁸⁰ As the QM cluster is a part of the infinite pre-optimized structure in the MM software environment (modelled using the GULP code), to keep all calculations consistent, it is reasonable to use the corresponding dielectric constants calculated by the same force field (which shows good agreement with experiments in previous study⁸⁰).

3 Results and discussion

3.1 Ionisation potentials

The ionisation potential of a material is the minimum energy required to excite an electron to the vacuum state. In semiconductor materials, this energy corresponds to the electron transition from the top of valence band (valence band maximum, VBM) to the vacuum level.¹⁰⁹ With the QM/MM embedded cluster method, we can compute the ‘bulk’ value, that is the ionisation potential neglecting any particular surface effects,¹¹⁰ by taking the energy difference between an ionized cluster (losing an electron and allowing electronic degrees of freedom to relax, E^ionized) and a perfect cluster (E^perfect):

E^IP = E^ionized − E^perfect + E^vertical_Jost (17)

E ^vertical_Jost (−0.38 eV) represents the polarisation correction of the transition from 0 to +1 charge state. As shown in Table 1, our calculated ionisation potentials (the negative of which is the VBM level relative to vacuum) are overall deeper than previous values obtained from surface/interface structures using PBC-based methods,^112,113 which are limited by the ambiguity of their vacuum level and finite-size effects. The ionisation potential obtained from our force field model⁸⁰ is slightly smaller (so that the VBM is closer to the vacuum level) than our QM/MM levels, but the discrepancy of at most 1 eV indicates that the force field model provides a useful initial estimate. In addition, the B97-2 VBM level is higher than that calculated using BB1K, which can be attributed to the different electronic localisation performance of the two hybrid functionals. A similar trend is also found in GaN⁹⁸ using the same QM/MM method (c.f.Table 1). Importantly, we will see in later sections that the different electronic localization characteristics from different functionals and their different resultant VBM levels have a strong impact on defect formation energies and charge transition energies, as they are key terms determining these quantities.

Table 1 Ionisation potentials (eV) calculated using our QM/MM method, classical force field (FF)⁸⁰ and comparison with GaN levels calculated using the same QM/MM method¹¹¹ and previous DFT values^112,113

	QM/MM			FF	Prev. DFT
	B97-2	PBE0	BB1K
AlN	7.63	7.98	8.38	7.35	5.8–7.5
GaN	6.63	N/A	7.34		5.7–6.7

---

To the best of our knowledge, there is no experimental ionisation potential reported for AlN, but we can compare our calculated levels indirectly with measured electron affinities. We should note that electron affinities of AlN cannot be obtained directly using our QM/MM method, as high electron delocalisation of the conduction band states cannot be modelled using the embedded cluster method. Here, we obtain an estimation of the electron affinity of AlN by subtracting the experimental bandgap (6.2 eV) from the VBM levels. The electron affinities of AlN are: 1.23 eV (B97-2), 1.78 eV (PBE0), and 2.18 eV (BB1K). No negative electron affinity (NEA) is found for the pure AlN cluster from our calculations,¹¹⁴ which is in good agreement with experimental observations where electron affinities are measured to be between 0.3–2.6 eV.^115–119 The large variation of the experimental CBM values is mainly due to the surface composition, orientation, and dipole moment of the measured samples, as proved by many theoretical studies.^65,120,121 The experimentally measured CBM is affected by either upwards or downwards surface bending depending on the surface condition, so it is expected that their values differ from our “bulk ionisation potential”⁸⁶ as discussed by Zhang et al.¹¹⁰

3.2 Formation of intrinsic point defects

Fig. 1 presents our calculated formation energies for the intrinsic defects, and their charge transition levels in the band gap region for N-rich and Al-rich conditions. We investigated the N and Al vacancies (V_N and V_Al), Al interstitial (Al_i), N interstitial, and N antisite (N_Al). For the N interstitial, two different configurations are investigated explicitly: the “split-interstitial” configuration (N_i-split) and the interstitial located close to the centre of the octahedral chamber of the wurtzite lattice (N_i-oct). The experimental band gap of 6.2 eV is used. The Al antisite defect is not considered for three reasons: first, the Al³⁺ ion is three times smaller than the N³⁻ ion; second, Al is not likely to be stable in an anionic state; third, the B97-2 formation energy of the neutral Al antisite defect calculated in ChemShell is around 14.21 eV (Al-rich condition), much higher than all our considered defect species.

Fig. 1 The defect transition diagram showing formation energies of intrinsic defects in AlN as a function of Fermi energy level (functionals used, from left to right: B97-2, PBE0, and BB1K). The diagrams start from 1 eV below the valence band (VB) maximum (x = 0) on the left of the x-axis and end at 1 eV above the conduction band (CB) minimum on the right of the x-axis. The slopes of the lines indicate different charge states of the defect (see examples in the top, left subplot). The dots represent the transition levels where the defect charge changes.

From the results calculated using B97-2 and PBE0, the more energetically favourable intrinsic point defects are V_N, N_i-split, and V_Al in N-rich condition (N_i-split is excluded in Al-rich condition) for Fermi energies within the band gap region. These more stable defect species are in agreement with a previous study using the supercell approach and the HSE functional.⁷⁵ However, we find Al_i becomes more favourable than the most favourable V_N as the proportion of HF exchange in the functional increases. For Fermi energies in the lower half of the band gap, the formation energy of Al_i is only 0.20 eV higher than V_N using the PBE0 functional, and lower by 0.87 eV using the BB1K functional. An early report found Al_i had the lowest formation energy in zinc blende AlN in the band gap region closer to VBM using the LDA density functional,⁶¹ but the authors state its formation energy is over 2 eV higher in the wurtzite structure counterpart. Zhang et al. calculated the formation energy of Al_i in the wurtzite structure, which they found to be significantly higher than that of V_N using the GGA-level PBE functional.⁶⁴ However, later calculation⁷⁵ using the hybrid HSE functional with spin polarisation taken into account found the formation energy of Al_i almost negligibly higher than that of V_N. In our calculation, a higher fraction of HF exchange deepens the VBM level (Table 1), which in turn increases the energies needed for electrons to localize on (or ionise from) the charged defects and significantly changes their formation energies.

Alkauskas et al. have identified the problem of formation energies of charged defect species in relation to the shifts of band edge levels caused by different levels of theory.¹²² They proposed a band edge correction scheme, where inaccurate band edge energy levels from lower-level DFT theory simulations are shifted according to the more accurate band edge energies calculated with a hybrid functional. Inspired by the theory and the later implementation in GaN,¹²³ another GGA-PBE simulation, with band edges corrected according to HSE energy levels, finds Al_i formation energy lower than V_N in the lower band gap region,⁷⁹ similar to our BB1K results. Nevertheless, our calculations and the wide range of previous results show the predicted stability of Al_i is closely dependent on the level of theory used and the band edge energy levels.

Although our results on the relative stability of each defect species are consistent with previous calculations, our formation energies of positively and negatively charged defect states are systematically lower. Taking V³⁺_N as an example (summarised in Table S3, ESI†), our highest formation energy (from the B97-2 functional) is 3.81 eV lower than the lowest value from previous simulations at the VBM level. One reason could be the choice of DFT functional. There is a significant energy difference among all the formation energies of V³⁺_N calculated at different DFT levels (Table S3†), with the lowest value obtained from simulations using a band-edge-corrected scheme at the GGA-PBE level.

Furthermore, our formation energies of charged defects may differ from values calculated using the PBC-based methods, as our method employs a different treatment of polarization and ionization processes. Nevertheless, both factors should not result in discrepancies in the energies of the neutral charge states. To illustrate this hypothesis, the formation energies of V⁰_Al calculated in our QM/MM environment using different density functionals from the LDA level are compared with previous calculations (summarised in Table S4 (ESI)†). Our QM/MM method reproduces the results of previous reports at the LDA and the GGA level within reasonable accuracy. At the hybrid level, our BB1K formation energy of the V⁰_Al shows the best agreement compared with previous HSE results, and our B97-2 and PBE0 counterparts are closer to previous GGA-PBE results (we have not directly compared with the HSE functional results, as the functional is not currently implemented in NWChem).

When comparing defect formation energies between our three functionals (taking vacancies and interstitials as examples, shown in Fig. S2, ESI†), BB1K results are generally lower than B97-2 and PBE0 results for positively charged defects, and higher for neutral and negatively charged defects, except for V¹⁺_Al, for which the BB1K formation energy is the higher. For neutral defects, the formation energies of most defect types are similar across the three functionals, with one exception in the case of V⁰_Al. For the negative charge states, all formation energies calculated with BB1K surpass significantly those calculated with B97-2 and PBE0. Accurate electronic localization performance is crucial in the calculation of thermochemical formation energies, which can provide other important information such as ionisation potential energies. We have shown that our calculated different ionisation potential energies (Table 1) of electron/hole localisation from the vacuum level are key terms in the formation energies of all charged species. Among the three functionals, BB1K clearly simulates the most stable hole state as it calculates the largest ionization energy, but the B97-2 ionisation energy is also large in comparison to previous DFT simulations. These deeper VBM levels are the primary reason for our discrepancies to the results from previous studies.

A direct consequence of our lower formation energies of charged defects is the shift to higher energies relative to the VBM of all our charge transition levels.⁷⁵ Table S5 (ESI)† summarises the charge transition energies of all the defect species discussed and the comparison with previous calculations. All our charge transition energy levels are closer to the CBM than previous calculated levels obtained using PBC-based methods. From previous calculations, such a shift to the CBM for the energy levels can also be seen when a “higher-rung” (more accurate) DFT functional is used (except for V_N, as the direct transition of (+3|+1) is found in the lower-rung DFT calculations), which arises from the improved band edge levels. These transition levels are important for the optical transition processes, as they are often linked to experimental measurements. Previously deep acceptor V_Al(−2|−3) was assigned to optical emission signals at 2.78 and 3.2 eV,²⁰ and V_Al(−1|−2) to the 2 eV signal.¹²⁴ The statements from these two examples were drawn from previous DFT simulations ranging from LDA to hybrid levels of theories. From our calculations, the 2.78/3.2 eV signal can be related to a deep level of V_Al(−1|0), and the 2 eV signal to acceptor levels of V_Al(−2|−3). However, these transition levels represent charge transitions between fully relaxed defect structures and optical transition occurs much quicker than relaxation. Therefore, it is more suitable to associate experimental signals to calculated vertical transition energies between different charge states using configuration-coordination diagrams. Nevertheless, our simulations highlight the importance of ionization energy level for accurate defect energetics.

3.3 Defect structural and electronic properties

We now discuss the structural and electronic properties of each defect type in all their charge states and compare our results with available simulations and experiments. In particular, the localisation of electrons and holes at each charge state will be explicitly discussed, which, for AlN systems, have only been reported for very few intrinsic defect types^56,61,69 without a systematic investigation. All our simulations start by creating defect species in their formal charges and adding/removing electrons. Most of our results reported below, including the defect structures and electronic information, are qualitatively consistent for all three functionals, unless otherwise stated.

3.3.1 N vacancy (V_N). All seven charge states (+3 to −3) can be structurally stable in AlN. For B97-2 and PBE0, four different charge states (+3 to 0) are stable within the band gap for V_N, whereas the V⁰_N is not favourable within the gap region for BB1K (orange in Fig. 1). We find that V_N behaves more like a donor, an observation agreeing with other theoretical results about the donor/acceptor behaviour of the vacancy. If an N³⁻ ion is removed from the lattice, the surrounding Al ions have no electron density in their valence 3s and 3p orbitals, and V³⁺_N pushes the three neighbouring basal Al ions outwards by 21% and the axial neighbouring Al ion by 26% of the perfect AlN bond lengths. If an electron is added to the system, it is localised at the vacancy site instead of at the surrounding cations. As shown in Fig. 2a, the trapped electron forms a hydrogen-like s orbital state around the centre of the vacancy site; the defect is clearly an “F-centre”, where stable electronic states are trapped by the vacancy defect. The F-centre can trap one more electron, forming V¹⁺_N, where it forms a closed-shell configuration by coupling two electrons. The third additional electron cannot localize at the F-centre; instead, the population of the new spin is shared between the two surrounding Al–Al bonds evenly (Fig. 2b). In the neutral charge state, the spin density wraps around the vacancy centre. This electronic structure is formed by the outermost diffuse atomic orbitals of the surrounding Al³⁺ ions. From experiments, the existence of F-centre type defects has been linked to V_N in AlN as a donor species using the electron paramagnetic resonance (EPR) and the electron-nuclear double resonance (ENDOR) techniques.^12–17 Our results support such an assignment.

Fig. 2 N vacancy in charge states (a) +2 and (b) 0 in AlN. Electron localisation is shown with spin density of (yellow isosurfaces of 0.01 a.u., 0.0075 a.u., and 0.005 a.u.).

Although the local geometry of differently charged vacancies is in reasonable agreement with previous simulations, we did not find the stability of the negatively charged V_N within the band gap, contrary to almost all previous simulations. As more electrons are added to the defect, the surrounding Al³⁺ ions gradually approach the defect centre, leading to a reduction in the bond lengths. When the charge state becomes more negative, the V_N–Al bond lengths are shorter than those in the perfect lattice. At this point, the vacancy defect starts to attract the surrounding positively charged Al ions, due to the negative central Madelung potential. However, as shown in Fig. 1, the charge transition levels for the negatively charged V_N are within the conduction band, indicating auto-ionization and instability for these species.

3.3.2 Al vacancy (V_Al). V_Al can also stabilize structurally in seven charge states (−3 to +3) as an acceptor species, and all charge states are favourable within the band gap for all three functionals (blue in Fig. 1). Our finding is in contrast with the previous work where it was found that the positively charged V_Al was not stable in the gap region^61,62,64 or stable in charge state +1,^{68,71,73–75,79} not in charge states +2 and +3. As V³⁻_Al is created, the surrounding N ions repel each other with an excess of electrons localised on the lone pairs. While the charge state of V_Al becomes less negative, we observe that the holes are localized on the surrounding N ions. Our results show that it is more energetically favourable that the first three electrons are lost successively from the three N ions on the same basal plane (Fig. 3a–c) and the fourth electron from the axial N ion (Fig. 3d). With four hole states localised at the four N ions separately, a spin quadruplet is formed. For each localised hole, the corresponding N ion is pushed further outwards by 2–4% compared to its respective AlN distance between N ion and the V³⁻_Al site centre. Upon further ionisation, the fifth and sixth electrons are lost in the same sequence from the basal N ions (Fig. 3e and f). The presence of negatively charged V_Al in AlN was identified by the positron annihilation spectroscopy technique in previous experiments.¹⁸ The defect species has been linked to optical bands from spectroscopy signals by matching DFT results.^20,125

Fig. 3 Al vacancy in different charge states: (a–f) correspond to −2 to +3 charge states respectively. The spin density (yellow isosurfaces of 0.1 a.u., 0.05 a.u., and 0.025 a.u.) shows hole(s) localisation.

3.3.3 Al interstitial (Al_i). We found Al_i stabilizes at the centre of the c-channel in the wurtzite lattice (i.e. the octahedral void). Our QM/MM calculations find that Al_i at the tetrahedral site has a higher formation energy than that at the octahedral site, consistent with a previous DFT investigation⁶¹ and our prediction using classical force field techniques.⁸⁰ Al_i is structurally stable in four states from +3 to 0, but it is more favourable in charge +3 and +1 within the band gap (green in Fig. 1). The electron can still localise at Al²⁺_i forming an sp-hybrid orbital, which can contribute to an optical transition state. As the positive charge decreases, the defect starts to shift out of the perfect centre, and the surrounding ions are displaced in response. Al²⁺_i and Al¹⁺_i can be the source of a deep donor optical state. This conclusion agrees with previous computational reports,⁷⁵ and there is one experimental report assigning the Al_i²⁺ as a deep donor by EPR,¹²⁶ but no optical band has yet been assigned to it.

3.3.4 N interstitial. Previous DFT calculations have reported two common configurations of N interstitials in AlN, the octahedral-centric interstitial (N_i-oct) and the split-interstitial (N_i-split). In our simulations, we find three stable configurations of the N interstitial defect: an octahedral-centric interstitial (N_i-oct(cen)), a lower-symmetry reconstructing octahedral interstitial (N_i-oct(recon)), and a split interstitial (N_i-split). These configurations, as shown in Fig. 4 and discussed in detail below, show the very complex nature of N interstitial defect. In the first configuration (Fig. 4a–c), a nitrogen ion stabilises at the centre of the octahedral void of the wurtzite structure, around the same position as the Al_i discussed above. As with the Al_i, we did not find the tetrahedral-centric N interstitial to be stable. In our previous defect study using force field method, we identified a new stable interstitial configuration where the interstitial N³⁻ ion is displaced from the ideal octahedral-centric site towards one of the three nearest neighbouring on-lattice N³⁻ ions (termed as “interstitialcy configuration” in our previous study⁸⁰), which has not been reported before in III–V nitride systems. Here, by using our QM/MM DFT method, we also obtain the same stable configuration (N_i-oct(recon) in this work) and we now continue investigating its geometry and electronic structures in different charge states.

Fig. 4 The configurations and the spin densities of (a–c) N octahedral-centric interstitial defect in charge states −2, −1, and 0 respectively, (d–f) N octahedral-reconstruction interstitial defect in charge states −2, −1, and 0 respectively, and (g–i) N split-interstitial defect in charge states 0, +1, and +2 respectively (yellow isosurfaces are 0.1 a.u., 0.075 a.u., and 0.05 a.u.).

Fig. S3 (ESI)† summarises the stability of three types of N interstitials. Both N_i-oct(cen) and N_i-oct(recon) are structurally and energetically stable in four different charge states (−3 to 0) in the band gap. In the upper gap region, N_i-oct(recon) is more thermodynamically favourable than N_i-oct(cen) by 0.4 eV (for charge −3) and 0.74 eV (for charge −2) from B97-2 simulations (0.35 eV/0.64 eV from PBE0 and 0.39 eV/0.60 eV from BB1K). As the Fermi level shifts towards the mid-gap region, the energies of the two defects become closer in charge state −1 and N_i-oct(cen) becomes slightly more favourable than N_i-oct(recon) for the neutral charge state. On removing the first electron from charge state −3, the hole state localises on the introduced N species for N_i-oct(cen), whereas for N_i-oct(recon) the population of the spin is shared between the defect and the nudged-off N ion (Fig. 4d). The lower-symmetry configuration and delocalisation of the hole in the latter case can stabilise the local environment further around the defect. However, when ionising the second and the third electron local configuration, hole localisation can only occur on the introduced N ion (even if we forcefully ionise the nudged-off N ion), as shown in Fig. 4e and f. This factor could explain the largest energy difference found in the −2 charge state. The special case of the spin delocalization N²⁻_i-oct(recon) indicates the electronically unstable nature of negative N ions in the wurtzite octahedral void. Thus, the two defect types are considered as interchangeable cases for the same N_i-oct defect and their combined lowest energies are presented in the single red line in Fig. 1.

Removing three electrons from the 2p orbitals of the N³⁻ ion in turn forms corresponding high spin states (up to 3 lone spins) of the defect. On further ionisation, we cannot find any stable configuration of the N_i-oct in any of the positive charge states, and further relaxation forces the N defect to bind with another on-site N ion (the split-interstitial). Compared with all other types of defect states in the band gap, forming N_i-oct requires the highest formation energy. Therefore, N_i-oct is predicted to be a high energy species in AlN, so it is unlikely to be present in significant concentrations.

Fig. 4g–i present the configuration and spin densities for the N split-interstitial (N_i-split), where the interstitial N ion is bonded with another on-site N ion forming an N–N dimer structure, which is the most stable interstitial configuration across the whole gap region (see Fig. 1 and 4j–o). The species are stable in the charge states from +3 to −2 in the band gap. Some studies on this defect type did not find +3 charge states stable in cubic⁶¹ or wurtzite^69,75 AlN, while others did^73,79 (Chen et al. did not find the +2 and +1 states stable⁷⁹). In the +3 charge state, the bond length of the N–N dimer is close to the bond length of the N₂ gas molecule (1.09 Å). As more electrons are trapped by the defect, the N–N dimer is stretched by 6% (+2 charge state), 15% (+1), 24% (0), and 33% (−1). These geometries are found to be in good agreement with previous hybrid DFT studies.^69,73 Meanwhile, the N–N dimer rotates around the same N atom site of the perfect AlN lattice along the (110) lattice plane, as the charge decreases. Going from the charge state of −1 to −2, no significant change is found in the bond length and N–N orientation. The added electron is delocalised across the whole QM cluster region, resulting in bond breakage between the two nitrogen atoms. The formation energy of N²⁻_i-split is higher than N²⁻_i-oct by 0.78 eV (from b97-2, 0.83 eV from PBE0, 1.18 eV from BB1K), indicating a possible structural transformation to the latter defect configuration in this charge state.

We find the most stable configuration of N¹⁺_i-split is a triplet (Fig. 4h). In the charge state of 0 or +2, a singlet is found in the most stable configuration (Fig. 4g and i). Furthermore, in the −1 and +3 charge states, the defect is most stable in the closed-shell electronic configuration. The spin densities of different charge states prove that the electronic configuration of the split interstitial is analogous to that of an N₂ molecule. In the charge states of 0, +1, and +2, the spin density is as expected from the occupation of an anti-bonding π orbital. The N²⁻_i-split does not form deep states under CBM, although a stable configuration can be found after geometry optimisation. The N_i-split defect is calculated to be amphoteric (Fig. 1). The optical transition energies are those of a deep donor in the lower half of the band gap region.

3.3.5 N antisite (N_Al). N_Al is found to be stable in three different configurations, two of which are stabilised by the substitutional N bonding with another neighbouring N atom. In the first configuration, the N defect (termed here “config-A”, or N_Al-A) bonds with a basal N atom (Fig. 5a and b). N_Al-A is attracted towards this basal N atom, reducing the N–N bond distance. N_Al-A is calculated to be stable in charge states +3 to −2 in the band gap region. In charge state −2, the N–N bond distance is shorter than the perfect Al–N bond length by 0.57 Å (30%) from B97-2 simulations, close to the 29% value calculated previously for cubic AlN.⁵⁹ As the charge becomes more positive, the N–N bond distance becomes shorter, and the shortest value is reached in the neutral charge state at 1.17 Å from B97-2 simulations (1.14 Å using BB1K), close to the bond length of an N₂ gas molecule. The spin density in charge states −2 and −1 accords with an anti-bonding π orbital (Fig. 5a and b), which shows the electronic similarity to an N₂ molecule.

Fig. 5 Optimized geometry of N antisite defect: “config-A” in the charge states (a) −2 and (b) −1, and “config-B” in the charge states of (c) +1 and (d) +3. Electron localisation is shown with spin density (yellow isosurfaces of 0.1 a.u., 0.05 a.u., and 0.025 a.u.).

We found another structurally stable configuration for charge state −1 N_Al, where the antisite defect can bond with the axial N atom. This configuration is less stable than the N_Al-A in the same charge state, but it can transform into another antisite configuration (“config-B”, or N_Al-B) as more electrons are removed from the point defect. We found N_Al-B is stable in charge states from +3 to −1 in the band gap, and it is more stable than N_Al-A in all positive charge states from +1 to +3. In the neutral charge state, the N axial dimer moves towards the two closest basal N atoms and reconstructs itself into a stable trigonal configuration (Fig. 5c and d). In charge state +1, the spin density shows a p-orbital-like characteristic (Fig. 5c) and localises at the central N of the trigonal structure. No localised hole is found in the system in charge state +2. In charge state +3, spin density shows a hole localises on the basal N atom outside the trigonal structure.

There have been several theoretical studies of N_Al in AlN, but only a few of them discuss the two configurations of the defect. Most of the early work investigated the defect only in the cubic phase of AlN,^57,59,61 and claimed that there is a small difference between the defect in the two structural phases.⁶¹ The EL2 geometry of the N antisite in cubic AlN in these reports is the closest configuration to our config-A and config-B types. More recent studies of the N_Al are either unclear about the geometry,⁷⁵ or only about the substitutional type,⁷³ which is our proposed third configuration of N_Al. However, the substitutional N_Al is generally less stable than the other two configurations, so we do not present a detailed discussion. Experimentally, there is still no identification of the N antisite defect in AlN.

3.4 Diffuse states of intrinsic defects

There are limitations in our QM/MM method, with which it is difficult to calculate a highly diffuse defect state. The diffuse electronic part of the defect will be confined due to the nature of the embedded system, which might cause electron over-localisation and an incorrect electronic exchange–correlation energy. However, with a sufficiently large QM region, useful results can still be obtained on diffuse states, such as the stability of any shallow donor/acceptor states. To complement our results, we consider the possible cases where “compact” defect states can trap one or more electrons/holes in “diffuse” atom-like states of large effective radius, and compute the “attaching energies” of diffuse carriers attached to compact defects. The theory and the calculated attaching energies are summarised in Note S1 and Tables S6–S8 (ESI).†

Table S9–S15 (ESI)† list the calculated formation energies of all the diffuse states of defects in different charge states with the corresponding compact states, where their formation energies can be compared directly, to determine which type is more favourable in each charge states. From these results, we found most defects have lower energies in their compact state. V_Al, N_i-oct, N_i-split, N_Al-A, and N_Al-B show no preference for the diffuse states in any of their charge states, as their compact state energies are always lower than the diffuse variants. For Al_i and V_N, the energies of diffuse states are lower than the compact neutral charge state (Tables S9 and S10, ESI†). Calculations using B97-2 and PBE0 functionals show that Al_i⁺ with a diffusely-bound electron has a lower energy than the compact Al⁰_i (BB1K calculation shows the diffuse states of one, two, and three electrons are all more favourable than compact Al⁰_i), which is a direct consequence of the higher formation energy of Al⁰_i at the CBM (Fig. 1). The favourable diffuse states can also be seen in other charge states of Al_i from the BB1K calculations, while the other functionals indicate otherwise. For V_N, only the BB1K calculations predict the defect to be more favourable than its diffuse state of one electron compared with V⁰_N. The BB1K calculation shows another lower energy diffuse electron state for V¹⁻_N, while PBE0 calculation shows the diffuse state is more favourable by only 0.02 eV and B97-2 calculation shows compact V¹⁻_N is more favourable. In conclusion, from these calculations, we found the diffuse states are more likely to appear in Al interstitials, where Al⁰_i is a possible shallow donor source. The prediction as to the nature of V⁰_N and V¹⁻_N varies depending on the choice of hybrid functional, which has a large impact on defect property prediction.

3.5 Defect and charge carrier concentrations

From our defect formation energies, we can calculate the equilibrium concentrations of defect species and charge carriers. At a given temperature, the equilibrium concentration of a defect can be determined by its defect formation energy employing a Boltzmann distribution. Meanwhile, the equilibrium concentrations of the charge carriers (i.e., [n₀] for electrons and [p₀] for holes) in a semiconductor material can be calculated by integrating the electronic density of states with the Fermi–Dirac distribution formula. The Fermi energy (E_F), which is a common term in both concentrations, can be determined by a self-consistent method in the condition of charge neutrality in the system:

[n₀] + [acceptors] = [p₀] + [donors] (18)

Here, we apply the code SC-FERMI¹²⁷ to compute the self-consistent Fermi energies and the corresponding defect and charge carrier concentrations in the range of temperature 300 K ≤ T ≤ 2800 K (covering common temperatures of bulk AlN crystal growth using physical vapour transport techniques^3,128–130). The density of states of AlN is calculated by using the Vienna Ab initio Simulation Package (VASP)^131–133 and PBE0 functional (more detail of the computational method is presented in Note S2, ESI†), as B97-2 and BB1K functionals are not implemented in VASP. We obtained a band gap energy of 6.11 eV, close to the experimental gap value. Our defect concentration results are therefore calculated using only the PBE0 defect formation energies (Fig. 1).

As shown in Fig. 6 (left), in N-rich conditions, the Fermi level lies in the mid-gap region, slightly varying from 3.16 to 3.66 eV as the temperature increases. The highest equilibrium concentration among all the intrinsic defects is N_i-split and V_N (the latter only about half an order of magnitude lower), indicating V_N and N_i-split are the two dominant intrinsic defect species in AlN under N-rich chemical conditions. In the range of the self-consistent Fermi energy, V¹⁺_N and N¹⁺_i,split are the dominant charge states. In N-poor conditions (Fig. 6, right), due to the lower formation energy of V_N, the equilibrium concentration of V_N is significantly higher than all other defects. At 2500 K, the V_N concentration is higher by 3 orders of magnitude than its counterpart in N-rich conditions. The concentrations of holes and other defect species are much lower due to the rise of the equilibrium Fermi energy. The Fermi energy ranges from 4.73 to 5.16 eV, closer to the CBM. In turn, the electron concentration (overlapped with the green curve in Fig. 6) increases. Close to the Fermi level, V¹⁺_N is still the dominant charge state in N-poor conditions. This indicates the negative charge carrier is compensated by the positive charge V¹⁺_N in N-poor chemical condition to maintain charge neutrality in the material.

Fig. 6 The calculated self-consistent Fermi energy and equilibrium concentrations of charge carriers and defect species as a function of temperature.

As the defect formation energy is a function of chemical potential μ, which can be determined by the temperature and the partial pressure of the system, it is possible to compare our results directly to experiment. Here we follow the guidelines developed by Reuter and Scheffler:¹³⁴ the N chemical potential (μ_N) as a function of temperature (T) and N partial pressure (P) can be expressed as:

in which

The zero state (μ_N(0 K, P⁰) = 1 eV) is defined when P⁰ = 1 atm. The enthalpy (H_N2) and the entropy (S_N2) of N₂ at different temperatures are obtained from the thermochemical experiment database.¹³⁵ The μ_N(T, P) with respect to temperature and N₂ partial pressure is summarised in Table S16 and Fig. S4 (ESI).† Given the condition of −3.296 eV ≤ μ_N ≤ 0 eV for forming stable AlN, all the defect formation energies and, hence, the equilibrium concentrations of all the defect species can be determined by a specifying temperature and pressure. The equilibrium concentrations of defect species and charge carriers with respect to the N₂ partial pressure (10⁻² to 10² atm) at 2400 K are presented in Fig. 7 (only the four highest defect species are shown). Our results indicate that regardless of the N₂ partial pressure (during crystal growth), V_N remains the dominant defect species in pristine bulk AlN. As the concentration of electron carriers is several orders of magnitude higher than that of hole carriers, AlN exhibits n-type semiconductor characteristics. There appears to be no direct experimental evidence as to the intrinsic n-type nature of AlN (although it is nominally believed to be the case¹⁶), but n-type conductivity in undoped GaN has been observed and the carrier concentration is reported¹³⁶ to be of the order of 10¹⁸ cm⁻³. In the same report, V_N is identified as the primary native donor accounting for its n-type conductivity. Our result is qualitatively consistent with the conclusion of the GaN experiments and predicts the dominant concentration of V_N results in intrinsic n-type conductivity in AlN.

Fig. 7 The calculated equilibrium concentrations of charge carriers and defect species as a function of N₂ partial pressure at 2400 K.

To investigate further the self-compensation characteristics in the material, a fully ionised donor or acceptor is introduced to determine how the equilibrium Fermi energy and species concentrations respond in the system. Our procedure allows us to add a fixed concentration of mass-less charged species to bias the charge neutrality of the equilibrium system without introducing any extrinsic defect species.¹²⁷ In response to the bias, Fermi energies and concentrations of intrinsic defects and charge carriers are then readjusted self-consistently to meet the condition of charge neutrality. In the case of donors, we introduce a fix-concentration (1 × 10¹⁹ cm⁻³) of a fully-ionised donor (represented by a +1 charge species) to represent the n-type doping scenario. In both N-rich and N-poor conditions (see Fig. 8a and b), the negative charge carrier concentration increases significantly due to the presence of the ionised donor. In N-rich conditions, the ionised donor is primarily compensated by the negative charge V_Al (3 × 10¹⁸ cm⁻³) due to its low formation energy (Fig. 1) in the range of the equilibrium Fermi level of 5.00–5.76 eV above the VBM. At higher temperatures, the positive charge V_N starts to contribute more to the compensation, but it still leaves a significant number of uncompensated donors. Hence, n-type conductivity is restricted by the presence of negatively charged V_Al. In N-poor conditions, our calculations show that the concentration of electrons remains at the level of 10¹⁹ cm⁻³ and most ionised donors are not compensated. The equilibrium Fermi level decreases from 6.41 eV (0.21 eV above CBM) to 6.12 eV (0.08 eV below CBM), which indicates significant quantities of free electrons. In a previous report of Tuomisto et al. using positron annihilation spectroscopy,¹⁸ the concentration of negative charge V_Al was estimated analytically at about 1 × 10¹⁷ cm⁻³ in bulk AlN crystal (with O impurities concentration at the order of 10¹⁸ cm⁻³, one order of magnitude lower than our fixed donor concentration). Additionally, cation vacancies were identified as the major defect species in the Si-doped n-type Al_xGa_1−xN.⁹ Our results are in good agreement with these reports.

Fig. 8 The calculated self-consistent Fermi energy and equilibrium concentrations of charge carriers and defect species as a function of temperature, as a fully ionised donor (a and b) or a fully ionised acceptor (c and d).

On introducing fixed-concentration (1 × 10¹⁹ cm⁻³) fully-ionised acceptors (represented by a −1 charge species), we found different behaviour. In N-rich conditions (Fig. 8c) the ionised acceptors are compensated primarily by positively charged V_N and Al_i. As the temperature rises, N_i-split starts to make an increasing contribution to the compensation because of the low formation energies of positively charged species. Although the Fermi levels are in the lower region of the band gap, between 1.74–2.36 eV above VBM, the concentration of hole charge carriers (p₀) remains several orders of magnitudes smaller than the introduced acceptors and other defect species. Only at high temperatures, we see a significant concentration of hole carriers. In N-poor condition (Fig. 8d), only V_N contributes significantly to the compensation compared to other defect species. As the Fermi levels are at the upper region of the band gap between 3.14–4.62 eV above VBM, the concentration of electron charge carriers is much higher than hole carriers, leading to weakly n-type semiconductor characteristics. Enhancing p-type conductivity has always been a challenge in Mg-doped p-type AlN, a characteristic that we can also see from our relatively low concentration of p₀ in both chemical conditions. It is widely established that suppressing the formation of V_N can effectively achieve higher p-type conduction in Mg-dope AlN by choosing N-rich growing condition¹⁰ or with a high V/III ratio.¹¹ These experiment reports have not interpreted their findings in terms of Al_i as another important self-compensating species, as it was ruled out by the high Al_i formation energy calculated by early simulations.⁶¹

Since our investigation of fixed concentration charged species did not introduce any new defect species to the system, our discussion did not consider the formation of extrinsic defects, such as oxygen, carbon, etc., and their complex defect structures formed with the intrinsic defects. These impurity-incorporated defect complexes could have an impact on the defect formation energies, the charge transition levels, and, in turn, the equilibrium Fermi level positions. Among all, the Al-vacancy-O-substitution defect complex (V_Al–n(O_N)) is the most discussed defect complex, whose formation energy was calculated (using supercell approaches) to be over 2 eV lower than that of V_Al at the HSE level.^68,75 The lower formation energies could result in more defect species (V_Al in Fig. 8a) compensating the ionised donors and possibly a higher n-type conductivity.

The identification of deep/shallow donor and acceptor states in most experiment reports commonly relies on photoluminescence (PL) signals, with assignments made based on previous theoretical results, specifically the charge transition levels and the configuration-coordinate diagrams. Using the QM/MM embedded cluster method, we have demonstrated a shift of the charge transfer levels towards CBM compared to the levels obtained using supercell methods (refer to the detailed discussion in Section 3.2). Future work will explore both the effects of dopants and optical transition energies.

4 Summary and conclusions

We have systematically investigated the stability and electronic structures of all intrinsic defect species at the dilute limit using our advanced embedded cluster QM/MM approach. While the relative stability of all defect species is consistent with previous DFT studies, our work provides significant new insights. Notably, the aluminium interstitial, Al_i, emerged as the most favourable species in the lower region of the band gap from the hybrid BB1K calculation, indicating a significant presence of Al_i in the material. Furthermore, we confirm the presence of the novel lower-symmetry nitrogen interstitial defect, calculated in our previous study using force field methods, and explore the electronically unstable nature of the defect species. We have found our formation energies of the charged defects to be generally lower than previous values calculated using supercell methods. As a result, shifts towards the CBM were observed in our charge transition energies compared to those from previous calculations. These differences may stem from either the choice of DFT functionals or the nature of our QM/MM method, which omits the finite-size effect from the supercell-based methods, enabling the simulation of accurate ionization potential levels against the true vacuum reference level. In addition to the structural properties of the defects, the electron localisation characteristics of all the defects have been systematically investigated for every charge state. The diffuse states of the defect species are compared with the compact counterparts, highlighting Al⁰_i as a potential shallow donor source. The shallow nature of V_N⁰ and V_N¹⁻ varies depending on the exchange fraction of the hybrid functional.

Concentration calculations indicate that AlN exhibits intrinsic n-type behaviour throughout the temperature range of 300 to 2800 K. With a concentration of 1 × 10¹⁹ cm⁻³ of fully-ionised donor species, a substantial amount of compensating negative charge V_Al (on the order of 10¹⁸ cm⁻³ in N-rich condition) is evident. When the same concentration of ionised acceptors is introduced into the material, the p-type conductivity is hindered by compensation from the positive V_N in both chemical potential limits and from Al_i as a secondary contribution in N-rich condition. While these results reaffirm the dominant role of cation and anion vacancies in determining the defect properties of AlN, we emphasize that cation interstitials could also play a key role. Our implementation of the advanced QM/MM technique and the new insights into the nature of intrinsic defects in AlN will benefit not only the III–V nitride semiconductor science but also a broader spectrum of defect chemistry.

Data availability

The data supporting this article have been included as part of the ESI.† All open-source codes used in the study, including fit_my_ecp and SC-FERMI are referenced in the paper and available on GitHub.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like to acknowledge the usage of several HPC services for the completion of this work. The usage of THOMAS/YOUNG is via the UK Materials and Molecular Modelling Hub for computational resources, MMM Hub, which is partially funded by EPSRC (EP/L000202, EP/T022213). The usage of ARCHER/ARCHER2 UK National Supercomputing Service (http://www.archer2.ac.uk/) is via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202, EP/R029431, EP/T022213, EP/X035859). The co-author Qing Hou acknowledges funding from Shanghai Pujiang Program (grant 22PJ1411300).

References

1. R. Hahn, V. Glaw, A. Ginolas, M. Töpfer, K. Wittke and H. Reichl, Microelectron. Int., 1999, 16, 21–26.
2. T. Takano, T. Mino, J. Sakai, N. Noguchi, K. Tsubaki and H. Hirayama, Appl. Phys. Express, 2017, 10, 31002.
3. M. Kneissl, T. Y. Seong, J. Han and H. Amano, Nat. Photonics, 2019, 13, 233–244.
4. Y. Taniyasu, M. Kasu and T. Makimoto, Nature, 2006, 441, 325–328.
5. Y. Xue, H. Wang, N. Xie, Q. Yang, F. Xu, B. Shen, J. Shi, D. Jiang, X. Dou, T. Yu and B. Sun, J. Phys. Chem. Lett., 2020, 11, 2689–2694.
6. J. R. Weber, W. F. Koehl, J. B. Varley, A. Janotti, B. B. Buckley, C. G. de Walle and D. D. Awschalom, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 8513–8518.
7. W. Werdecker and F. Aldinger, IEEE Trans. Compon., Hybrids, Manuf. Technol., 1984, 7, 399–404.
8. T. Koppe, H. Hofsäss and U. Vetter, J. Lumin., 2016, 178, 267–281.
9. A. Uedono, K. Tenjinbayashi, T. Tsutsui, Y. Shimahara, H. Miyake, K. Hiramatsu, N. Oshima, R. Suzuki and S. Ishibashi, J. Appl. Phys., 2012, 111, 013512.
10. Y. Wu, D. A. Laleyan, Z. Deng, C. Ahn, A. F. Aiello, A. Pandey, X. Liu, P. Wang, K. Sun, E. Ahmadi, Y. Sun, M. Kira, P. K. Bhattacharya, E. Kioupakis and Z. Mi, Adv. Electron. Mater., 2020, 6, 2000337.
11. T. Kinoshita, T. Obata, H. Yanagi and S. Inoue, Appl. Phys. Lett., 2013, 102, 012105.
12. K. Atobe, M. Honda, N. Fukuoka, M. Okada and M. Nakagawa, Jpn. J. Appl. Phys., 1990, 29, 150–152.
13. S. M. Evans, N. C. Giles, L. E. Halliburton, G. A. Slack, S. B. Schujman and L. J. Schowalter, Appl. Phys. Lett., 2006, 88, 62112.
14. S. Nakahata, K. Sogabe, T. Matsuura and A. Yamakawa, J. Am. Ceram. Soc., 1997, 80, 1612–1614.
15. P. M. Mason, H. Przybylinska, G. D. Watkins, W. J. Choyke and G. A. Slack, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1937–1947.
16. V. A. Soltamov, I. V. Ilyin, A. A. Soltamova, E. N. Mokhov and P. G. Baranov, J. Appl. Phys., 2010, 107, 113515.
17. V. A. Soltamov, I. V. Ilyin, A. A. Soltamova, D. O. Tolmachev, N. G. Romanov, A. S. Gurin, V. A. Khramtsov, E. N. Mokhov, Y. N. Makarov, G. V. Mamin, S. B. Orlinskii and P. G. Baranov, Appl. Magn. Reson., 2013, 44, 1139–1165.
18. F. Tuomisto, J.-M. Mäki, T. Yu. Chemekova, Y. N. Makarov, O. V. Avdeev, E. N. Mokhov, A. S. Segal, M. G. Ramm, S. Davis, G. Huminic, H. Helava, M. Bickermann and B. M. Epelbaum, J. Cryst. Growth, 2008, 310, 3998–4001.
19. A. Uedono, S. Ishibashi, S. Keller, C. Moe, P. Cantu, T. M. Katona, D. S. Kamber, Y. Wu, E. Letts, S. A. Newman, S. Nakamura, J. S. Speck, U. K. Mishra, S. P. DenBaars, T. Onuma and S. F. Chichibu, J. Appl. Phys., 2009, 105, 54501.
20. A. Sedhain, L. Du, J. H. Edgar, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2009, 95, 262104.
21. T. Koyama, M. Sugawara, T. Hoshi, A. Uedono, J. F. Kaeding, R. Sharma, S. Nakamura and S. F. Chichibu, Appl. Phys. Lett., 2007, 90, 241914.
22. T. Hoshi, T. Koyama, M. Sugawara, A. Uedono, J. F. Kaeding, R. Sharma, S. Nakamura and S. F. Chichibu, Phys. Status Solidi C, 2008, 5, 2129–2132.
23. B. Bastek, F. Bertram, J. Christen, T. Hempel, A. Dadgar and A. Krost, Appl. Phys. Lett., 2009, 95, 32106.
24. S. F. Chichibu, K. Hazu, Y. Ishikawa, M. Tashiro, T. Ohtomo, K. Furusawa, A. Uedono, S. Mita, J. Xie, R. Collazo and Z. Sitar, Appl. Phys. Lett., 2013, 103, 142103.
25. M. Bickermann, B. M. Epelbaum, O. Filip, P. Heimann, S. Nagata and A. Winnacker, Phys. Status Solidi B, 2009, 246, 1181–1183.
26. A. Sedhain, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2012, 100, 221107.
27. Q. Hu, S. Tanaka, T. Yoneoka and T. Noda, Nucl. Instrum. Methods Phys. Res., Sect. B, 2000, 166–167, 70–74.
28. U. Vetter, S. Müller, M. Brötzmann, H. Hofsäss and J. B. Gruber, Diamond Relat. Mater., 2011, 20, 782–784.
29. J. A. Freitas, J. Cryst. Growth, 2005, 281, 168–182.
30. K. B. Nam, M. L. Nakarmi, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2005, 86, 222108.
31. N. Nepal, M. L. Nakarmi, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2006, 89, 92107.
32. S. Bellucci, A. I. Popov, C. Balasubramanian, G. Cinque, A. Marcelli, I. Karbovnyk, V. Savchyn and N. Krutyak, Radiat. Meas., 2007, 42, 708–711.
33. A. Sedhain, N. Nepal, M. L. Nakarmi, T. M. Al tahtamouni, J. Y. Lin, H. X. Jiang, Z. Gu and J. H. Edgar, Appl. Phys. Lett., 2008, 93, 41905.
34. I. A. Weinstein, A. S. Vokhmintsev and D. M. Spiridonov, Diamond Relat. Mater., 2012, 25, 59–62.
35. L. Shen, N. Wang and X. Xiao, Mater. Lett., 2013, 94, 150–153.
36. Z. Wu, W. Zhang, H. Hu, S. Zuo, F. Wang, P. Yan, J. Wang, R. Zhuo and D. Yan, Mater. Lett., 2014, 136, 95–98.
37. W.-Y. Wang, P. Jin, G.-P. Liu, W. Li, B. Liu, X.-F. Liu and Z.-G. Wang, Chin. Phys. B, 2014, 23, 87810.
38. J. Pastrňák, S. Pačesová and L. Roskovcová, Czech. J. Phys., B, 1974, 24, 1149–1161.
39. R. A. Youngman and J. H. Harris, J. Am. Ceram. Soc., 1990, 73, 3238–3246.
40. M. Benabdesselam, P. Iacconi, D. Lapraz, P. Grosseau and B. Guilhot, J. Phys. Chem., 1995, 99, 10319–10323.
41. Q. Hu, S. Tanaka, T. Yoneoka and V. Grismanovs, Radiat. Eff. Defects Solids, 1999, 147, 283–292.
42. S. Schweizer, U. Rogulis, J.-M. Spaeth, L. Trinkler and B. Berzina, Phys. Status Solidi B, 2000, 219, 171–180.
43. Y. G. Cao, X. L. Chen, Y. C. Lan, J. Y. Li, Y. P. Xu, T. Xu, Q. L. Liu and J. K. Liang, J. Cryst. Growth, 2000, 213, 198–202.
44. B. Berzina, L. Trinkler, J. Sils and E. Palcevskis, Radiat. Eff. Defects Solids, 2001, 156, 241–247.
45. A. Dadgar, A. Krost, J. Christen, B. Bastek, F. Bertram, A. Krtschil, T. Hempel, J. Bläsing, U. Haboeck and A. Hoffmann, J. Cryst. Growth, 2006, 297, 306–310.
46. B. E. Gaddy, Z. Bryan, I. Bryan, J. Xie, R. Dalmau, B. Moody, Y. Kumagai, T. Nagashima, Y. Kubota, T. Kinoshita, A. Koukitu, R. Kirste, Z. Sitar, R. Collazo and D. L. Irving, Appl. Phys. Lett., 2014, 104, 202106.
47. X. Tang, F. Hossain, K. Wongchotigul and M. G. Spencer, Appl. Phys. Lett., 1998, 72, 1501–1503.
48. K. Kornitzer, W. Limmer, K. Thonke, R. Sauer, D. G. Ebling, L. Steinke and K. W. Benz, J. Cryst. Growth, 1999, 201–202, 441–443.
49. S. Tungasmita, P. O. Å. Persson, L. Hultman and J. Birch, J. Appl. Phys., 2002, 91, 3551–3555.
50. T. Nagashima, Y. Kubota, T. Kinoshita, Y. Kumagai, J. Xie, R. Collazo, H. Murakami, H. Okamoto, A. Koukitu and Z. Sitar, Appl. Phys. Express, 2012, 5, 125501.
51. M. Bickermann, B. M. Epelbaum, O. Filip, B. Tautz, P. Heimann and A. Winnacker, Phys. Status Solidi C, 2012, 9, 449–452.
52. K. B. Nam, M. L. Nakarmi, J. Li, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2003, 83, 878–880.
53. M. L. Nakarmi, N. Nepal, J. Y. Lin and H. X. Jiang, Appl. Phys. Lett., 2009, 94, 091903.
54. A. Pandey, X. Liu, Z. Deng, W. J. Shin, D. A. Laleyan, K. Mashooq, E. T. Reid, E. Kioupakis, P. Bhattacharya and Z. Mi, Phys. Rev. Mater., 2019, 3, 053401.
55. P. Wang, D. Wang, S. Mondal, M. Hu, J. Liu and Z. Mi, Semicond. Sci. Technol., 2023, 38, 043002.
56. P. Boguslawski, E. Briggs, T. A. White, M. G. Wensell and J. Bernholc, MRS Online Proc. Libr., 1994, 339, 693–698.
57. T. Mattila and R. M. Nieminen, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 16676–16682.
58. I. Gorczyca, A. Svane and N. E. Christensen, Solid State Commun., 1997, 101, 747–752.
59. I. Gorczyca, A. Svane and N. E. Christensen, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 8147–8157.
60. A. Fara, F. Bernardini and V. Fiorentini, J. Appl. Phys., 1999, 85, 2001–2003.
61. C. Stampfl and C. G. Van de Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 155212.
62. H. Seo, M. Govoni and G. Galli, Sci. Rep., 2016, 6, 20803.
63. Y. Osetsky, M.-H. Du, G. Samolyuk, S. J. Zinkle and E. Zarkadoula, Phys. Rev. Mater., 2022, 6, 094603.
64. Y. Zhang, W. Liu and H. Niu, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 35201.
65. P. G. Moses, M. Miao, Q. Yan and C. G. Van de Walle, J. Chem. Phys., 2011, 134, 84703.
66. L. Silvestri, K. Dunn, S. Prawer and F. Ladouceur, Appl. Phys. Lett., 2011, 99, 122109.
67. L. Gordon, J. L. Lyons, A. Janotti and C. G. de Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 85204.
68. Q. Yan, A. Janotti, M. Scheffler and C. G. Van de Walle, Appl. Phys. Lett., 2014, 105, 111104.
69. A. Szállás, K. Szász, X. T. Trinh, N. T. Son, E. Janzén and A. Gali, J. Appl. Phys., 2014, 116, 113702.
70. J. B. Varley, A. Janotti and C. G. Van de Walle, Phys. Rev. B, 2016, 93, 161201.
71. J. S. Harris, J. N. Baker, B. E. Gaddy, I. Bryan, Z. Bryan, K. J. Mirrielees, P. Reddy, R. Collazo, Z. Sitar and D. L. Irving, Appl. Phys. Lett., 2018, 112, 152101.
72. D. Alden, J. S. Harris, Z. Bryan, J. N. Baker, P. Reddy, S. Mita, G. Callsen, A. Hoffmann, D. L. Irving, R. Collazo and Z. Sitar, Phys. Rev. Appl., 2018, 9, 54036.
73. Y. Gao, D. Sun, X. Jiang and J. Zhao, J. Appl. Phys., 2019, 125, 215705.
74. P. C. Bowes, Y. Wu, J. N. Baker, J. S. Harris and D. L. Irving, Appl. Phys. Lett., 2019, 115, 52101.
75. I. A. Aleksandrov and K. S. Zhuravlev, J. Phys.: Condens. Matter, 2020, 32, 435501.
76. A. Kyrtsos, M. Matsubara and E. Bellotti, J. Phys.: Condens. Matter, 2020, 32, 365504.
77. S. Lany and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 235104.
78. A. Walsh, npj Comput. Mater., 2021, 7, 72.
79. Y. Chen, L. Wu, D. Liang, P. Lu, J. Wang, J. Chen, H. Cao and L. Han, Comput. Mater. Sci., 2021, 188, 110169.
80. L. Zhu, C. R. A. Catlow, Q. Hou, X. Zhang, J. Buckeridge and A. A. Sokol, J. Mater. Chem. A, 2023, 11, 15482–15498.
81. P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. R. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjøvoll, A. Fahmi, A. Schäfer and C. Lennartz, J. Mol. Struct.: THEOCHEM, 2003, 632, 1–28.
82. J. Kästner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander and P. Sherwood, J. Phys. Chem. A, 2009, 113, 11856–11865.
83. Y. Lu, M. R. Farrow, P. Fayon, A. J. Logsdail, A. A. Sokol, C. R. A. Catlow, P. Sherwood and T. W. Keal, J. Chem. Theory Comput., 2019, 15, 1317–1328.
84. E. Aprà, E. J. Bylaska, W. A. de Jong, N. Govind, K. Kowalski, T. P. Straatsma, M. Valiev, H. J. J. van Dam, Y. Alexeev, J. Anchell, V. Anisimov, F. W. Aquino, R. Atta-Fynn, J. Autschbach, N. P. Bauman, J. C. Becca, D. E. Bernholdt, K. Bhaskaran-Nair, S. Bogatko, P. Borowski, J. Boschen, J. Brabec, A. Bruner, E. Cauët, Y. Chen, G. N. Chuev, C. J. Cramer, J. Daily, M. J. O. Deegan, T. H. Dunning, M. Dupuis, K. G. Dyall, G. I. Fann, S. A. Fischer, A. Fonari, H. Früchtl, L. Gagliardi, J. Garza, N. Gawande, S. Ghosh, K. Glaesemann, A. W. Götz, J. Hammond, V. Helms, E. D. Hermes, K. Hirao, S. Hirata, M. Jacquelin, L. Jensen, B. G. Johnson, H. Jónsson, R. A. Kendall, M. Klemm, R. Kobayashi, V. Konkov, S. Krishnamoorthy, M. Krishnan, Z. Lin, R. D. Lins, R. J. Littlefield, A. J. Logsdail, K. Lopata, W. Ma, A. V. Marenich, J. Martin del Campo, D. Mejia-Rodriguez, J. E. Moore, J. M. Mullin, T. Nakajima, D. R. Nascimento, J. A. Nichols, P. J. Nichols, J. Nieplocha, A. Otero-de-la-Roza, B. Palmer, A. Panyala, T. Pirojsirikul, B. Peng, R. Peverati, J. Pittner, L. Pollack, R. M. Richard, P. Sadayappan, G. C. Schatz, W. A. Shelton, D. W. Silverstein, D. M. A. Smith, T. A. Soares, D. Song, M. Swart, H. L. Taylor, G. S. Thomas, V. Tipparaju, D. G. Truhlar, K. Tsemekhman, T. Van Voorhis, Á. Vázquez-Mayagoitia, P. Verma, O. Villa, A. Vishnu, K. D. Vogiatzis, D. Wang, J. H. Weare, M. J. Williamson, T. L. Windus, K. Woliński, A. T. Wong, Q. Wu, C. Yang, Q. Yu, M. Zacharias, Z. Zhang, Y. Zhao and R. J. Harrison, J. Chem. Phys., 2020, 152, 184102.
85. J. Buckeridge, C. R. A. Catlow, D. O. Scanlon, T. W. Keal, P. Sherwood, M. Miskufova, A. Walsh, S. M. Woodley and A. A. Sokol, Phys. Rev. Lett., 2015, 114, 16405.
86. J. Buckeridge, C. R. A. Catlow, M. R. Farrow, A. J. Logsdail, D. O. Scanlon, T. W. Keal, P. Sherwood, S. M. Woodley, A. A. Sokol and A. Walsh, Phys. Rev. Mater., 2018, 2, 54604.
87. Q. Hou, J. Buckeridge, A. Walsh, Z. Xie, Y. Lu, T. W. Keal, J. Guan, S. M. Woodley, C. R. A. Catlow and A. A. Sokol, Front. Chem., 2021, 9, 780935.
88. P. J. Wilson, T. J. Bradley and D. J. Tozer, J. Chem. Phys., 2001, 115, 9233–9242.
89. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982–9985.
90. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
91. Y. Zhao, B. J. Lynch and D. G. Truhlar, J. Phys. Chem. A, 2004, 108, 2715–2719.
92. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
93. B. G. Dick and A. W. Overhauser, Phys. Rev., 1958, 112, 90–103.
94. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
95. J. D. Gale, J. Chem. Soc., Faraday Trans., 1997, 93, 629–637.
96. A. Logsdail, Software to fit ECPs for the boundary on a QM/MM calculation with ChemShell, https://www.github.com/logsdail/fit_my_ecp, accessed February 22, 2022.
97. Z. Xie, Y. Sui, J. Buckeridge, C. R. A. Catlow, T. W. Keal, P. Sherwood, A. Walsh, D. O. Scanlon, S. M. Woodley and A. A. Sokol, Phys. Status Solidi A, 2017, 214, 1600445.
98. Z. Xie, Y. Sui, J. Buckeridge, C. R. A. Catlow, T. W. Keal, P. Sherwood, A. Walsh, M. R. Farrow, D. O. Scanlon, S. M. Woodley and A. A. Sokol, J. Phys. D: Appl. Phys., 2019, 52, 335104.
99. A. A. Sokol, S. A. French, S. T. Bromley, C. R. A. Catlow, H. J. J. van Dam and P. Sherwood, Faraday Discuss., 2007, 134, 267–282.
100. D. J. Wilson, A. A. Sokol, S. A. French and C. R. A. Catlow, MRS Proc., 2004, 848, 371–376.
101. C. R. A. Catlow, A. A. Sokol and A. Walsh, Chem. Commun., 2011, 47, 3386–3388.
102. A. J. Logsdail, C. A. Downing, T. W. Keal, P. Sherwood, A. A. Sokol and C. R. A. Catlow, Phys. Chem. Chem. Phys., 2016, 18, 28648–28660.
103. X. Zhang, L. Zhu, Q. Hou, J. Guan, Y. Lu, T. W. Keal, J. Buckeridge, C. R. A. Catlow and A. A. Sokol, Chem. Mater., 2023, 35, 207–227.
104. A. A. Sokol, S. T. Bromley, S. A. French, C. R. A. Catlow and P. Sherwood, Int. J. Quantum Chem., 2004, 99, 695–712.
105. F. A. Kröger and H. J. Vink, in Solid State Physics, ed. F. Seitz and D. Turnbull, Academic Press, 1956, vol. 3, pp. 307–435.
106. J. R. Rumble, CRC Handbook of Chemistry and Physics, CRC Press/Taylor & Francis, Boca Raton, FL, USA, 2020.
107. CRC Handbook, CRC Handbook of Chemistry and Physics, CRC Press/Taylor & Francis, Boca Raton, FL, USA, 2021.
108. W. Jost, J. Chem. Phys., 1933, 1, 466–475.
109. H. Jiang and Y.-C. Shen, J. Chem. Phys., 2013, 139, 164114.
110. X. Zhang, T. Liu, L. Zhu, J. Guan, Y. Lu, T. W. Keal, J. Buckeridge, R. Catlow and A. A. Sokol, Angew. Chem., Int. Ed., 2023, 62, e202308411.
111. Z. Xie, Y. Sui, J. Buckeridge, C. R. A. Catlow, T. W. Keal, P. Sherwood, A. Walsh, M. R. Farrow, D. O. Scanlon, S. M. Woodley and A. A. Sokol, J. Phys. D: Appl. Phys., 2019, 52, 335104.
112. C. G. Van de Walle and J. Neugebauer, Nature, 2003, 423, 626–628.
113. P. G. Moses, M. Miao, Q. Yan and C. G. Van de Walle, J. Chem. Phys., 2011, 134, 84703.
114. M. C. Benjamin, C. Wang, R. F. Davis and R. J. Nemanich, Appl. Phys. Lett., 1994, 64, 3288–3290.
115. C. I. Wu, A. Kahn, E. S. Hellman and D. N. E. Buchanan, Appl. Phys. Lett., 1998, 73, 1346–1348.
116. C. I. Wu and A. Kahn, Appl. Surf. Sci., 2000, 162–163, 250–255.
117. T. Kozawa, T. Mori, T. Ohwaki, Y. Taga and N. Sawaki, Jpn. J. Appl. Phys., 2000, 39, L772–L774.
118. S. P. Grabowski, M. Schneider, H. Nienhaus, W. Mönch, R. Dimitrov, O. Ambacher and M. Stutzmann, Appl. Phys. Lett., 2001, 78, 2503–2505.
119. P. Reddy, I. Bryan, Z. Bryan, J. Tweedie, S. Washiyama, R. Kirste, S. Mita, R. Collazo and Z. Sitar, Appl. Phys. Lett., 2015, 107, 91603.
120. M. S. Miao, A. Janotti and C. G. Van De Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 1–9.
121. P. Strak, P. Kempisty, K. Sakowski and S. Krukowski, J. Vac. Sci. Technol., A, 2017, 35, 21406.
122. A. Alkauskas, P. Broqvist and A. Pasquarello, Phys. Rev. Lett., 2008, 101, 046405.
123. J. L. Lyons and C. G. Van de Walle, npj Comput. Mater., 2017, 3, 12.
124. I. A. Aleksandrov, V. G. Mansurov, V. F. Plyusnin and K. S. Zhuravlev, Phys. Status Solidi C, 2015, 12, 353–356.
125. I. A. Aleksandrov, T. V. Malin, D. S. Milakhin, B. Y. Ber, D. Y. Kazantsev and K. S. Zhuravlev, Semicond. Sci. Technol., 2020, 35, 125006.
126. S. B. Orlinskii, P. G. Baranov, A. P. Bundakova, M. Bickermann and J. Schmidt, Phys. B, 2009, 404, 4873–4876.
127. J. Buckeridge, Comput. Phys. Commun., 2019, 244, 329–342.
128. G. A. Slack and T. F. McNelly, J. Cryst. Growth, 1977, 42, 560–563.
129. M. Bickermann, B. M. Epelbaum and A. Winnacker, J. Cryst. Growth, 2004, 269, 432–442.
130. R. Yu, G. Liu, G. Wang, C. Chen, M. Xu, H. Zhou, T. Wang, J. Yu, G. Zhao and L. Zhang, J. Mater. Chem. C, 2021, 9, 1852–1873.
131. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
132. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
133. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
134. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 65, 035406.
135. M. Chase, NIST-JANAF Thermochemical Tables, American Institute of Physics-1, 4th edn, 1998.
136. M. Ilegems and H. C. Montgomery, J. Phys. Chem. Solids, 1973, 34, 885–895.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ta04335a

---