Hybrid density functional theory has become a standard method for calculations of defects in semiconductors. The majority of work in this field is done using hybrid functionals tuned to reproduce the experimental bandgap of the host material. This approach usually yields results in reasonable agreement with the experiment. Alternatively, hybrid functional can be tuned to fulfill the generalized Koopmans' condition for defect orbitals, which cancels self-interaction energy and restores the linear behavior of energy with respect to electron occupation. Here, we investigate the methods of hybrid functional tuning, which both satisfy the generalized Koopmans' condition and reproduce the experimental bandgap, using one of the most well-studied defects in GaN, carbon acceptor. We test different charged defect correction schemes, the influence of Ga3*d*-electrons, and compare the results with accurate photoluminescence measurements. We find that using different charged defect correction methods can lead to substantially different hybrid functional parametrizations. However, the calculated optical properties of the carbon acceptor are found to be weakly dependent on specific parameters.

## I. INTRODUCTION

Gallium nitride is a wide bandgap semiconductor with numerous applications in light-emitting devices,^{1} lasers,^{2} solar cells,^{3} and high-power electronics.^{4,5} A wide variety of devices depend on the ability to control and take advantage of the properties of defects in GaN. Point defects in GaN can be detrimental for device operation because they lead to lowering of the breakdown voltage in high-power devices. On the other hand, intentional doping by impurities, such as magnesium, leads to the creation of *p*-type materials, indispensable in solid-state lighting. One of the most common impurities in GaN is carbon. It is often found as a contaminant in GaN samples grown by a variety of methods. Carbon is present in undoped, Fe-doped, and Si-doped GaN grown by metalorganic vapor phase epitaxy (MOVPE), undoped GaN grown by hydride vapor phase epitaxy (HVPE), and undoped GaN grown by molecular beam epitaxy (MBE). Carbon occupying a nitrogen site in the GaN lattice (C_{N}) acts as an acceptor and is energetically favorable in *n*-type and high resistivity GaN samples with Fermi energy above the middle of the bandgap.^{6,7} Early density functional theory (DFT) calculations, which used local density approximation (LDA), predicted a relatively shallow acceptor level of 0.2 eV above the valence band maximum (VBM) for C_{N}.^{6} However, more recent theoretical calculations using hybrid density functional theory consistently yield a deep acceptor 0/– level at 0.9–1.0 eV above the VBM for C_{N}. They also predict C_{N} to be a source of the yellow luminescence (YL or YL1) band with a maximum at about 2.2 eV in GaN.^{7–13}

The attribution of the YL1 band to the isolated carbon acceptor is firmly supported by recent experiments. In particular, Zvanut *et al.*^{14} have obtained *p*-type GaN by doping with carbon. The temperature dependence of the hole concentration revealed the ionization energy of the C_{N} acceptor of about 1 eV. Electron paramagnetic resonance (EPR) experiments also showed the presence of C_{N} acceptors that can be ionized with 1 eV light.^{14} In photoluminescence (PL) studies, the concentration of defects responsible for the YL1 band matches the concentration of carbon estimated from secondary ion mass spectrometry analysis.^{15} Finally, a second transition level (0/+) at about 0.3 eV above the valence band maximum (VBM) for the C_{N} defect, predicted theoretically,^{8,9} has been detected experimentally.^{15,16} The YL1 band is shown in Fig. 1.

