Charge transition levels and stability of Ni- and Ir-doped β-Ga₂O₃: a comprehensive hybrid functional study

Abstract

In this study, the optimized hybrid functional HSE(0.26,0.0) is employed to investigate the incorporation of nickel (Ni) and iridium (Ir) dopants in β-Ga₂O₃. The formation energies and charge transition levels of Ni and Ir at gallium sites are calculated. The results show that Ni prefers substitution at the octahedral (Ga₂) site, with a formation energy approximately 1 eV lower than at the tetrahedral (Ga₁) site. Ni at the Ga₂ site (Ni_Ga2) exhibits both donor and acceptor behaviors, with charge transition levels at (+1/0) 1.0 eV and (0/−1) 2.24 eV above the VBM, respectively. Ir similarly favors the octahedral site, displaying donor behaviors with charge transition levels at (+2/+1) 1.04 eV and (+1/0) 3.15 eV above the VBM. Our computational findings for the charge transition levels of Ni and Ir ions are in good agreement with recent experimental measurements, and they explain the correlation between Ni³⁺ and Ir⁴⁺ ion concentrations observed in electron paramagnetic resonance studies. Additionally, the calculated vertical transitions at 2.56 eV and 4.25 eV for Ni_Ga2, and at 2.91 eV and 4.62 eV for Ir_Ga2, below the conduction band minimum, are in good agreement with optical absorption results, confirming the presence of Ni and Ir substitutions at the Ga₂ site in β-Ga₂O₃. These computational results provide a detailed understanding of the behavior of Ni- and Ir-doped β-Ga₂O₃, highlighting the potential applications of Ni- and Ir-doped β-Ga₂O₃ for optoelectronic devices.

---

1. Introduction

The monoclinic phase of gallium oxide (β-Ga₂O₃), which is the most thermodynamically stable phase of Ga₂O₃, has a large band gap of approximately 5 eV.^1,2 In recent years, β-Ga₂O₃ has received a lot of attention and is considered a new competitor to 4H–SiC and GaN materials. In addition to its high thermal and chemical stability, the wide to ultra-wide band gap of β-Ga₂O₃, resulting in a high breakdown field of around 8 MV cm⁻¹, makes it a potential candidate for applications in high power electronics.^3,4 The availability of high-quality single crystals and thin films of β-Ga₂O₃ holds significant promise for applications in gas sensors,^5,6 electroluminescent devices,^7,8 and photocatalysis.^9–11

Unintentional dopants such as Si and H are often present in β-Ga₂O₃ crystals, leading to a moderate n-type doping concentration of around mid-10¹⁷ cm⁻³.^12–14 Additionally, Sn, Ge, and Nb impurities are also used to create n-type β-Ga₂O₃ crystals.^15–18 However, fabricating p-type crystals of β-Ga₂O₃ has proven to be very challenging, which limits its potential applications. Recently, considerable attention has been directed towards identifying suitable p-type dopants for β-Ga₂O₃. Despite extensive efforts, p-type single crystals of β-Ga₂O₃ have not yet been successfully produced. Various metal acceptors, such as Zn, Mg, and Ni, have been investigated as potential dopants in β-Ga₂O₃, but the results suggested deep acceptor behaviors of the dopants. Among these dopants, Ni-doped β-Ga₂O₃ samples show promise for producing semi-insulating substrates for lateral power devices. Gustafson et al.¹⁹ have grown a Ni-doped β-Ga₂O₃ crystal by the Czochralski method using an iridium crucible and radio frequency heating. As a result, the presence of Ir impurity is unavoidable. In fact, by using electron paramagnetic resonance and optical absorption, the presence of Ni and Ir substitution at gallium sites was reported by Gustafson et al.¹⁹ The optical absorption peaks at 303 nm and 442 nm were assigned to Ni at the octahedral site. When irradiated at 275 nm, the concentration of Ni³⁺ increases and the Ir⁴⁺ concentration decreases. After increasing the temperature above 375 °C, the concentration of Ni³⁺ and Ir⁴⁺ is restored. In addition, the Ni acceptor level is reported at 1.4 eV above the VBM. In the report of Seyidov et al.²⁰ the acceptor level of Ni was reported at 1.9 eV and a deep donor level at 1.1 eV above the VBM, which was not reported by Gustafson et al.¹⁹ In the case of unintentional doping with Ir, Ir_Ga was reported to be an n-type donor, the EPR study indicated nonmagnetic (S = 0) Ir³⁺, while Ir⁴⁺ was reported as S = 1/2 in ref. 19 and 21. In ref. 22, Ir exhibited both double and single donor behaviors, whereas in the study by Zachinskis and co-workers,²³ only the single donor state of Ir at the Ga site was observed.

