Computational study of native defects and defect migration in wurtzite AlN: an atomistic approach †

We derive an empirical, lattice energy consistent interatomic force ﬁ eld model for wurtzite AlN to predict consistently a wide range of physical and defect properties. Using Mott – Littleton techniques, we calculate formation energies of vacancies and interstitials, which show good agreement with previous ab initio calculations at the edge of the band gap. A novel N 3 − interstitialcy con ﬁ guration is proposed to be of lower energy than the octahedral-channel-centred counterpart. With the assistance of the QM/MM method, our potential can predict a VBM level ( − 7.35 eV) comparable to previous experimental measurements. We further investigate the migration mechanisms and energy barriers of the main intrinsic defects. For the vacancy migration, the axial migration barrier is found to be lower than the basal migration barrier, in contrast to previous calculations. Two interstitialcy migration mechanisms for the interstitial defects are proposed, the “ knock-out ” mechanism for Al interstitial and the “ hand-over ” mechanism for N interstitialcy defects. The new force ﬁ eld model proposed here demonstrates that the empirical two-body interatomic potential is still e ﬀ ective for the study of defect properties, electronic states, and other extended systems of III/V semiconductors and further can be employed with QM/MM embedded techniques.


Introduction
Aluminium nitride (AlN), together with gallium nitride (GaN) and indium nitride (InN) and in solid solution with various compositions of group-III nitrides (AlGaN, etc.), is one of the key materials for solid-state light sources.Having had the least attention among the three nitrides, AlN has started to gain more interest owing to its application in deep ultraviolet (wavelength shorter than 300 nm) light-emitting diodes (LEDs), [1][2][3] due to its very wide band gap (6.2 eV, equivalent of 200 nm).The state-of-the-art AlN-based deep-UV lab LED was successfully made with the lowest possible wavelength of 210 nm (5.9 eV), 3 the potential of which has stimulated great interest in the material.Some additional advantages, such as good thermal conductivity (285 W m −1 K −1 ) and high thermal stability, [4][5][6][7] make AlN increasingly important for energy and electrical devices. 8Recently, single photon emission was observed in an AlN thin lm, 9 showing AlN to be amongst the most promising materials for applications in quantum technology. 10efects plays a key role in the optical properties of AlN.Early experimental work by Slack and McNelly showed that defects, which are mainly formed during crystal growth, show a strong correlation with the optical properties of AlN. 113][24][25][26][27][28][29][30][31] The source of the observed optical properties in AlN has now become clearer owing to these computational studies, as the accuracy of the calculations has improved, by the usage of higher level of exchange-correlation energy functionals from local density approximation (LDA) to hybrid functionals.
However, these ab initio techniques are not favourable for large systems of thousands of atoms, due to their computational expense and scaling performance.Although our future work will explore the use of such methods, particularly employing quantum mechanical/molecular mechanical (QM/MM) methods, less computationally demanding methods based on interatomic potential models are essential to explore much larger systems and a wider range of congurations.3][34][35] They allow us to establish the overall nature of the defect structure of the materials, with the subsequent (QM/MM) focusing on electronic and spectroscopic properties of key defects using our newly derived interatomic potential in the hybrid QM/MM approaches.
Moreover, a robust interatomic potential model is very benecial in investigating disordered systems and the solid solutions.
Since the rst reported interatomic potential model for wurtzite AlN, 36 there have been a number of evolutions of the models.The earliest pair-wise Coulomb and Buckingham potential implementations were based on the ionic model of AlN; 36,37 such models require only a small number of empirical parameters and allow the prediction of defect energetics and the physical properties including lattice constants and dielectric constants. 37Later studies shied focus to lattice dynamics, with more emphasis on the covalent components of the bonding of the material, with molecular mechanical rather than an ionic approach to the development of the potential model.These models, including bond order Tersoff potential models, [38][39][40][41][42][43] three-body potential models 44,45 and other many-body models, 46 have allowed accurate modelling of the physical properties of the material.However, the more complex interatomic potentials require more parameters; moreover, their use is difficult in modelling defect and surface properties of the material.More recently, a machine learning interatomic potential model was developed for silicon, reproducing wide range of DFT results including elastic properties, surface energies, defect formation energies and defect migration barriers. 47However, the data-driven model lacks proper description of long-range interactions and suffers from its intrinsic feature of high dimensionality in exchange of transferability (the so-called "curse of dimensionality"). 47,48A simple empirically tted interaction potential model still proves to be effective in investigating defect energetics, and transferable to a wide range of structures and physical properties.
In this work, we propose a new pairwise and a 3-body interatomic potential models for wurtzite AlN, from which we will present a comprehensive study of intrinsic defect structure of the material.The new interatomic potential can accurately predict a wide range of physical properties including lattice constants, elastic constants and dielectric constants.The interatomic potential models presented are designed to be lattice energy consistent for reproducing the energies of electronic states accurately.The Mott-Littleton method is implemented for the electron/hole ionization energy and intrinsic defects in wurtzite AlN, focusing on the formation of vacancies and interstitials, which have been previously studied using DFT techniques.While maintaining high calculation efficiency for systems of thousands of atoms, our force eld is the rst AlN interatomic potential model that can be used to calculate defect formation energies accurately.A novel conguration of a highly charged nitrogen "interstitialcy" defect is found to be more stable than the channel-centred counterpart.With further investigation into the defect migration, our results give new insights into the vacancy migration barriers and the interstitialcy migration mechanisms of Al and N interstitial defects.Our results are of relevance to all III-V nitride systems and indeed to all the wurtzite structured materials.