The following properties of the YL1 band and C_{N} acceptor have been established.^{15,17} At *T* < 20 K, the YL1 band is caused by transitions from shallow donors to the 0/− level of C_{N}, the so-called donor–acceptor pair (DAP) recombination. At *T* > 50 K, electron transitions from the conduction band to the same acceptor replace the DAP transitions. The zero-phonon line (ZPL) of the YL1 band is found at 2.57 eV at *T* < 20 K and 2.59 eV at *T* ≈ 50 K in unstrained GaN. Analysis of ZPL in a large number of GaN samples reveals that the 0/− level of the C_{N} is located at 916 ± 3 meV above the VBM. The maximum of this broad band is at 2.17 ± 0.05 eV, and the full width at half maximum is about 0.44 eV at low temperature. The hole-capture and electron-capture coefficients for the C_{N} acceptor are (3.7 ± 1.6) × 10^{−7} and (1.1 ± 0.1) × 10^{−13} cm^{3}/s, respectively, and their temperature dependence is insignificant. While YL in GaN is sometimes also attributed to other defects, such as V_{Ga}, V_{Ga}-3H, and V_{Ga}O_{N}-2H,^{18} we do not observe any other YL band except for the YL1 band with the above-listed parameters and properties. The reliable identification and accurately determined parameters of the YL1 band make it an excellent case-study object for comparison with theoretical calculations.

The most widespread theoretical approach for calculating the optical properties of defects in semiconductors is hybrid density functional theory. Among different hybrid functional methods, for calculations of defects in semiconductors, Heyd–Scuseria–Ernzerhof (HSE)^{19,20} hybrid functional is often used, with functional parameters fit to reproduce the host semiconductor bandgap. Recently, it has also been suggested that correct electronic interactions can be reproduced in HSE if the generalized Koopmans' condition for the defect orbitals is fulfilled instead of (or in addition to) tuning HSE to the experimental gap.^{13,21–23} In this approach, the HSE parameters are adjusted to cancel the self-interaction energy by wavefunction relaxation with the addition (or removal) of an electron to (from) the defect state orbital. This procedure restores the correct linear behavior of the total energy *E*(*N*) with electron occupation in HSE. Usually, this is achieved by finding the HSE parameters, fraction of exact exchange *α* and range separation parameter *μ*, for which the electron addition energy is equal to the defect state Kohn–Sham eigenvalue $E(N+1)\u2212E(N)=ei(N)$.^{22} However, within HSE, the set of parameters that fulfills this condition is not unique, and a range of *α* and *μ* can satisfy it for a given defect orbital. Therefore, recently, it was suggested that this allows one to tune HSE to both obtain the correct bandgap and enforce the generalized Koopmans' condition for the defect state.^{13} The tuning procedure involves computing total energies of the defect upon addition/removal of an electron, and in periodic supercells, it depends on the correction scheme for the spurious interaction between periodic images of charged defects. The two most commonly used correction schemes are due to Lany–Zunger (LZ)^{24} and Freysoldt–Neugebauer–Van de Walle (FNV).^{25,26} Both approaches yield similar results, within less than 0.1 eV for ∼100 atom supercells,^{27} but it is not yet clear which method is more appropriate for tuning of the HSE since small total energy differences (less than 0.1 eV) can lead to substantially different values of *α* and *μ* that appear to satisfy the generalized Koopmans' condition. Here, we test HSE tuning to enforce the generalized Koopmans' condition using LZ and FNV charged defect correction schemes. We also test a commonly used parametrization of HSE which only reproduces the experimental bandgap, as well as the inclusion of the semicore Ga3*d*-electrons in the valence band.

## II. METHODS

All calculations were performed using Vienna Ab initio Simulation Package (VASP)^{28} using the projector augmented-wave method. HSE hybrid functional calculations were carried out with a varying fraction of exact exchange *α* and range separation parameter *μ* in 128 atom hexagonal (4 × 4 × 2) supercells at the Γ point, with plane-wave energy cutoffs of 500 eV. During the initial tuning procedure, lattice parameters were kept at values (*a* = 3.21 Å, *c* = 5.20 Å) relaxed for a common gap-tuned HSE (*α* = 0.31, *μ *= 0.2). During initial tuning, unrelaxed defect geometries were used, which allow more accurate calculations of charged defect corrections in the FNV scheme, circumventing complications due to ionic screening effects. Electron addition energies were computed as total energy differences of the negatively charged and neutral C_{N} acceptors. Total energies of the charged defects were corrected using both LZ and FNV correction schemes. In the LZ model, the artificial interactions between charged defects are calculated as the electrostatic energy due to periodic images. It contains the well-known Madelung energy of the periodic array of point charges and the energy of interaction of the localized charges with the compensating background, which also includes screening effects. Here, the LZ scheme was also applied in combination with the potential alignment procedure.^{27} In the FNV model, the defect induced potential is obtained directly from ab initio calculations and is used to match the potential of the model charge density. The correction for the charged defect total energy is then calculated from the Coulomb interaction of the periodic array with model charge density. A detailed review of both correction schemes is given in Ref. 27. In both correction schemes, a high-frequency (ion-clamped) dielectric constant $\epsilon \u221e=5.0$ was used for unrelaxed lattice structures. This value was calculated for GaN using the self-consistent polarization response to the applied electric field. The high-frequency dielectric constant was found to be nearly isotropic, 4.95 in-plane and 5.08 along the wurtzite *c*-axis, and practically independent of HSE parametrization, within the parameter range used in this work.