While standard density functional theory (DFT) and its modified versions, such as DFT + U and DFT + U + V, have played a significant role in understanding the electrical and optical properties of semiconductors.^24–27 They are not enough when applied to defects in wide band-gap materials.^28,29 In recent years, hybrid functionals have emerged as a promising alternative for studying wide band-gap semiconductors, particularly screened hybrids^30,31 with the HSE06 functional developed by Heyd, Scuseria, and Ernzerhof being especially notable^32,33 More recently, the optimized hybrid functional calculation method has demonstrated its advantages in studying defects in various wide-bandgap semiconductors, including GaN, CuGaS₂, and cubic-ZnS.^34–36 For β-Ga₂O₃, the optimized hybrid functional with parameters (α = 0.26, μ = 0.00) has already shown its advantages in reproducing and explaining experimental results both with intrinsic defects and,^13,37,38 with the substitution of metals at the Ga site.^39,40 The procedure for obtaining the optimized parameters that reproduce the gap-optimized, Koopmans' theorem-compliant hybrid functional for β-Ga₂O₃ has been discussed in the section III of ref. 37. Although, numerous studies have investigated Ni and Ir in β-Ga₂O₃, a comprehensive understanding of these two impurities considering various aspects such as formation energy, spin states, and the Fermi-level remains incomplete. In this research, the optimized hybrid functional will be employed to investigate Ni- and Ir-doped β-Ga₂O₃. In this paper, HSE(0.26,0.00) will be used to investigate Ni- and Ir-doped β-Ga₂O₃. By comparing our calculated results with previous experimental and simulation studies, this research aims to provide a complete understanding of Ni and Ir impurities in β-Ga₂O₃. This includes examining formation energies, the most favorable spin states of defects, and correlating Ni and Ir concentrations with EPR spectra, optical absorption, and predicting luminescence emission peaks associated with Ni and Ir substitution in β-Ga₂O₃.

2. Computational details

In this report, the gap-optimized, Koopmans' theorem-compliant hybrid functional HSE, with 26% Hartree–Fock exchange (α = 0.26) and a conventional hybrid functional (μ = 0.00), where exchange interactions are treated uniformly, was employed. Calculations were carried out using the Vienna Ab initio Simulation Package (VASP 5.4.4) along with the projector augmented wave (PAW) method.^41–43 The semi-core d electrons of gallium were treated as part of the valence shell. A plane-wave basis set cutoff was set at 420 eV, with a kinetic energy cutoff for augmentation charges at 840 eV. In this study, Ni and Ir are doped into the 3D bulk structure of β-Ga₂O₃. A 160-atom supercell, constructed as a 1 × 4 × 2 expansion of the base-centered monoclinic unit cell, was used for defect calculations. The supercell has dimensions of approximately 12.2 Å × 9.1 Å × 11.6 Å, and the lattice parameters adopted for this work are a = 12.25 Å, b = 3.04 Å, c = 5.82 Å, and β = 103.8°.³⁷ The self-consistent electronic energy was converged to 10⁻⁴ eV, with ion relaxation continuing until forces on each ion were reduced to below 0.02 eV Å⁻¹. Calculations employed the Γ-approximation for Brillouin zone sampling, proven to be convergent within 0.1 eV for this supercell size in the previous work.⁴⁴ To mitigate artificial interactions between repeated images of charged defects, the charge correction method from ref. 45 and a high-frequency dielectric constant of ε_∞ = 3.55 (ref. 2) have been used in our calculations.