Interatomic potential
The rst interatomic potentials for AlN were reported by Cormack, 36 using a shell model with formal charges.Later models include those of Chisholm 37 which employed partial charges.In both studies, only a Buckingham potential was used whereas in our work, we employ a more complex multi-functional model, giving a more accurate and smooth potential prole, which has been successfully applied to oxides lately. 33The parameterisations of the potential models in this work were completed using the General Utility Lattice Program (GULP) code. 49he interatomic potential model in this work is based on the Born model 50 of the ionic solid, as both the wurtzite structure adopted by AlN, and the dielectric properties of the material are consistent with the use of the ionic model.To obtain a precise and smooth prole of interatomic potentials across short-range and long-range distances, we tted a combination of pairwise long-range Coulomb, and short-range Buckingham, Lennard-Jones, and polynomial potentials to experimental data including lattice constants, elastic constants, dielectric constants, piezoelectric constants and phonon frequencies.
Regarding the short-range potentials, for the interaction between the cation and anion, here, we employ the Buckingham potential with respect to internuclear distance r: where A, r and C are potential parameters.Separate Buckingham potentials are applied to the 1st and the 2nd atomic neighbouring shells, so that the tting can be controlled separately.The polynomials are applied to interpolate between different Buckingham potentials and extrapolate to zero aer the 2nd neighbouring shell: where n is the index of summation and C n are polynomial coefficients.
A repulsive r −12 potential (implemented as a Lennard-Jones type without the attractive term) is applied covering the whole interatomic short range between cation and anion, to prevent atoms from accessing unphysically short separations: with C 12 also used as a parameter.Between anions, we employ the 4-range Buckingham ("Buckingham-4") potential of separated Born-Meyer (for 1st neighbouring anions between r min and r cut1 ) and dispersion (for 2nd neighbouring anions between r cut3 and r max ) terms, where concept of separately treating the 1st and the 2nd neighbouring anions adopted here is achieved in a condense form.A, r, C 12 , and the polynomial a n , b n , C n constants are all variable parameters that are adjusted during the tting process.
The chosen values for these parameters is discussed below.
In this model, we rst apply the pair potential approximation, which is commonly used for ionic solids.In addition, polarisation of the nitride ions is included via the shell model, 51 which provides a simple mechanical model of polarizability in which the ion is divided into a core linked by an harmonic spring to a shell, and which has proved effective in modelling dielectric properties of ionic materials.In the harmonic approximation, the polarizability of an ion in vacuo is given by: where Y is the charge of the shell and k is the spring constant.
The energy from such spring interaction can be obtained using a harmonic square term which is complemented here by an anharmonic quartic potential that in turn will reduce the polarisability and dampen the shell displacement in the external eld: As noted, in our model, the aluminium ion is treated as unpolarizable, i.e., core only (charge +3), whereas the nitride ion is polarisable with both a core (charge +1.4) and a shell (charge (Y) −4.4).Two main short-range interactions are considered, i.e., "Al-core-N-shell" and "N-shell-N-shell".Between "N core" and "N shell", the interactions are controlled solely by harmonic spring (including square and quartic terms).
For many ionic materials, the Buckingham potential still proves to be an effective model.The parameters of the potential for "Al-core-N-shell" interaction are all tted empirically against experimental properties.As the Al ion is considered here to have a negligible polarizability, the dispersion term of "Al-core-N-shell" interaction is omitted (C = 0).Parameters of the repulsive "N-shell-N-shell" interaction, A and r in the Buckingham potential, are transferred from a widely used "O-shell-O-shell" potential, 52 which successfully predicts defect properties and defect diffusion energies for a range of oxide materials. 53The transferability of the potential can be justied by the isoelectronic Ne 10 structure of nitride and oxide ions and the minor role of this potential in our simulations.In particular, this potential affects the calculated defect properties only slightly and serves mostly in preventing unphysical congurations with short anion-anion separation distances.For the attractive dispersion parameter C in the "N-shell-N-shell" Buckingham potential, we applied the theory rst introduced in the work by Pyper 54 and made our approximation previously applied for GaN. 55According to the Slater-Kirkwood formula, 56 the dispersion parameter for anion-anion interaction can be determined from the frozen in-crystal polarizability, which is calculated by taking the spring constants tted to reproduce the high-frequency dielectric constants.
The coefficients of the polynomial functions are calculated to connect all the regions in the Buckingham potentials, under the condition that the 1st and 2nd derivatives at both end points must be matched.The parameters of zero-degree polynomial functions control the vertical offsets of the Buckingham potential, which is adjusted according to the lattice energy.These zero-gradient offsets do not affect the physical properties but can tune the lattice energy, a key term for obtaining our electron affinities and in turn calculating the defect formation energies.The cut-off distances of the short-range potentials are 10 Å, which extends potentials beyond nearest and next-nearest neighbours where the interactions are still signicant but allows for a smooth and rapid transition to zero thus reducing the computational costs.The tted parameters are reported in Table 1.
The resulting structural parameters, physical properties and lattice energy are obtained using the GULP code and shown in Table 2: the lattice constants and dielectric constants agree well with experimental and other ab initio computational results, although there are signicant discrepancies in the off-diagonal elastic constants.9][40][41][42][43][44][45][46] All these models, however, employ partial ionic charges, which can be problematic in defect calculations.We also note that although the previous two-body interatomic potential model presented a slightly better tting of elastic constants, the lattice energy was particularly low, 37 as a result of the choice of partial charge leading to a much lower Madelung potential.A reliable lattice energy value is crucial in our calculation of defect energy as presented in the following sections.
To demonstrate that the 3-body interatomic potential improves substantially the prediction of elastic constants while still using formal charges, we t another set of parameters by adding a 3-body potential to the same 2-body ones.Here the 3body interatomic potential takes the form of: where K and r 3B are the parameters tted empirically, R 3B 0 is the average Al-N bond length of perfect lattice, and r 12 and r 13 are the interatomic distances between Al and different neighbouring N atoms.Here we assume that the "Al core" interacts equally with the two neighbouring "N shells" in both axial and basal directions, so the two values of r are xed to be the same.On introducing these additional terms in the potential, the parameters of our 2-body interatomic potential for "Al core" and "N shell" interaction must be retted.The 2-body Buckingham potential of "N shell" and "N shell" remains unchanged.Table 3 presents all the parameters of the potential, and Table 2 shows the prediction of the properties using the new 3-body potential model.Now all the properties agree well with experimental and ab initio results, except for one of the piezoelectric constants.However, as we will discuss in the following section, the adopted 3-body model raises problems for defect calculations without some additional compensating terms in the force eld.Such a compensation can in turn reduce the accuracy of the predicted physical properties.Therefore, for our defect study, the 2-body interatomic potential is primarily applied.
Table 1 The 2-body short-range interatomic potential parameters of wurtzite AlN