Figure 2 shows the planar averaged defect induced electrostatic potentials of negatively charged and neutral C_{N} acceptors (solid lines). The defect is located at the origin, and the first periodic image is 4*a* away at 12.8 Å. The electrostatic potential of the negative acceptor behaves as expected, with the model potential for FNV correction method^{25,26} calculated using the above dielectric constant (dashed line) reproducing *ab initio* calculation very well. However, for the neutral C_{N}, the electrostatic potential is not constant between defects, indicating interactions between periodic images. In previous works, interactions between formally neutral defects were analyzed using an effective defect charge.^{29,30} Similarly, here these interactions are estimated by adopting an effective defect charge of *q *= 1 in the FNV model and finding the model potential (dashed line) that approximates the *ab initio* potential (solid line). This estimation suggests that the energy of spurious interactions of neutral defects is roughly 3.6 times lower than that of negatively charged defects. Therefore, during HSE tuning, FNV and LZ corrections were applied to both negatively charged $CN\u2212$ and neutral $CN0$ defects (for the latter, using an effective dielectric constant of 18).

Optical transitions via C_{N} were computed from defect energies following the established procedure.^{31,32} For each parametrization of HSE, lattice constants of bulk GaN were relaxed minimizing forces to 0.01 eV/Å or less. Atoms in supercells containing defects were relaxed to 0.05 eV/Å. The PL maximum is calculated as

where $E(CN0)$ is the total energy of the neutral acceptor, $E\u2217(CN\u2212)$ is the total energy of the negatively charged acceptor in neutral acceptor lattice geometry, and *E _{g}* = 3.5 eV is the experimental bandgap of GaN. ZPL is calculated as

where $E(CN0)$ and $E(CN\u2212)$ are total energies of neutral and negatively charged defects in their respective relaxed lattices. All total energies were corrected applying LZ and FNV electrostatic corrections as outlined above. One important caveat is that high-frequency dielectric constant $\epsilon \u221e$ is used for corrections of $E\u2217(CN\u2212)$ (negatively charged defect in neutral lattice geometry) when calculating the PL maxima, since in the neutral defect lattice, ionic screening effects are negligible. On the other hand, the low frequency (relaxed-ion) dielectric constant $\epsilon 0$ was used for the relaxed negatively charged defect, when calculating ZPLs to account for ionic screening. This is in line with the analysis of defect corrections given in Ref. 27. Our tests show that $\epsilon 0=8$ produces the most reasonable FNV model potentials, albeit it is slightly below the experimental value of 9.5.^{33}

## III. RESULTS AND DISCUSSION

Figure 3 shows the electron addition energies calculated using the LZ correction scheme (solid black lines), FNV correction scheme (solid red lines), and unoccupied single-particle eigenvalues (dashed lines) of the neutral carbon acceptor for a range of HSE parameters *α* and *μ*. Higher values of *α* indicate larger fractions of exact exchange, which increases the value of the bandgap and moves the electron addition energy and eigenvalues with respect to each other. When the generalized Koopmans' condition is fulfilled, the electron addition energy is