3. Results and discussion

In β-Ga₂O₃, there are two nonequivalent gallium sites: one with tetrahedral coordination and one with octahedral coordination. The formation energies of nickel and iridium incorporation on the 2 inequivalent gallium sites in different charge states are shown in Fig. 1, as a function of the Fermi-level position between the calculated valence band maximum (VBM) and the conduction band minimum (CBM), using the following expression:where E_tot[X^q] and E_tot[Ga₂O₃] represent the total energies of the structures containing the substitution X (Ni and Ir in the current research) at charge state q and the pristine β-Ga₂O₃ supercell, respectively. n_i indicates the number of atoms added to (n_i > 0) or removed from (n_i < 0) the pristine supercell to form the defect. E^q_corr is the correction for the finite-size supercell approach as described in the computational details. μ_i represents the chemical potential of the atoms in the reservoir corresponding to the crystal growth conditions. The E_F is the Fermi-level position between the VBM and CBM. The chemical potential of μ_Ni is bounded by the formation of NiGa₂O₄ in the O-rich limits (ΔH[NiGa₂O₄] = −13.79 eV, the calculated heat of formation of NiGa₂O₄) and Ga₃Ni₂ in the Ga-rich limit (ΔH[Ga₃Ni₂] = −2.24 eV). While in the case of iridium the chemical potential μ_Ir is set by Ga₉Ir₂ in the Ga-rich condition (ΔH[Ga₉Ir₂] = −4.95 eV) and IrO₂ in the O-rich condition (Δ[IrO₂] = −2.92 eV).

Fig. 1 Formation energy (eV) of Ni_Ga and Ir_Ga under as a function of the Fermi-level position.

In Fig. 1, the red lines represent the charge transition levels of Ni incorporating on Ga sites (Ni_Ga). The results indicate that Ni substitution at the octahedral site (Ni_Ga2) is energetically preferred over the tetrahedral site (Ni_Ga1) across all charge states. This result is in good agreement with previous studies.^19,20 Ni_Ga show deep acceptor behaviors with low formation energies for both sites under n-type conditions. The calculated acceptor levels (0/−1) are 2.02 eV and 2.30 eV above the VBM for Ni_Ga2 and Ni_Ga1, respectively. Our calculated results are higher in energy compared to the theoretical charge transition levels in the work of Seyidov et al.²⁰ For the donor states (0/+1), the transition levels are found at 1.26 eV and 1.48 eV above the VBM for Ni_Ga2 and Ni_Ga1, respectively. The differences in the charge transitions levels in our calculations and those in ref. 20 are likely from the choice of using the parameters in hybrid functional and dielectric constant used in charge correction.

In this study, it was also found that substitutional Ni at both Ga sites is preferred in the spin-triplet states in the acceptor state, and the spin-singlet state is more favorable in the neutral charge state. This finding aligns well with the research of Seyidov and co-workers.²⁰ In the donor state, Ni_Ga2 has the lowest formation energy in the singlet spin state, while Ni_Ga1 is stable in the spin-triplet state compared to the non-magnetic singlet spin state which is 0.65 eV higher in formation energy. The spin-triplet state of Ni_Ga1 in the donor charge state is likely caused by the difference in ionic radii between tetrahedral Ga³⁺ (47 pm) and Ni²⁺ (69 pm). The larger substitutional Ni ion at the tetrahedral Ga site causes a distortion that likely triggers the Jahn–Teller effect which removes the degeneracy and lowers the system's energy, thus achieving a more stable configuration in the triplet spin state. The singlet spin state of Ni_Ga2 in the +1 charge state aligns with the findings of ref. 20, but the donor state of Ni_Ga1 site was reported to be stable in spin-quartet state in that report.²⁰