Mott-Littleton and QM/MM calculations
Our potential is implemented in Mott-Littleton 68 calculations of defect energies.The method which is discussed in detail in ref.
49 divides the lattice around the central defect into 3 regions (Fig. 1).In the inner region (region I), displacements of the atoms are calculated by explicit energy minimisation.Secondly, the displacement of atoms in the middle region (region IIa) is calculated within the harmonic approximation for interactions with respect to the central defect, with explicit calculation of the interaction energy with region I. Response by the ions of the outer region (IIb) to the Coulomb eld of the central defect is treated using a dielectric continuum model but the lattice summation performed explicitly.In this work, the radii of region I and region IIa are 21 Å (around 3600 atoms) and 36 Å (around 15 000 atoms), respectively; with these values all defect energies have converged to within 0.1 eV or better with respect to subsequent expansion.All the calculations were again undertaken using the GULP code. 69n this study, we also make use of calculated ionisation energies obtained using the hybrid QM/MM embedded cluster method as implemented in the ChemShell soware. 70,71The QM/MM method employed here contains three levels of theory.In the central cluster region where the electron ionization occurs, the DFT method is applied with the hybrid functional PBE0 72,73 and Def2-TZVP 74 basis sets.The QM cluster region is surrounded by a large MM environment controlled by the same 2-body interatomic potential.Between the QM and MM region, there is a transition layer balancing the different forces from both sides by the tted electron localising pseudopotentials (PP). 75Compared to the DFT method using the periodic boundary condition, the QM/MM method can calculate the ionization potential as it has unambiguous denition of the   vacuum level.This theory is explained in more detail in other work. 70,71,76,77Further detail of the settings for QM/MM calculations can be found in the ESI.†