where *E*(*N*) and *E*(*N *+ 1) are the total energies of the system with *N* and *N *+ 1 electrons, i.e., the neutral and negatively charged C_{N} defect, respectively; and *e _{i}*(

*N*) is the (unoccupied) defect state eigenvalue in the neutral state. In Fig. 3, the intersection points (filled symbols) show the parameter sets for which the generalized Koopmans' condition is fulfilled. Using local approximations to the DFT, such as LDA, leads to the violation of this condition, which can be expressed as

^{22}

where $\Pi i$ is the positive self-interaction energy of the electron in the orbital *i* fixed at the initial unfilled state, and $\Sigma i$ is the negative wavefunction relaxation energy upon addition of an electron. In Hartree–Fock theory (exact exchange is unscreened), the self-interaction energy is zero, and from expression (4) it follows that due to the negative wavefunction relaxation energy, the eigenvalue *e _{i}*(

*N*) is higher than

*E*. On the other hand, in LDA, where self-interaction energy $\Pi i$ is significant, the situation is reversed, and the eigenvalue

_{add}*e*(

_{i}*N*) is lower than

*E*. Therefore, increasing the fraction of exact exchange from zero to one will move the defect state eigenvalue and

_{add}*E*in opposite directions. Similarly,

_{add}*E*increases and

_{add}*e*(

_{i}*N*) decreases with

*μ*because the effective screening length

*l*of the exact exchange changes as $ls=2/\mu $. At the points where

_{s}*E*is equal to

_{add}*e*(

_{i}*N*), we obtain parameters where self-interaction is canceled by wavefunction relaxation (Fig. 3), and Eq. (3) is restored.

In the case of the FNV correction scheme, the electron addition energies are shifted upward by 0.1 eV, resulting in intersections with the single-particle eigenvalues that are about 0.06 eV higher. This shifts the HSE parameters that satisfy the generalized Koopmans' condition by a significant amount. The two sets of parameters obtained using the two correction schemes are compared in Fig. 4(a). With one parameter fixed (*α* or *μ*), another parameter varies by about 0.05 for the two correction schemes. For instance, depending on the preferred electrostatic correction method, one may conclude that the set either *α* = 0.3, *μ* = 0.17 (using FNV corrections) or *α* = 0.3, *μ* = 0.22 (using LZ corrections) satisfies the generalized Koopmans' condition. The same could be claimed about the uncertainty in the fraction of exact exchange *α*. This uncertainty also leads to significant variations in the computed bandgap.

Figure 4(b) shows the dependence of the bulk bandgap of GaN for the set of parameters that satisfy the generalized Koopmans' condition (using LZ and FNV corrections) for C_{N}. HSE hybrid functionals are commonly tuned to the experimental low-temperature bandgap of 3.50 eV for GaN. However, recently it was suggested that since in experiment the bandgap is renormalized by electron–phonon interactions, and calculations are performed in the static lattice, the correct electron interactions are reproduced for the HSE sets of parameters that yield the non-renormalized bandgap.^{13} The zero-temperature gap renormalization energy was estimated to be 0.17 eV for bulk GaN.^{34} Therefore, the target bandgap value, to tune the HSE functional, is 3.67 eV, shown as a dashed line in Fig. 4(b). Assuming that the renormalization affects the valence and conduction band edges similarly for calculations of optical transitions, Eqs. (1) and (2), HSE computed *E*_{VBM} is adjusted upward by half of the renormalization energy. Figure 4(b) illustrates how sensitive the Koopmans' and bandgap tuned HSE parameters are to the adopted electrostatic correction scheme. Using the LZ correction scheme, the HSE parameters that both satisfy the Koopmans' condition and yield the correct bandgap are *α* = 0.5, *μ* = 0.428. Similarly, using the FNV correction scheme, we obtain a significantly different HSE with parameters *α* = 0.36, *μ* = 0.236. Interestingly, if neutral defect energies are not corrected for the (small) electrostatic interaction, we also obtain quite different HSE parameters (*α* = 0.4, *μ* = 0.297) and (*α* = 0.3, *μ* = 0.126) using LZ and FNV corrections, respectively.