Fig. 1 also displays the formation energy of Ir incorporating at Ga sites (Ir_Ga), by the blue lines. Similar to Ni substitution, Ir substitution at the octahedral Ga site is energetically preferred compared to the tetrahedral Ga site due to the ionic radius of Ir (68 pm). The octahedral Ga ion (62 pm) provides more space to accommodate the larger ionic radius of Ir compared to the tetrahedral sites. Ir substitution at the Ga₂ site (Ir_Ga2) minimizes lattice strain and leads to more stable dopant incorporation. Ir_Ga2 shows donor behavior, with charge transition levels at 0.94 eV (+2/+1) and 2.9 eV (+1/0). For substitution at the tetrahedral site, the donor states have charge transition levels at 1.35 eV and 3.03 eV for (+2/+1) and (+1/0), respectively. It was also found that Ir_Ga1 has an acceptor state with a charge transition (0/−1) occurring in the band gap at 4.61 eV above the VBM, the compensated acceptor level of Ir_Ga1 was also reported by Ritter et al.²² Our calculated results for Ir_Ga charge transition levels and formation energies are in good agreement at the (+2/+1) transition, but show deeper energy at (+1/0) compared to other studies.^22,23

It should be noted that while Ritter and coworkers²² reported all charge states, the authors did not specify the spin state of the substitution. On the other hand, Zachinskis et al.²³ demonstrated that the high spin state S = 4/2 has the lowest formation energy for Ir_Ga2 in the neutral charge state, and S = 3/2 is energetically favorable for Ir_Ga1 in the singly donor state. However, no information on the double donor state was provided. In our calculation, for the neutral charge state, Ir_Ga2 adopts the singlet spin state, while Ir_Ga1 has the lowest formation energy in the spin-triplet state, which is 0.25 eV lower than that of the singlet spin state. It is due to the distortion caused by the larger Ir dopant at the Ga₁ site and the Jahn–Teller effect, as discussed above with Ni_Ga1 in the donor charge state. In the singly donor states, Ir_Ga at both substitution sites is most stable in the S = 1/2 spin state. The most favorable state for Ir_Ga2 in the doublet spin state is in good agreement with the observation of Ir⁴⁺ in EPR measurements.^19,21

Experimental studies have shown that Ir is often present in β-Ga₂O₃ as an unintentional dopant, especially in Ni-doped β-Ga₂O₃ crystals grown by the Czochralski method.^22,46,47 EPR measurements by Gustafson et al.¹⁹ have identified the signals of Ni and Ir substitution at Ga sites and their concentration correlation. Specifically, when a Ni-doped sample is irradiated with 275 nm light, the concentration of Ni³⁺ ions increase, while the concentration of Ir⁴⁺ ions decrease. The concentration of Ni²⁺ can be restored by heating above 375 °C.¹⁹ The above EPR results can be explained as follows: when the sample is doped with an acceptor dopant like Ni, the correlation between unintentional n-type dopants (such as Si and H) and Ni pins the Fermi-level around the middle of the band gap, similar Fermi-level range was also reported in ref. 20. If one considers the Fermi-level as in the shaded orange row of Fig. 2, Ni_Ga2 and Ir_Ga2 are energetically favorable in the −1 (Ni²⁺) and +1 (Ir⁴⁺) charge states, respectively. When the sample is irradiated with 275 nm light, electrons are excited to the conduction band, leading to Ni_Ga2 and Ir_Ga2 are identified at 0 (Ni³⁺) and +2 (Ir⁵⁺) charge states. Thus, the EPR measurement shows an increase of Ni³⁺ ion and a decrease in Ir⁴⁺ ion concentration after irradiation. When the sample is heated, electrons recombine with holes, that witness the restoration of Ni²⁺ and Ir⁴⁺ concentrations.

Fig. 2 The schematic band diagram illustrates the defect-level energies for Ni_Ga and Ir_Ga defects. The shaded orange row indicates the estimated Fermi-energy range.