Results
In this section, we apply our interatomic potential model and the Mott-Littleton method to calculate defect energies and the energies of electron and hole formation.We rst compare the defect energies calculated by the 2-body and the 3-body models.
We then obtain estimates of electron and hole formation energies, which facilitate the calculation of formation energies for both vacancy and interstitial species.Lastly, we model defect migration processes.

Defect energies from Mott-Littleton calculations
The intrinsic point defect energies of AlN are presented in Table 4.As we use the Mott-Littleton approach, all point defect energies are calculated with their formal charges.We note that the energies are given with respect to the perfect lattice and ions at innity.For both defect species and charge carriers we use the Kröger-Vink notation.
There are two distinct high-symmetry sites, or voids, for interstitials in the wurtzite structure adopted by AlN: the tetrahedral (TE) and octahedral (OC) sites (Fig. 2).An interstitial at the TE site has four cations and four anions as the nearest neighbours, while that at the OC site has six nearest neighbours (three cations and three anions).We may also consider the OC interstitial sites from the perspective of the z-axis, where the interstitial at the OC site is situated at the centre of the "channel".Previous calculation found the Al interstitial at the OC site is more stable than at the TE site in AlN, 26 with analogous results found for GaN. 28,78In our work, when the Al 3+ interstitial is initially located at a TE site it moves during the relaxation towards the OC site.Hence, only the results for this conguration of Al 3+ interstitial defect are presented here.
For the N 3+ interstitial, two distinct stable congurations are found: the octahedral-central (OC-central) conguration and the "interstitialcy" conguration.Our calculations show the defect energy of OC-central type N interstitial is close to its cation counterpart with less than 1% difference in energy (Table 4).Such a conguration is consistent with the previous ab initio simulations of GaN 79,80 and InN 81the only reports obtaining the stable OC-central N interstitial in charge state −3 in wurtzite structured nitride systems.
The N 3− defect is found to have another stable "interstitialcy" conguration in the system.The new conguration is obtained by rst relaxing the OC-central conguration defect, and then adjusting the positions of the defect and its surrounding neighbouring N 3− ions.Aer complete relaxation, the introduced N 3− ion (N1 in Fig. 3(a)) is displaced from the ideal OC site (by 0.34 Å) in the lattice basal plane towards one of the three nearest neighbouring N 3− ions.The on-site N 3− ion (N2 in Fig. 3(a)) is then nudged from its original position (by 0.62 Å) towards the axial direction along the line interconnecting the defect and the on-site ion (Fig. 3(b) and (c)).The two N ions comprise a special "interstitialcy" or split-interstitial conguration, in which two N ions straddle the same N lattice site at the distance of 2.40 Å, forming an interstitial-vacancyinterstitial defect complex structure.We calculated the defect energy of such conguration to be 0.24 eV lower than the OCcentral type (Table 4), suggesting that the OC-central N interstitial is a local energy minimum conguration in the system.
The latter type of N interstitial defect has been widely overlooked by other studies of III-V nitride system.The term "splitinterstitial" was used in another conguration reported in many previous ab initio studies, 26,31,65,82 where their N split-interstitial defects have a form of N 2 molecule positioned at the N vacancy site.The bonding in those N-N structure shows a strong covalency nature, which were found to be stable in the charge state from +3 to −1 due to the electronic occupancy at their antibonding-pp-like orbitals. 82We note that our novel ionic N 3− interstitialcy defect is a separate conguration from those covalent split-interstitial discussed previously.In the case of our N interstitialcy defect in charge state −3, the two N 3− ions are rearranged at the lowest energy conguration as a result of the Madelung potential (electrostatic interactions) and the short-range interatomic potential, which is purely an ionic response.The defect structure of anion interstitial in formal charge is not new in other material system (see examples in alkali halides 83 and oxides 84 ).In particular, our N 3− interstitialcy  defect shows the closest resemblance to the O 2− interstitial defect in UO 2 , 52 where the O 2− interstitial distorted from bodycentre position by 0.2 Å and 0.27 eV lower in energy compared to the symmetric counterpart.Similar electrostatically-driven reconstruction of interstitial point defect has recently been discussed in another report. 85he authors stated and proved that high-symmetry relaxation methodology adopted by many computational defect studies could easily lead to metastable state (local minima) rather than the global minimum, 86 which can be found using their novel symmetry-breaking technique with standard DFT calculations.They found the defect energy of Te 2− split-interstitial is 0.2 eV lower than the higher symmetric tetrahedral Te 2− interstitial in CdTe, which is consistent with our lower symmetric N 3− interstitialcy defect conguration.Different congurations of interstitialcy defects, including OC-OC, TE-TE, and OC-TE types, have also been explored in our system, but all of them are unstable in our simulations.Our N 3− interstitialcy conguration has also been conrmed by the QM/MM calculation, which will be presented in full detail in our subsequent paper.
The same defect calculation proves to be problematic for our 3-body potential model.We found that the Mott-Littleton defect calculation cannot converge for all four intrinsic defects without an additional correction in the force eld.The preexponential parameter K is tted empirically to be negative (Table 3), which indicates the 3-body potential acts as an attractive force.As a defect is introduced, the cluster is likely to contract due to the additional attractive force, which leads to abnormally high potential gradients at energy minimum position.As a result, a repulsive term must be added in the force eld to counteract the force for convergence (see ESI † for more detail).However, the compensation term results in the physical properties deteriorating.Therefore, the 2-body potential is chosen to be more effective for the defect calculation, as the defect energies are close to the ones from our 3-body potential (Table 4).All subsequent calculations use solely the 2-body potential.