It should be noted that tuning HSE to the generalized Koopmans' condition and the bandgap may still not be optimal. For example, it has been suggested that tuning HSE to the Koopmans' condition and bulk ionization potential would yield better optical transition energies relative to the valence band.^{21} This is because in this case, the VBM should be correct with respect to a physical reference, such as the vacuum level. The problem, however, is that the values of the ionization potential in GaN obtained from experiments are scattered between 6.2 and 7.2 eV (see Ref. 21 and references therein), making reliable tuning of HSE to experimental ionization potential difficult. Fitting HSE to fulfill the generalized Koopmans' condition for a well-localized defect state, such as the deep 0/– level of the gallium vacancy, and the value of GaN ionization potential of 6.7 eV (an average value of available experimental results) yields HSE with parameters *α* = 0.25, *μ* = 0.161.^{21} There is significant uncertainty in whether this parametrization accurately reproduces the bulk ionization potential. Nevertheless, it is illustrative to compare this approach to the parametrizations outlined above.

The results of the calculations of optical properties of C_{N} are shown in Table I. Interestingly, even though different approaches in electrostatic corrections lead to significantly different sets of HSE parameters, the resulting overall calculated energies of optical transitions are relatively close to each other. Tuning HSE to satisfy the generalized Koopmans' condition using the LZ correction scheme produces results slightly closer to experiment, by roughly 0.1 eV, compared to those obtained by tuning HSE using the FNV corrections. Including a small correction for the interactions of neutral defects does not affect the computed ZPL's and PL maxima, while values of *μ* and *α* are significantly different. On the other hand, values of PL maxima are all systematically underestimated by 0.1–0.35 eV. This indicates overestimated computed lattice relaxation energies. Tests for HSE (*α* = 0.25, *μ* = 0.161) using 300 atom supercells yield very similar results (ZPL = 2.58 eV, PL_{MAX} = 2.02 eV), indicating that this overestimation is unlikely due to the supercell size. We also note that a commonly used tuning of *α* = 0.31, *μ* = 0.2 produces similar reasonable results. Tests show that for this parametrization, the electron addition energy (using LZ corrections) is 0.07 eV away from the neutral defect eigenvalue, i.e., it is close to being Koopmans' compliant.

HSE α, μ (Å^{−1})
. | ZPL (eV) . | PL_{MAX} (eV)
. | 0/– level (eV) . | E (eV)
. _{g} | a (Å)
. | c (Å)
. |
---|---|---|---|---|---|---|

0.5, 0.428 (LZ) | 2.52 | 1.89 | 0.98 | 3.66 | 3.193 | 5.181 |

0.36, 0.236 (FNV) | 2.46 | 1.82 | 1.03 | 3.67 | 3.200 | 5.193 |

0.4, 0.297 (LZ)* | 2.53 | 1.91 | 0.97 | 3.68 | 3.197 | 5.188 |

0.3, 0.126 (FNV)* | 2.40 | 1.81 | 1.095 | 3.68 | 3.205 | 5.201 |

0.25, 0.161 (LZ)^{21} | 2.59 | 2.06 | 0.91 | 3.22 | 3.212 | 5.215 |

0.26, 0.1 (Ga-d) (LZ)^{13} | 2.56 | 2.05 | 0.94 | 3.63 | 3.181 | 5.171 |

0.31, 0.2 (LZ) | 2.50 | 1.93 | 1.00 | 3.49 | 3.205 | 5.202 |

Experiment | 2.587 | 2.17 | 0.916 | 3.67 | 3.189 | 5.185 |

HSE α, μ (Å^{−1})
. | ZPL (eV) . | PL_{MAX} (eV)
. | 0/– level (eV) . | E (eV)
. _{g} | a (Å)
. | c (Å)
. |
---|---|---|---|---|---|---|

0.5, 0.428 (LZ) | 2.52 | 1.89 | 0.98 | 3.66 | 3.193 | 5.181 |