Fig. 3 illustrates the localization of electrons at the spin state S = 1/2 of Ni_Ga2⁰ and Ir_Ga2¹⁺, corresponding to Ni³⁺ and Ir⁴⁺, as addressed in the experimental studies.¹⁹ The wave functions of electrons in this spin state are not only localized at the substitutional site but also at several oxygen sites around the substitutional sites. This spread of wave function localization explains the broad EPR lines of Ni³⁺ and Ir⁴⁺ ions in the experimental observation.^19,21 In the analysis of Ni-doped β-Ga₂O₃ crystals, Gustafson et al.¹⁹ also found broad optical absorption bands at 442 nm (2.81 eV) and 303 nm (4.09 eV), which were attributed to Ni substitutions. Experimental data of optical absorption peaks are often interpreted as vertical electron transitions from the defect to the CBM.^29,37 In our calculations, the vertical charge transitions of Ni at (+1/0) and (0/−1) are 2.56 eV and 4.25 eV below the CBM, respectively. This is in good agreement with the experimental results in ref. 19. In the case of Ir_Ga2, the vertical transition levels are calculated at 2.92 eV and 4.54 eV for the neutral and +1 charge states, respectively. The vertical transition of the neutral charge state is in line with the study by Ritter et al.,⁴⁶ which reports an absorption at 2.8 eV. Further details on the comparison of the optical absorption and the calculated vertical transition energies for Ni and Ir substitutions are provided in Table 1.

Fig. 3 Diagram of (a) Ni and (b) Ir substitution at Ga₂ site. The green and red balls represent Ga, O, respectively. The blue lobes show the localization of electron (a) Ni_Ga2⁰ and (b) Ir_Ga2¹⁺ in the spin-singlet state.

Table 1 Charge transition levels with respect to the calculated VBM and CBM (eV), in comparison with values in recent studies in both experimental and HSE calculations

Defects	Transition levels	Experimental transition	Recent HSE calculations	This work HSE(0.00,0.26)	Optical absorption	Recent HSE calculations	This work HSE(0.00,0.26)
NiGa2	+1/0	X	VBM +1.09 (ref. 20)	VBM +1.26	CBM −4.09 (ref. 19)	CBM −3.8 (ref. 20)	CBM −4.25
					CBM −2.81 (ref. 19)
	0/−1	VBM +2.0 (ref. 20)	VBM +1.88 (ref. 20)	VBM +2.02		CBM −3.16 (ref. 20)	CBM −2.56
		VBM +1.4 (ref. 19)
IrGa2	+2/+1	CBM −2.2 (ref. 46)	VBM +1.01 (ref. 22)	VBM +0.94	CBM −2.9, −3.5, −4.4 (ref. 48)	X	CBM −4.54
	+1/0		CBM −2.25 (ref. 22) (VBM +2.6)	CBM −2.1 (VBM +2.9)	CBM −2.8 (ref. 22)	CBM −2.99 (ref. 23)	CBM −2.92
			CBM −2.58 (ref. 23) (VBM +2.15)

---

We also predicted the emission energy, which corresponds to the emitted wavelength from the recombination of an electron with a hole. The recombination was calculated considering a limited-size supercell, resulting in an accuracy of 100 meV for the photoluminescence energy. The difference between an electron in a shallow donor state and a conduction band electron bound by the trapped hole has not been considered. Therefore, an electron is considered to be at the CBM. The equilibrium geometry of the bound exciton was relaxed under the constraint of the orbital occupations, and the spin multiplicity was kept constant during the entire process. As explained above, the correlation between unintentional dopants and Ni acceptors places the Fermi-level in the middle of the gap. Therefore, under irradiation, the electron of donor Ni_Ga2 is excited to the CBM and then recombines with the hole. In the case of Ir_Ga2 this process occurs in the +1 charge state. The calculated recombination energy between a trapped hole at Ni²⁺ and an electron at the CBM was found to be 2.23 eV, while for Ir⁴⁺, it is 2.55 eV. The luminescence of Ir_Ga2 has yet to be confirmed, but it is likely due to the inadvertent Ir impurity.

4. Conclusion