Electron and hole formation
Computational models based on interatomic potentials cannot provide accurate estimates of electronic energies, but a useful estimate of hole and electron energies is still possible.We assume that the conduction band minimum (CBM) can be represented by an electron localising on a cation (Al) and the valence band maximum (VBM) by a hole on an anion (N), i.e., the electron state is thus approximated as a localised Al 2+ state, and the hole state as a localised N 2− state.(Material's electron affinity and ionisation potential determine the absolute position of CBM and VBM below vacuum level set at zero.)By using our interatomic potential model and Mott-Littleton method in GULP, defect energies of Al 2+ and N 2− are calculated.To model optical excitations, only shells are relaxed during these calculations, i.e., there is no relaxation of the cores.In addition, we must include appropriate intra-atomic terms.Here, we make another simplication that such intra-atomic energies are equal to their respective ionisation energies as free ions.The third ionisation potential of Al is available from experiment.The third electron affinity of N is not a measurable quantity but has been calculated using the QM/MM techniques. 70All the energy terms are collected in Table 5.
The calculated VBM is 0.85 eV higher than that based on the experimental ionization potential on the AlN surface of 8.2 eV. 88e obtain a higher band gap energy (10.60 eV) than reported experimentally (6.23 eV at low temperatures 87 ).The calculated large positive value for the CBM is clearly inaccurate, which can be attributed to our assumption of a localised Al 2+ state representing a CBM state.The experimental value of the electron affinity of AlN is uncertain with estimates, ranging from 0.25-1.9eV, 88-90 but our calculated CBM value is far from these values.As our approach appears to be more reliable for estimating a hole state, we can estimate the revised CBM value by adding the experimental band gap energy at low temperatures to our VBM value.The revised CBM now falls within the range of experimental values.

Defect reaction energies
The reaction forming an Al vacancy ðV 00 0 Al Þ in the AlN lattice is represented by the following equations depending on the Al and N chemical potentials: N-rich/Al-poor: N-poor/Al-rich: The "N-rich/Al-poor" or "N-poor/Al-rich" conditions indicate the respective reactions occurring in the case of either "N 2(g) excess/Al (s) deciency" or "N 2(g) deciency/Al (s) excess".Considering thermodynamic equilibrium in the Al-N binary systems, the formation of solid AlN imposes a bound on their respective chemical potentials: N 2ðgÞ /AlN ðsÞ whose energy is DH AlN(s) with the experimental value of −3.296 eV. 91he neutral aluminium vacancy created will initially accept three electrons resulting in forming three holes in the system.The energy of forming electrons and holes in the following analysis is determined by the electron chemical potential, or Fermi level, E F , whose values in the material would typically lie in the band gap, between VBM and CBM, −E h # E F # −E e (see Table 5), although the heavily doped system can also fall outside this range.

Al is the defect energy of Al vacancy in
Applying the same approach to the N vacancy, Al interstitial, and N interstitial, the following defect reactions can be formulated: N vacancy ðV N Þ: (N-rich/Al-poor): (N-poor/Al-rich): Al interstitial ðAl i Þ: (N-rich/Al-poor): (N-poor/Al-rich): N interstialcy ðN 00 0 int-cy Þ: (N-rich/Al-poor): (N-poor/Al-rich): And the respective formation energies are as following: where A 1−3 N is the sum of the three electron affinities of N atom.The experimental electron affinity energies of nitrogen atom are not readily available and would, in fact, depend on the stabilising environment of the triply negatively charged species, but we can estimate the sum of three electron affinities via the Born-Haber cycle, for AlN formation, using our calculated 2body lattice energy in Table 2. Our resulting estimate of the sum of the energies of the three electron affinities of N atom is 27.04 eV.The experimental value of the dissociation energy of N 2 , D N 2 , used for these calculations is 9.79 eV. 91able 6 summarises our calculated formation energies.For vacancy defects, formation energy V N is higher in N-rich condition, while the formation energy of V 00 0 Al is higher in Alrich condition, which is as expected and agrees with other studies using rst-principles methods. 23,27,30,65,92,93The comparison with earlier electronic structure calculations is not straightforward as the latter are typically reported as a function of the position of the Fermi level and only in a graphical form, so, in Table 6, we compare reported calculated formation energies for Fermi levels at the band edges.The energies for the N vacancy and Al vacancy in either environment condition agree well with previous DFT calculations at the same Fermi level positions.It was predicted in the previous ab initio calculations that the V N is the most stable defect close to the VBM, 26,30,65 which indicates its strong correlation to donor species in AlN.For interstitial defects, the formation energy of Al i in Al-rich condition is, however, not as close to previous DFT results.The greater value in N-rich condition can be attributed to a lower enthalpy of AlN formation used in the other work. 94Further, there is no available data for N 00 0 int-cy in the literature, as the formally charged N interstitial defect is reported to be not as stable as lower charged states of the N interstitial at the CBM.
The results demonstrate our interatomic potential model can obtain accurate defect formation energies, which reproduce well those are calculated using rst-principles methods.Note our model requires a much lighter computational load (more detail about the benchmarking is presented in ESI †).

Formation of Frenkel and Schottky defect pairs
From the energies reported in (Table 4) we can calculate energies of the Schottky and Frenkel pairs listed below: Schottky defect: Anti-Schottky defect: Anion Frenkel defect: Cation Frenkel defect: All the energies are summarised in Table 7.
For AlN, Schottky disorder has the lowest energy, suggesting that Al and N vacancies will be the dominant disorder type in the material.Cormack 36 and Chisholm 37 arrived at the same conclusion using interatomic potential methods, but they did not report results for the "Anti-Schottky" energy.The considerably higher Schottky defect energy here compared with that reported by Chisholm et al. reects the choice of interatomic potential, especially the differences in ionic charges.The higher anti-Schottky and N Frenkel defect energies reported by Cormack can be attributed to the different N interstitial structure, as discussed above.So far, by accurately predicting the physical properties, the electron/hole states, the defect energetics/structures, and the defect formation, we have tested the robustness of our interatomic potential model.The potential model can provide a reliable embedded environment for our QM/MM system, where a full scope of AlN point defect investigation will be carried out.We should note that our interatomic potential model has not been designed to simulate point defects in other charge states or other covalent defects, which may be a drawback for particular applications, and caution should be applied when modelling, e.g., N diffusion using molecular dynamics.

Point defect migration
The motion of native defects and surrounding ions controls atomic transport and is fundamental to the understanding of material growth and degradation.However, there is very little information on such processes in AlN, despite extensive interest in III-V nitride semiconductor materials over the last few decades.Experimental studies have been limited to oxygen 95 or dopant diffusion, 96,97 with none on intrinsic self-diffusion.][100][101][102] But there has not been a detailed computational investigation on intrinsic defect migration in AlN, only some comparisons with GaN. 31,101Here we will present our computational results of a systematic study of defect migration in AlN and compare with currently available data from studies on GaN and InN.
Our approach is to perform a comprehensive search of the energy landscape for the migrating species.Along each path, 10-20 points with equal distance are marked and rational function optimization calculations are conducted at each point until a suitable transition state is identied.The transition states are further conrmed by identifying a single imaginary vibrational mode frequency related to the defect.A nal check involves a steepest descent optimisation with a small step size from each transitional state in both directions (to initial and nal states).Hence, we dene a valid migration path as that involving a continuous route which can be successfully connected by steepest descent minimisation via one or multiple transitional states.We note the need for the Mott-Littleton defect cluster centre, to avoid bias, to be at the mid-point of the path so that the activation energies are accurate and reliable.
3.5.1.Vacancy migration.The migration of a vacancy involves a neighbouring ion moving into the vacancy.We investigate the migration of V 00 0