0.36, 0.236 (FNV) | 2.46 | 1.82 | 1.03 | 3.67 | 3.200 | 5.193 |

0.4, 0.297 (LZ)* | 2.53 | 1.91 | 0.97 | 3.68 | 3.197 | 5.188 |

0.3, 0.126 (FNV)* | 2.40 | 1.81 | 1.095 | 3.68 | 3.205 | 5.201 |

0.25, 0.161 (LZ)^{21} | 2.59 | 2.06 | 0.91 | 3.22 | 3.212 | 5.215 |

0.26, 0.1 (Ga-d) (LZ)^{13} | 2.56 | 2.05 | 0.94 | 3.63 | 3.181 | 5.171 |

0.31, 0.2 (LZ) | 2.50 | 1.93 | 1.00 | 3.49 | 3.205 | 5.202 |

Experiment | 2.587 | 2.17 | 0.916 | 3.67 | 3.189 | 5.185 |

Here, we also compare HSE parametrization used in Ref. 13 (*α* = 0.26, *μ* = 0.1), obtained by tuning HSE to Koopmans' condition and the bandgap, including semicore Ga3*d*-electrons in the valence band. In their respective relaxed primitive cells, the bandgap computed with this parametrization of HSE including Ga3*d*-electrons is 3.63 eV, while removing them reduces the bandgap to 3.49 eV. The Ga3*d*-electrons push the valence band up, which changes the position of the VBM; however, the changes in the exact exchange energy also increase the bandgap. Including Ga3*d*-electrons does not change the position of the defect state eigenvalue above the VBM. However, it moves *E _{add}* down somewhat, which means that without Ga3

*d*-electrons

*E*is higher than the defect state eigenvalue, indicating an excess of self-interaction energy. In other words, such HSE parametrization without Ga3

_{add}*d*-electrons is slightly closer to the local approximations of DFT.

Better agreement with the experimental results is obtained with HSE parameters from Refs. 13, 15 and 21. This is especially surprising for HSE (*α* = 0.25, *μ* = 0.161), given the fact that this parametrization was obtained without tuning to the bandgap, i.e., it yields *E _{g}* = 3.22 eV, significantly lower than the experimental (non-renormalized) value of 3.67 eV. Nevertheless, optical transitions obtained by this HSE parametrization are very similar to those obtained using HSE tuned with Ga3

*d*-electrons, which is a significantly more computationally demanding calculation. This suggests that reproducing the experimental bandgap is less important than satisfying the generalized Koopmans' condition.

## IV. CONCLUSIONS

Using LZ and FNV electrostatic correction schemes, we obtained several parametrizations of the HSE hybrid functional, which satisfy the generalized Koopmans' condition for the C_{N} acceptor and reproduce the experimental bandgap of GaN. Optical properties of the C_{N} acceptor calculated using these parametrizations of HSE are compared with experimental data: (1) The HSE parameters are sensitive to the electrostatic correction scheme used during tuning. Nevertheless, all Koopmans' tuned HSE parametrizations yield reasonable results. (2) The results that agree better with experiment are obtained using Koopmans' gap-tuned HSE with Ga3*d*-electrons. (3) Good results can also be obtained without Ga3*d*-electrons, by tuning HSE to the Koopmans' condition and ionization potential. (4) The results of the latter HSE parametrization (*α* = 0.25, *μ* = 0.161) for C_{N} are in close agreement with the experiment. This parametrization was also shown to reproduce the experiment for the magnesium acceptor Mg_{Ga} and fulfill the Koopmans' condition for the gallium vacancy V_{Ga} in GaN.^{21} This suggests that Koopmans' tuned HSE parametrization is transferrable and can be used for a variety of defects in GaN. This is also in agreement with earlier studies that suggested that the same set of hybrid functional parameters satisfies the Koopmans' condition for a series of defects.^{23}

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (NSF) (No. DMR-1904861). The calculations were performed at the VCU Center for High Performance Computing.