This study shows a detailed theoretical investigation of Ni and Ir dopants in β-Ga₂O₃ using the optimized hybrid functional HSE(0.26,0.00), which closely agree with experimental results. The calculations results reveal that both Ni and Ir exhibit a strong preference for substitution at the octahedral gallium site. Ni dopant demonstrates dual donor and acceptor behaviors at charge transition levels (+1/0) 1.0 eV and (0/−1) 2.24 eV above the VBM, Ir impurity shows donor behaviors with charge transition levels (+2/+1) at 1.04 eV and (+1/0) at 3.15 eV above the VBM. These results align well with recent experimental measurements, including the correlation between Ni and Ir ion concentrations observed in electron paramagnetic resonance studies. Moreover, the calculated vertical transition levels for Ni_Ga2 (2.56 eV and 4.25 eV below the CBM) and Ir_Ga2 (2.91 eV and 4.62 eV below the CBM) correspond closely with experimental optical absorption peaks. This comprehensive analysis advances the understanding of Ni- and Ir-doped β-Ga₂O₃ and underscores their potential applications in optoelectronic devices, particularly in high-resistivity and semi-insulating substrates.

Data availability

The data that supports the findings of this study are available within the article.

Author contributions

Quoc-Duy Ho: conceptualization, investigation, data curation, methodology, formal analysis, visualization, writing – original draft. K. Dien Vo, Nguyen Thanh Tien, Huynh Anh Huy, and Duc-Quang Hoang: investigation, data curation, formal analysis, visualization, writing – editing. Duy Khanh Nguyen: supervision, investigation, data curation, validation, formal analysis, manuscript review and editing.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Duy Khanh Nguyen acknowledges the supports of Van Lang University (VLU), Ho Chi Minh City, Vietnam.

References