Al and V
N in the wurtzite AlN lattice in the axial direction (out-of-plane) and basal direction (in-plane, perpendicular to c axis, Fig. 4).In Table 8 we present our calculated activation energies and other currently available data from the literature on AlN, GaN, and InN for comparison.
Our results show an energy barrier of 1.67 eV with Al migrating in the axial direction and 2.62 eV in the basal direction (Table 4).The initial ground state positions, transition state positions and the interatomic distances are presented in Fig. 5. Previous calculations of the vacancy defect migration present conicting conclusions.Our higher basal barrier of Al vacancy migration matches the results from the work of Aleksandrov et al., 101 where by using a rst-principles method the barrier for the Al vacancy migration was found to be higher than that for the Ga vacancy migration in GaN.However, Warnick et al. 99 found a lower barrier (1.6 eV) for Al vacancy migration in AlGaN alloy than in pure GaN, which agrees with our lower barrier results.We found a similar mechanism for N vacancy migration.Our calculated N vacancy migration barrier is 2.92 eV and 2.20 eV in the basal and axial directions respectively.Even though there are no available data for AlN with which to compare, our results are higher than the calculated activation energies in GaN, which agrees with the reported higher barrier trend of N vacancy in the AlGaN alloy 99 and the previous DFT calculation results. 31ompared with the ionic radii of Al 3+ and Ga 3+ ions, that of the N 3− ion is much bigger (by a factor of three compared with Al 3+ and of two for Ga 3+ ), so it should be expected that the migration  will be more hindered in AlN, with higher migration barriers.
Previous studies of vacancy self-diffusion in similar materials were reported to be isotropic. 78,99,100A more recent DFT defect calculation in InN shows an anisotropic migration behaviour, but the trend is inverse to ours with higher barrier in the axial direction. 103To our knowledge, no experimental data are available for direct comparison to date.
3.5.2.Interstitial migration.We rst discuss Al interstitial migration, and then that of the N interstitial; different approaches are needed for the two species, as explained below.
As stated in Section 3.3, Al interstitials are stable at the OC site.Therefore, the migration path between two nearest neighbour OC sites can either be along the hexagonal channel (along the c-axis) or horizontally passing the TE site.Here we apply the  same methodology as for the vacancies where the steepest descent search follows valid saddle point identication.The latter motion involves an "interstitialcy" migration mechanism in which the migrating Al interstitial knocks off and replaces a host Al ion, with the "knocked-out" Al ion becoming a new interstitial defect moving towards its closest neighbouring OC site (Fig. 6).This horizontal migration is energetically more favourable than the vertical migration through the hexagonal channel (Table 9).Our activation energy is very close to that of Ga interstitial migration in GaN, which again supports our reasoning in Section 3.5.1.
We nd that migration of the N interstitial occurs via a more complex mechanism.The N interstitial is easily pulled or dragged by other on-site N ions resulting in the distortion around the defect.Like the investigation of the Al interstitial migration, we assume that there could also be parallel and perpendicular motion for N interstitial migration with respect the c axis, so saddle point and initial/nal stable states searches are also undertaken here for both.A more detailed account of the procedure is given below.
For the migration through the hexagonal channel, at the saddle point three nearest neighbouring N ions are pushed away by the defect with different bond lengths-all three are shorter than the stable interstitialcy defect (Fig. 7).This interaction between three N ions makes it more difficult for the defect to return to the initial/nal stable conguration, inducing a higher migration barrier.
For the other migration mechanism, we found a path in a single conventional unit cell with a much lower barrier.Here we applied an iterative approach to nd the transition state.At rst, we located two close stable interstitialcy defect congurations at neighbouring lattice layers.These two congurations are rotationally symmetric around the 3-fold wurtzite c axis running through a hexagonal channel, so they are likely to be the initial and nal states of the path.Initially, the interstitial N defect (N1 in Fig. 8) was moved along the straight line from N2 to N3. N1 was kept xed on the path, as we relaxed only the surrounding ions to give us potential energy surface with respect to N1. Aerwards, we conducted Rational Functional Optimisation (RFO) transition state searches on each resulting conguration until a transition state (proved by the identication of a single imaginary vibrational mode frequency) was found, resulting in a new possible migration path of two straight lines connecting three points (N1 to N2 and N1 to N3).Subsequently, we again moved the defect along this new path and repeated the same process iteratively until the highest saddle point was found.Finally, as with all other types of migration, a steepest descent optimisation was conducted to verify the path.
The resultant migration motion is illustrated in Fig. 8 and is close to previous reports on GaN. 78The path has two rotationally symmetric interstitialcy defects as its initial and nal stable congurations.At the saddle point, another nearest neighbouring N ion (N3) is pushed away from its original position while the initial interstitialcy conguration is maintained with a stretched bond.As the defect approaches its nal position forming a new stable pair, the "le-over" N (N2) ion slowly moves back to the closest lattice site.This "hand-over" motion between three N ions is signicantly different from the "knock-on" mechanism of the Al interstitial, where the former involves three N ions and the latter one only involves two Al ions.
The migration barriers/activation energies for both mechanisms are presented in Table 9.Our result shows the interstitialcy migration mechanism is energetically more favourable than the hexagonal channel migration.This barrier is also comparable to the result of the same migration mechanism in GaN, which supports the earlier reasoning related to the N ion size.
A comprehensive atomistical investigation of point defect migration suggests that the energy barriers and migration paths of self-diffusion behaviour in AlN are overall close to their counterparts in GaN.The migration barrier of the Al vacancy is

