# Fileset

[APL26-AR-ICNS2026-02280.pdf](https://mdr.nims.go.jp/filesets/7b20c7c4-7c40-44b0-abcb-3d9a82a809de/download)

## Creator

[Hiroyuki Ishihara](https://orcid.org/0009-0002-5824-430X), [Masato Oda](https://orcid.org/0000-0002-8731-5122), [Tsuyoshi Miyazaki](https://orcid.org/0000-0003-3534-4404)

## Rights

This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in Hiroyuki Ishihara, Masato Oda, Tsuyoshi Miyazaki; Stable structures of (VGaVN)n defects in bulk GaN. Appl. Phys. Lett. 24 August 2026; 129 (8): 081908 and may be found at https://doi.org/10.1063/5.0331898.[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Stable structures of (VGaVN)                    <i>n</i>                    defects in bulk GaN](https://mdr.nims.go.jp/datasets/6a03e2bc-6b1f-40f4-bada-1b5a22ac8a12)

## Fulltext

Stable structures of (VGaVN)n defects in bulk GaN  Stable Structures of (VGaVN)𝑛 defects in bulk GaN  Hiroyuki Ishihara,1,2 Masato Oda,1, 2, a) and Tsuyoshi Miyazaki2  1Department of Materials Engineering, Wakayama University, Wakayama, 640-8510, Japan 2Research Center for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0044, Japan     Vacancy-type defects in GaN can migrate and aggregate during post-growth annealing. Here, we perform large-scale first-principles calculations to clarify the stability of (VGaVN)𝑛 aggregates (𝑛 = 1-6) in bulk GaN. For several structure models at each 𝑛, the defect formation energy of (VGaVN)𝑛  and the energy gain for the aggregation of (VGaVN)𝑛−1   and VGaVN  are calculated to investigate the stability during the stepwise growth. The calculated results show that the aggregation is always favorable and the energy gain for this process is significantly large at 𝑛 = 3 and 6. This result is consistent with the recently reported positron-annihilation experiments. The structural stability and the electronic structure of the various models are discussed with respect to the number of dangling bonds and lattice relaxations, depending on the morphology of (VGaVN)𝑛. ____________________________ a) Electronic mail:  moda@wakayama-u.ac.jp.             Gallium nitride (GaN) is a key wide-bandgap semiconductor for next-generation power and high-frequency electronics[1–6]. Despite rapid progress in device development, device performance and long-term reliability remain strongly influenced by defects. In particular, vacancy-type defects can act as carrier traps and efficient recombination centers. Such defects can be introduced during growth and subsequent processing (e.g., ion implantation) and can evolve during post-growth annealing for dopant activation and damage recovery. In recent device architectures, controlling defect evolution during thermal processing is essential for stable operation and reproducible device characteristics[7–9]. Upon annealing, vacancies can migrate, interact, and form larger complexes that are often more electrically detrimental than isolated point defects[10–12]. Thus, annealing steps for dopant activation and damage recovery can promote defect aggregation and the emergence of complex defects.  It has been suggested from experiments that vacancy aggregation in GaN can occur at technologically relevant annealing conditions. Positron annihilation spectroscopy (PAS) has provided direct evidence for vacancy aggregation in Mg-implanted GaN[13, 14]. In particular, vacancy aggregation starts already at moderate annealing temperatures (∼800 °C), leading to vacancy clusters such as (VGaVN)3, and after annealing above 1000 °C the measured (S,W) parameters suggest that (VGaVN)3 clusters become the major defects[14-16]. These PAS results indicate that VGaVN pairs can act as building blocks for thermally This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 2  activated aggregation and that the trimer complexes exist under realistic annealing conditions. Consistent with the PAS evidence for vacancy aggregation, divacancy-related complexes have also been discussed as plausible electrically active centers in GaN, which motivates an atomistic understanding of the structures and electronic states of (VGaVN)𝑛 aggregates. In addition, GaN sometimes exhibits larger void-like defects with pronounced anisotropy, such as hexagonal-based pyramidal cavities bounded by specific facets[17]. Such pyramidal and cavity-type defects have been reported particularly in heavily Mg-doped GaN and p-type GaN:Mg materials. In this respect, it is important to clarify how the stepwise aggregation of  (VGaVN)𝑛 can be realized at atomic scale.   Despite clear experimental evidence for vacancy aggregation, the atomic-scale structures and detailed growth pathways remain unclear. Previous first-principles studies have treated the GaN divacancy: Li et al. found short-distance configurations to be energetically favored [18], while Gohda discussed stable local configurations in terms of relaxation and electron transfer [19]. Here, we focus on neutral divacancy-based complexes as a representative model system, without claiming a unique or globally most stable aggregation pathway. In this work, we use first-principles calculations based on the density functional theory (DFT) [20, 21] to systematically examine the stability and electronic states of (VGaVN)𝑛 complexes in bulk GaN.  All calculations were performed using the large-scale DFT code "CONQUEST[22-24]" with the generalized gradient approximation (Perdew-Burke-Ernzerhof functional [25]). Pseudopotentials were generated by Optimized Norm-Conserving Vanderbilt Pseudopotential code [26, 27] with the parameters given in PseudoDojo[28], and PAO basis set of double-zeta plus polarization functions (DZP) consistent with the pseudopotentials were utilized[29,30]. Note that the 3d orbitals of Ga are included as the valence states. In this study, since our target systems contain more than 1000 atoms, the multisite support function method was used with the cutoff radius of 8.0 Bohr, to reduce the calculation cost[31]. We adopted a cutoff energy of 150 Ha for charge density grid, and used a 2×2×2 k-point mesh. The electronic self-consistency tolerance was set to 10^−8 e/Bohr^3 in the charge-density residual, resulting in the error of 10^-7 Ha in energy. All atomic positions were fully optimized until Hellmann-Feynman forces were smaller than 0.5 mHartree/Bohr.  To model large vacancy aggregation in the wurtzite GaN, we first define the defect formation energy for a (VGaVN)𝑛 complex as 𝐸form𝑛 = 𝐸𝐷 − (𝐸𝑃 − 𝑛𝐸Ga−N),       (1) where, 𝐸form𝑛  is the formation energy of the defect structure, 𝐸𝐷 is the total energy of the system with the defect,  𝐸𝑃 is the energy of a perfect wurtzite GaN crystal of the same size, and 𝐸Ga−N is the energy of one GaN unit in the perfect crystal. Using This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 3  Eq. (1), we assess the reliability of the model by examining the convergence of 𝐸form𝑛  with respect to supercell size. As test cases, we consider 𝑛 = 1 and 2 (Fig. 1, insets). Because vacancy-induced lattice relaxations can be long ranged, interactions with their own periodic images may affect the calculated energetics. Figure 1 shows that  𝐸form𝑛  converges as the supercell size approaches 1,000 atoms. Based on this convergence behavior, we adopt a 1600-atom wurtzite GaN supercell with the fixed lattice dimensions of 25.8885×28.0252×26.4373 Å for the configurational survey of larger n as shown in Fig. 2. The relaxed atomic coordinates of the representative 3A, 4D, and 6E structures are provided in XYZ format in the Supplementary Material. Using the basic supercell, we construct a set of (VGaVN)𝑛 models up to 𝑛 = 6. For each 𝑛, we consider several possibilities for the vacancy sites of Ga and N. Figure 3 shows examples of the candidate structures for 𝑛 = 3. To investigate the energetic stepwise growth, we define the energy gain by aggregation for 𝑛 ≥ 2 as 𝐸gain𝑛 = −{𝐸form𝑛 − (𝐸form𝑛−1 + 𝐸form1 )} .      (2) In Eq. (2), 𝐸form𝑛−1  refers to the formation energy of the lowest-energy (n-1)-mer obtained in the present survey. Accordingly, the reference structures used for the 𝑛 =2, 3, 4, 5, and 6 data points in Fig. 4(b) are 1C, 2D, 3A, 4D, and 5D, respectively. With this definition, a positive 𝐸gain𝑛  indicates that forming an 𝑛-mer from an (𝑛 − 1)-mer and an isolated VGaVN pair is energetically favorable. In the present study, Eqs. (1) and (2) are applied only to neutral (VGaVN)𝑛 complexes. For such complexes, the chemical-potential contribution reduces to 𝑛(𝜇Ga + 𝜇N) = 𝑛𝜇GaNbulk, and therefore the formation energy defined by Eq. (1) is independent of Ga-rich or N-rich growth conditions. In addition, because only neutral defects are considered, no explicit Fermi-level term appears in the present definition. Thus, Eqs. (1) and (2) should be regarded as restricted energetic measures within the present neutral stoichiometric model, not as fully general defect-thermodynamic quantities including arbitrary chemical potentials, charge states, and Fermi-level dependences. We examine the stability of neutral (VGaVN)𝑛 complexes over a wide range of candidate configurations by evaluating the formation energies 𝐸form𝑛  and energy gains by stepwise growth 𝐸gain𝑛 . Since the number of possible candidates increases rapidly with respect to 𝑛, our survey needs to be limited to the structural models which seem to have lower formation energies. In this respect, we mainly consider the vacancy models which have a small number of dangling bonds. Note that the reconstruction of surrounding lattice should be also important for the stability of (VGaVN)𝑛 complexes.  Figures 4(a) and 4(b) show 𝐸form𝑛  and 𝐸gain𝑛 , respectively, for candidate configurations considered for 𝑛  = 1–6. The calculated results show that, within the present restricted model, the stepwise aggregation gain is particularly large at 𝑛 = 3 and This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 4  𝑛 = 6, indicating enhanced binding for the addition of one extra VGaVN pair at these sizes. For n=1, we label the nearest-neighbor VGaVN divacancy structures by their orientation (c-axis or tilted) and by their local relaxation mode (A-like or B-like), where the A/B-like terminology refers to the presence/absence of pronounced outward breathing of Ga atoms neighboring the VN site[19]. For both orientations, two locally stable relaxed configurations are obtained. In the neutral state, the A-like configuration is lower in energy than the B-like one for both orientations. The two B-like configurations (c-axis and tilted) remain nearly degenerate (6.77 and 6.81 eV), whereas the A-like configurations are lower in energy (6.65 and 6.47 eV for the c-axis and tilted cases, respectively). The evaluated atomic displacements around the vacancy sites are summarized in Table 1. For n = 2, more configurations are possible and the formation energies vary within about 1 eV. The lowest-energy structure is 2D, and its formation energy is substantially lower than 2𝐸form1 , indicating that aggregation is already energetically favorable at this size. Here, 𝐸form1  denotes the formation energy of the lowest-energy 𝑛 = 1 divacancy used as the reference. When increasing 𝑛 by 1, the added Ga and N vacancy sites can be either nearest-neighbor (“bonded”) or separated (“non-bonded”). Notably, even among the three bonded configurations (2A–2C), the stability order depends on the supercell size (Fig. 1), indicating that local relaxation affects the relative energetics and that finite-size effects must be treated carefully even for small defects. Though the stabilization energy at 𝑛 = 2 is significant, the value for the aggregation from 𝑛 = 2 to 𝑛 = 3 is even larger. For 𝑛 = 3, the most stable trimer is 3A shown in Fig. 3(a) and its formation energy is much lower than the other three models, 3B, 3C and 3D, shown in Figs. 3(b)-(d), respectively. The remarkable stability of the 3A model comes from its small number of dangling bonds. The increase of the dangling bonds for 3A from 2D is five, while it is six for other models. Structurally, the 3A model is the first compact structure that can share the dangling bonds at two edges when one additional “bonded” Ga-N divacancy is introduced. Due to this great stability, we use 3A model as a starting configuration to generate the candidate structures beyond trimer (n > 3) and then add VGaVN unit one by one, as shown in Fig. 5 for 𝑛 = 4, and Fig. 6 for 𝑛 = 6. Considering the large stabilization energy at 𝑛 =3, together with the smaller energy gain at 𝑛 =4 and 5, it is expected that the trimer cluster (VGaVN)3 should be a major defect in Mg-implanted GaN after high-temperature annealing. The trap at 𝑛 =3 during the growth of vacancy clusters (VGaVN)𝑛  is also suggested from the kinetic effect: the population of VGaVN pair is expected to decrease as the aggregation proceeds, while the compact trimers should be less mobile than isolated divacancy pairs. The substantial population of (VGaVN)3 is consistent with the observation by PAS experiments. The stabilization energy at 𝑛 =6 is also large. Again, although the trimer–trimer coalescence (3+3→6) is energetically favorable (with an energy gain This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 5  of 6.34 eV), the actual growth pathway depends on kinetic factors such as cluster mobility and migration barriers. The 5+1 process provides a slightly larger energetic driving force (6.75 eV) in our dataset.  Figure 5 shows representative 𝑛 = 4 structures generated from 3A. Although 4A–4C have the same dangling-bond count, their formation energies differ substantially, indicating an important role of lattice relaxation. For 𝑛 = 5, the lowest-energy model, 5D, shown in Fig. 2(e), is obtained by adding the non-bonded unit to the remaining open edge of 4D. Its energy gain is close to that for 𝑛 = 4, despite the different change in dangling-bond number. Figure 6 shows the five candidate structures for 𝑛 = 6. Among these models, 6E has the lowest number of dangling bonds and is the most stable. The total number of dangling bonds in 6E is 18, which is six larger than that in 3A. Notably, 5D and 6E have the same total number of dangling bonds, despite the addition of one VGaVN unit. This small number of dangling bonds is probably the main reason why energy gain for this model is much larger than other cases. The different number of dangling bonds also results in the qualitative difference of the electronic structure between 6E and other models. The remarkable stability of the 6E model leads to the deduction that further growth after this hexamer may proceed by adding “bonded” VGaVN unit one by one to the three edges until 𝑛 = 9. During this process, the increase of the dangling bonds by one additional VGaVN is five and the energy gain should be larger. After that, only “non-bonded” VGaVN unit can be added up to 𝑛 =12, resulting in smaller 𝐸gain𝑛 s. Then, we can expect higher 𝐸gain𝑛 s again, after 𝑛 = 12. Such selective pathway may be relevant to the anisotropic void morphologies reported in GaN, including hexagonal-based pyramidal voids bounded by {0001} basal facets and {101-1} side facets. This extrapolated growth pathway toward n=9 and n=12, as well as its possible connection to pyramidal void formation, should be regarded as speculative. Configurational entropy, vibrational free-energy contributions, and other finite-temperature effects relevant to annealing at 800–1000 ℃ are not included in the present zero-temperature energetic analysis. Although the analysis based on the number of dangling bonds can roughly estimate the stability of vacancy clusters, such analyses do not consider the difference between the Ga- and N-related dangling bonds and the effect of local structure relaxation.  Therefore, we further analyzed the relaxed atomic configurations of representative structures 3A, 6E, and 6C. To quantify the local relaxation, we evaluated the displacement of atoms that are nearest neighbors (NN) of the vacancy sites in the perfect structure. The results are summarized in Table 1. Here, dangling bonds are counted as broken nearest-neighbor Ga–N bonds terminating at atoms remaining in the defective structure; bonds between two removed atoms are not counted as dangling bonds. In contrast, NN atoms used for the displacement analysis denote unique atoms that are NN of the vacancy sites in the corresponding perfect structure. The average displacement of vacancy-NN atoms is about 0.20 Å for 3A, 6E, and 6C, This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 6  which is not negligibly small. On the other hand, the change of bond lengths for non-NN atoms is less than 0.01 Å, indicating that the relaxation is strongly localized. The maximum NN displacements of the compact structures 3A and 6E are 0.358 and 0.378 Å, respectively, whereas 6C shows a larger maximum NN displacement of 0.626 Å. Since the average NN displacement of 6C remains close to those of 3A and 6E, this large maximum value reflects a localized distortion in the chain-like configuration rather than a uniform reconstruction of the defect. We also examined possible rebonding across the vacancy sites by checking local Ga–Ga and N–N distances around the defect core. For the compact structures 3A and 6E, the shortest distances remain 3.23–3.26 Å, close to the corresponding distance in the perfect lattice, 3.24 Å. Thus, no clear dimer-like Ga–Ga or N–N rebonding is found in these stable compact structures. For the 3A and 6E models, the inward/outward components of the NN displacements indicate that Ga atoms neighboring the vacancy sites show outward relaxation of about 0.33–0.34 Å on average, whereas the corresponding displacement of N atoms is much smaller. This outward relaxation of the neighboring Ga atoms is qualitatively consistent with the reported type-A breathing relaxation of the VGaVN divacancy[19]. Because the present systems are larger vacancy clusters, we use the type-A/type-B terminology only as a qualitative guide to the local relaxation, and we do not attempt a one-to-one assignment for every local environment in the clusters. Overall, the analysis based on the number of dangling bonds captures the main stability trend and the low-energy envelope, while morphology-dependent local relaxation and the distinct character of Ga- and N-related dangling-bond equivalents account for the remaining quantitative differences among compact structures. Figure 7 shows the total density of states (DOS) of (VGaVN)6 for the 6E and 6C models in Figs. 7(a) and 7(b), respectively. Here, we only discuss the qualitative difference in the electronic structure between 6E and 6C models. This difference can be qualitatively interpreted as follows. Since Ga and N vacancies are introduced in pairs throughout the present (VGaVN)𝑛 models, the numbers of Ga- and N-derived dangling bonds remain balanced, and the global electron-counting condition is satisfied for all the structures considered here. This condition alone, however, does not necessarily guarantee a gap-like electronic structure. In the compact 6E structure, local coordination and relaxation can facilitate electron transfer between Ga- and N-derived dangling-bond states, resulting in filled N-derived states and empty Ga-derived states, consistent with the gap-like DOS in Fig. 7(a). In contrast, in the chain-like 6C structure, dangling bonds are concentrated around the open-edge region, where such electron transfer and the resulting separation of occupied and unoccupied dangling-bond states are less effective. Consequently, finite DOS remains at the Fermi level, as shown in Fig. 7(b).The gap opening at n=6 is not unique to 6E: other hexamers, 6A and 6D, with bonded VGaVN pairs at the edges also show gap-like DOS, whereas the non-bonded or open-edge structures 6B and 6C exhibit gap closing. A similar relationship between local geometry and electronic structure is found for other (VGaVN)𝑛 clusters. As is well known, PBE does not accurately reproduce the absolute band gap or defect-level positions. Quantitative This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 7  analysis would require at least hybrid-functional calculations, which are beyond the scope of the present 1600-atom supercell study. The energy origin is set at the Fermi level, 𝐸F. Our results show that distinct configurations produce different electronic properties near 𝐸F. This suggests that high-resolution spectroscopic measurements might be able to identify and classify the dominant defect complexes, although the definite assignments may need the explicit treatment of charge states and excitonic effects. We performed large-scale first-principles calculations on neutral (VGaVN)𝑛 aggregates (𝑛 = 1–6) in bulk GaN. Aggregation is energetically favorable for all sizes, with particularly large energy gains at 𝑛 = 3 and 𝑛 = 6, consistent with PAS reports suggesting (VGaVN)3 as a major defect species after high-temperature annealing of Mg-implanted GaN. Compact structures tend to exhibit gap-like DOS, whereas chain-like structures can remain metallic-like, reflecting differences in local coordination and relaxation. These results provide a basis for discussing plausible growth pathways and motivate future kinetic studies of migration barriers between competing configurations. Although charged states may modify the absolute formation energies in p-type or Mg-implanted material, the enhanced stability at n=3 and n=6 is expected to remain qualitatively robust because it mainly originates from local coordination effects, such as dangling-bond sharing in compact geometries. In experiments, many vacancy-related defects are generated by Mg implantation.  However, it should be noted that a substantial fraction of the implanted Mg atoms can become electrically inactive by forming specific defect structures, such as pyramidal inversion domains[32]. Thus, the present neutral stoichiometric model still provides a useful reference for their intrinsic local stability.  Supplementary Material The supplementary material contains the relaxed atomic coordinates of the representative 3A, 4D, and 6E structures in XYZ format.  The present study was partially supported by JSPS KAKENHI Grant Number 23K04602 and 23H01871. The calculations were carried out on Numerical Materials Simulator at NIMS, and on the supercomputer Fugaku provided by the RIKEN Center for Computational Science (hp250292). This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 8   1H. Amano, Y. Baines, M. Beam, et al., J. Phys. D: Appl. Phys. 51, 163001 (2018). 2L. Nela, M. Xiao, Y. Zhang, and E. Matioli, Appl. Phys. Lett. 120, 190501 (2022). 3K. H. Teo, Y. Zhang, N. Chowdhury, S. Rakheja, R. Ma, Q. Xie, E. Yagyu, K. Yamanaka, K. Li, and T. Palacios, J. Appl. Phys. 130, 160902 (2021). 4M. Meneghini, C. D. Santi, I. Abid, M. Buffolo, M. Cioni, R. A. Khadar, L. Nela, N. Zagni, A. Chini, F. Medjdoub, G. Meneghesso, G. Verzellesi, E. Zanoni, and E. Matioli, J. Appl. Phys. 130, 181101 (2021). 5M. Haziq, S. Falina, A. A. Manaf, H. Kawarada, and M. Syamsul, Micromachines 13, 2133, (2022). 6G. Meneghesso, M. Meneghini, and E. Zanoni, Jpn. J. Appl. Phys. 53, 100211 (2014). 7T. Kachi, T. Narita, H. Sakurai, M. Matys, K. Kataoka, K. Hirukawa, K. Sumida, M. Horita, N. Ikarashi, K. Sierakowski, M. Bockowski, and J. Suda, J. Appl. Phys. 132, 130901 (2022). 8A. Uedono, H. Sakurai, T. Narita, K. Sierakowski, M. Bockowski, J. Suda, S. Ishibashi, S. F. Chichibu, and T. Kachi, Sci. Rep. 10, 17349 (2020). 9F. Roccaforte, F. Giannazzo, and G. Greco, Micro 2, 23 (2022). 10T. Tsujio, M. Oda, and Y. Shinozuka, Jpn. J. Appl. Phys. 56, 091001 (2017). 11M. Oda, Jpn. J. Appl. Phys. 58, SCCC16 (2019). 12J. Nakamura, M. Oda, and Y. Kangawa, Phys. Status Solidi B 261, 2400026 (2024). 13F. Tuomisto and I. Makkonen, Rev. Mod. Phys. 85, 1583 (2013). 14A. Uedono, S. Takashima, M. Edo, K. Ueno, H. Matsuyama, W. Egger, T. Koschine, C. Hugenschmidt, M. Dickmann, K. Kojima, S. F. Chichibu, and S. Ishibashi, Phys. Status Solidi B 255, 1700521 (2018). 15A. Uedono, H. Iguchi, T. Narita, K. Kataoka, W. Egger, T. Koschine, C. Hugenschmidt, M. Dickmann, K. Shima, K. Kojima, S. F. Chichibu, and Shoji Ishibashi, Phys. Status Solidi B 256, 1900104 (2019). 16A. Uedono, R. Tanaka, S. Takashima, K. Ueno, M. Edo, K. Shima, S. F. Chichibu, J. Uzuhashi, T. Ohkubo, S. Ishibashi, K. Sierakowski, and M. Bockowski, Phys. Status Solidi B 261, 2400060 (2024). 17A. B. Yankovich, A. V. Kvit, X. Li, F. Zhang, V. Avrutin, H. Y. Liu, N. Izyumskaya, Ü. Özgür, H. Morkoç, and P. M. Voyles, J. Appl. Phys. 111, 023517 (2012). 18H. Li, M. Huang, and S. Chen, J. Semicond. 41, 032104 (2020). 19Y. Gohda and A. Oshiyama, J. Phys. Soc. Jpn. 79, 083705 (2010). 20P. Hohenberg and W. Kohn, Phys. Rev. 135, B864 (1964). 21W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 22D. R. Bowler, T. Miyazaki and M. J. Gillan, J. Phys. Condens. Matter 14, 2781 (2002). 23T. Miyazaki, D. R. Bowler, R. Choudhury and M. J. Gillan, J. Chem. Phys. 121, 6186 (2004). 24A. Nakata, J. S. Baker, S. Y. Mujahed, J. T. L. Poulton, S. Arapan, J. Lin, Z. Raza, S. Yadav, L. Truflandier, T. Miyazaki, and D. R. Bowler, J. Chem. Phys. 152, 164112 (2020). 25J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 26D. R. Hamann, Phys. Rev. B 88, 085117 (2013). 27M. Schlipf and F. Gygi, Comput. Phys. Commun. 196, 36 (2015). 28M. J. van setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, Comput. Phys. Commun. 226, 39 (2018). 29A. S. Torralba, M. Todorović, V. Brázdová, R. Choudhury, T. Miyazaki, M. J. Gillan and D. R. Bowler, J. Phys.: Condens. Matter 20, 294206 (2008). This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 9  30D. R. Bowler, J. S. Baker, J. T. L. Poulton, S. Y. Mujahed, J. Lin, S. Yadav, Z. Raza and T. Miyazaki, Jpn. J. Appl. Phys. 58, 100503 (2019). 31A. Nakata, D. R. Bowler, and T. Miyazaki, Phys. Chem. Chem. Phys. 17, 31427 (2015). 32A. R. Persson, A. Papamichail, V. Darakchieva, and P. O. Å. Persson, Sci. Rep. 12, 17987 (2022).         FIG. 1. Convergence of the defect formation energy 𝐸form𝑛  for 𝑛 = 1 (1A) and 𝑛 = 2 (2A, 2B and 2C) with respect to the supercell size. Insets show representative geometries of the  (VGaVN)1 and (VGaVN)2 models used for the convergence tests.  This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 10    FIG. 2  (a) Atomic structure of 1600-atom wurtzite GaN supercell adopted as the basic model. (b–f) Lowest-energy structures of the (VGaVN)𝑛 complexes for 𝑛 = 2–6. Green and gray spheres denote Ga and N atoms, respectively. Red spheres indicate vacancy sites (VGa and VN).    This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 11   FIG. 3. Candidate structures of (VGaVN)3 viewed along the [11-20] direction.       This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 12   FIG. 4. (a) Formation energies 𝐸form𝑛  and (b) Energy gains 𝐸gain𝑛  for the (VGaVN)𝑛  complexes. Symbols represent distinct configurations at each 𝑛.       This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 13   FIG. 5. Candidate structures for (VGaVN)4 generated by adding one VGaVN pair to the trimer model (3A) of (VGaVN)3. Blue spheres indicate the added VGaVN pair, while red spheres denote the vacancies belonging to the original trimer. Panels (a)–(d) correspond to models 4A–4D, respectively.        This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 14   FIG. 6. Candidate structures for the 𝑛 = 6 complexes (VGaVN)6: (a) 6A, (b) 6B, (c) 6C, (d) 6D, and (e) 6E. In each panel, the trimer core (3A) is highlighted in red, and the three added VGaVN units are highlighted.       This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898 15   FIG. 7. Total densities of states (DOS) for the two structure models of (VGaVN)6: (a) 6E (compact, most stable) and (b) 6C (chain-like). Energies are referenced to the Fermi level (𝐸𝐹 = 0).  .Table 1. Average and maximum atomic displacements for the representative defect structures 1A-1D (1A and 1C: A-like; 1B and 1D: B-like), 3A, 6E, and 6C. The number of dangling bonds (DBs) is also listed. model No. of DBs avg. disp. of NN(Å) max disp. of NN(Å) avg. disp. of non-NN(Å) max disp. of non-NN (Å) 1A 6 0.220 0.337 0.0049 0.087 1B 6 0.138 0.167 0.0029 0.044 1C 6 0.247 0.407 0.0047 0.089 1D 6 0.089 0.166 0.0024 0.037 3A 12 0.197 0.358 0.0053 0.100 6E 18 0.199 0.378 0.0070 0.119 6C 22 0.202 0.626 0.0076 0.126    This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0331898