1. M. Orita, H. Ohta, M. Hirano and H. Hosono, Appl. Phys. Lett., 2000, 77, 4166.
2. J. Furthmüller and F. Bechstedt, Phys. Rev. B, 2016, 93, 115204.
3. S. J. Pearton, J. Yang, P. H. I. V. Cary, F. Ren, J. Kim, M. J. Tadjer and M. A. Mastro, Appl. Phys. Rev., 2018, 5, 011301.
4. J. Zhang, J. Shi, D.-C. Qi, L. Chen and K. H. L. Zhang, APL Mater., 2020, 8, 020906.
5. J. Zhu, Z. Xu, S. Ha, D. Li, K. Zhang, H. Zhang and J. Feng, Materials, 2022, 15, 7339.
6. A. Afzal, J. Mater., 2019, 5, 542.
7. M. M. Afandi, S. Jeong and J. Kim, Opt. Mater., 2023, 146, 114613.
8. T. Minami, T. Shirai, T. Nakatani and T. Miyata, Jpn. J. Appl. Phys., 2000, 39, L524.
9. Q. Li, et al., Catal. Commun., 2019, 120, 23.
10. Y. Hou, L. Wu, X. Wang, Z. Ding, Z. Li and X. Fu, J. Catal., 2007, 250, 12.
11. E. Filippo, M. Tepore, F. Baldassarre, T. Siciliano, G. Micocci, G. Quarta, L. Calcagnile and A. Tepore, Appl. Surf. Sci., 2015, 338, 69.
12. J. B. Varley, J. R. Weber, A. Janotti and C. G. Van de Walle, Appl. Phys. Lett., 2010, 97, 142106.
13. N. T. Son, Q. D. Ho, K. Goto, H. Abe, T. Ohshima, B. Monemar, Y. Kumagai, T. Frauenheim and P. Deák, Appl. Phys. Lett., 2020, 117, 032101.
14. P. D. C. King, I. McKenzie and T. D. Veal, Appl. Phys. Lett., 2010, 96, 062110.
15. S. C. Siah, et al., Appl. Phys. Lett., 2015, 107, 252103.
16. A. Mauze, Y. Zhang, T. Itoh, E. Ahmadi and J. S. Speck, Appl. Phys. Lett., 2020, 117, 222102.
17. F. Alema, G. Seryogin, A. Osinsky and A. Osinsky, APL Mater., 2021, 9, 091102.
18. W. Zhou, C. Xia, Q. Sai and H. Zhang, Appl. Phys. Lett., 2017, 111, 242103.
19. T. D. Gustafson, N. C. Giles, B. C. Holloway, J. Jesenovec, B. L. Dutton, J. S. McCloy, M. D. McCluskey and L. E. Halliburton, J. Appl. Phys., 2022, 132, 185705.
20. P. Seyidov, J. B. Varley, J.-X. Shen, Z. Galazka, T.-S. Chou, A. Popp, M. Albrecht, K. Irmscher and A. Fiedler, J. Appl. Phys., 2023, 134, 205701.
21. C. A. Lenyk, et al., J. Appl. Phys., 2019, 125, 045703.
22. J. R. Ritter, J. Huso, P. T. Dickens, J. B. Varley, K. G. Lynn and M. D. McCluskey, Appl. Phys. Lett., 2018, 113, 052101.
23. A. Zachinskis, J. Grechenkov, E. Butanovs, A. Platonenko, S. Piskunov, A. I. Popov, J. Purans and D. Bocharov, Sci. Rep., 2023, 13, 8522.
24. Y. Z. Abdullahi, R. Caglayan, A. Mogulkoc, Y. Mogulkoc and F. Ersan, JPhys Mater., 2023, 6, 025003.
25. A. Akyildiz, I. Ilgaz Aysan, Y. Z. Abdullahi, B. Akgenc Hanedar, Z. Demir Vatansever and F. Ersan, Phys. Chem. Chem. Phys., 2024, 26, 26566.
26. Y. Z. Abdullahi, Comput. Condens. Matter, 2021, 29, e00614.
27. Y. Z. Abdullahi and F. Ersan, RSC Adv., 2023, 13, 3290.
28. W. R. L. Lambrecht, Phys. Status Solidi B, 2011, 248, 1547.
29. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van de Walle, Rev. Mod. Phys., 2014, 86, 253.
30. S. J. Clark and J. Robertson, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 085208.
31. X. Zheng, A. J. Cohen, P. Mori-Sánchez, X. Hu and W. Yang, Phys. Rev. Lett., 2011, 107, 026403.
32. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207.
33. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
34. P. Deák, M. Lorke, B. Aradi and T. Frauenheim, Phys. Rev. B, 2019, 99, 085206.
35. M. Han, Z. Zeng, T. Frauenheim and P. Deák, Phys. Rev. B, 2017, 96, 165204.
36. Q. Duy Ho and M. Castillo, Comput. Mater. Sci., 2023, 216, 111827.
37. P. Deák, Q. Duy Ho, F. Seemann, B. Aradi, M. Lorke and T. Frauenheim, Phys. Rev. B, 2017, 95, 075208.
38. Q. D. Ho, T. Frauenheim and P. Deák, Phys. Rev. B, 2018, 97, 115163.
39. Q. D. Ho, D. K. Nguyen and H. A. Huy, Comput. Condens. Matter, 2022, 32, e00727.
40. Q. D. Ho, T. Frauenheim and P. Deák, J. Appl. Phys., 2018, 124, 145702.
41. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
42. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
43. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251.
44. P. Deák, Q. D. Ho, F. Seemann, B. Aradi, M. Lorke and T. Frauenheim, Phys. Rev. B, 2017, 95, 075208.
45. C. Freysoldt, J. Neugebauer and C. G. Van de Walle, Phys. Rev. Lett., 2009, 102, 016402.
46. J. R. Ritter, K. G. Lynn and M. D. McCluskey, J. Appl. Phys., 2019, 126, 225705.
47. Z. Galazka, et al., ECS J. Solid State Sci. Technol., 2017, 6, Q3007.
48. J. Jesenovec, J. Varley, S. E. Karcher and J. S. McCloy, J. Appl. Phys., 2021, 129, 225702.

---