Summary and conclusions
We have presented an interatomic potential model for AlN, which we have used to calculate both lattice and defect properties.Physical properties, including lattice constants, and dielectric constants can be reproduced and agree well with experimental data and other theoretical results; the results for the elastic constants and bulk moduli show poorer agreement.We improve the elastic properties by utilizing an additional 3body interatomic potential, but it performs poorly for the defect calculation without the additional correction, which in turn worsens our physical property prediction.The classical Mott-Littleton method has been employed to investigate formation energies of formally charged native vacancy defects, interstitial defects, and Frenkel and Schottky defect pairs.We gain generally good agreement with previous classical and electronic structure calculations, although in some cases there are signicant differences.We propose a new ground state N 3− interstitialcy conguration, in which an on-site N atom is nudged away from its original position as a response to the Madelung potential and the short-range force from the N 3− defect.A comprehensive investigation of point defect migration suggests that the energy barriers of self-diffusion in AlN are very close to the counterparts in GaN.We found lower migration barriers for vacancy migration in the axial direction than in basal direction.This result contrasts with previous DFT calculations.We propose two interstitialcy migration mechanisms for the interstitial defects: the "knock-out" mechanism for Al interstitial and the "hand-over" mechanism for N interstitialcy defects, both of which are more favourable than c-axis channel migration.The detailed mechanisms of defect migration uncovered in this work would play a signicant role across the whole class of wurtzite structured material.Furthermore, with the assistance of the QM/MM method, the model can predict the VBM (−7.35 eV) accurately compared with experimental data, which shows the approach here still has use for electronic structure calculations.Future work will employ the hybrid QM/ MM method to investigate other problems including defect formation energy with different electronic charge states and their donor/acceptor behaviour.funded by EPSRC (EP/P020194 and EP/T022213).The co-author Qing Hou acknowledges funding from Shanghai Pujiang Program (grant 22PJ1411300).

Fig. 2
Fig. 2 Tetrahedral (TE) and octahedral (OC) sites for interstitial defect in wurtzite AlN.The z axis is oriented along the c crystallographic direction.

Fig. 3
Fig. 3 (a) The N 3− interstitialcy defect after relaxation (N1 is the introduced interstitial defect, N2 is the host atom) and the comparison to the local configuration of N 3− OC-central interstitial (in transparent) along the (b) x-axis and (c) z-axis.

Fig. 4
Fig. 4 Axial and basal directions of vacancy migration in wurtzite AlN (note that these arrows are not actual paths, merely an illustration of migration directions).

Fig. 7
Fig. 7 Atomic visualisation of the transition state of N Interstitial defect migration through hexagonal channel: see from x-axis (left); see from zaxis (right).

1. 67
eV, at least 0.23 eV lower than the calculated Ga vacancy migration barrier in GaN.Contrary to previous calculations, the vacancy migration is not found to be isotropic in our analysis.Different heights of migration barriers are assigned to the axial and the basal directions.For the interstitial defect migration, the energy barrier of c-channel Al interstitial migration is 1.42 eV, which is 1.58 eV lower than that for the GaN counterpart.Different "interstitialcy" migration mechanisms are seen in cation and anion interstitial migration.For the Al interstitial: the on-site Al atom can be "knocked-out" and replaced by the approaching interstitial defect.And for the N interstitial, to maintain the stable interstitialcy conguration, the interstitial rst bonded with one on-site N atom is handed over to the other closest on-site N atom.The study of point defect migration here provides further conrmation of the integrity of our new interatomic potential model.

Fig. 8
Fig. 8 Atomic visualisation of the "hand-over" mechanism of N Interstitial defect migration (looking from x-axis): (a) initial state; (b) transitional state (saddle point); (c) final state.

Table 4
Intrinsic defect energies in eV of wurtzite AlN obtained using the method of interatomic potentials

Table 4 ,
DH Sub(Al) heat of sublimation of Al solid (3.43 eV (ref.91)) and I 1−3 Al the ionisation energy of Al atom to Al 3+ ion.The experimental ionisation energies of Al can be extracted from the database

Table 5
Electron and hole formation energy in wurtzite AlN

Table 7
Frenkel and Schottky defect formation energies calculated here and compared with results from previous studies using interatomic potential methods

Table 8
Activation energies (eV) of vacancy migrations in AlN compared to the literature data

Table 9
78,100tion energies (eV) of interstitial point defect migration in AlN compared to literature data for GaN This work GaN78,100