GaN bandgap bias caused by semi-core treatment in pseudopotentials analyzed by the diffusion Monte Carlo method

The electronic structures of group III nitride semiconductors have been intensively studied due to their potential applications in Light Emitting Diodes (LEDs) and solar-cell batteries.^1–5 For example, the bandgaps of AlN and gallium nitride (GaN) were reported in the 1970s as ∼6.28 eV^6,7 and 3.4 eV,^8–10 respectively. InN^11,12 is plagued by difficulties in crystal growth due to the thermal instability and the large difference between the atomic radii of In and N.¹³ Therefore, the experimental value of the InN bandgap has varied from initial reports of 1.89 eV¹⁴ to ∼0.7 eV.^15–19 

Theoretical estimates of bandgaps are predominantly achieved through ab initio calculations based on density functional theory (DFT),^20,21 which typically produces underestimates^22,23 that are sensitive to the choice of exchange-correlation (XC) functional.^24,25 The incomplete cancellation of the self-interaction is a major reason for this underestimation, which prompted the development of the Self-Interaction Correction (SIC) scheme.^26–28 DFT calculations with XC functionals that more precisely consider the exchange part are also useful for solving this underestimation problem; previous examples include B3LYP,²⁹ HSE,³⁰ TB-mBJ,^31–33 and DFT + U.^34,35 Moreover, the GW method³⁶ is an important framework based on the many-body theory that aims to achieve reliable estimates of the bandgap.³⁷ However, GW evaluations are overestimated when the screened Coulomb potential W and the dielectric constants are both updated in the self-consistent cycle.³⁸ Several studies have attempted to understand the origin of bandgap underestimates for group III nitride semiconductors.^5,39,40 For example, previous research on a GaN semiconductor treated the d-using a local density approximation (LDA) method that accounted for the Coulomb interaction potential (LDA + U),⁴¹ which was further applied to the inner s and p.⁴² An earlier LDA study reported the gapless (metallic) band structure of InN with In-4d treated as valence electrons.⁴³ Subsequent DFT research on InN^44,45 reproduced the insulating band structure using the SIC scheme. The GW method has also been applied to InN,⁴⁶ revealing a gap of 0.69 eV, which is close to the most recent experimental value of 0.7 eV.¹⁹ Some of these values are summarized in Table I.

The diffusion Monte Carlo (DMC) method is an electronic structure method that does not rely on the XC approximation⁵⁰ responsible for bandgap underestimation. Although the justification for using DMC for treating excited states (e.g., tiling theorem) is still debated,^51,52 several previous studies have reported fairly accurate estimates of the bandgap^53–58 when compared to experimental measurements and GW results.⁵⁹ In DMC, the initial shape of a many-body wavefunction is numerically relaxed to allow it to become closer to the exact solution.⁵⁰ The shape refers to its amplitude and nodal positions, the latter of which is known to be difficult to update using the DMC framework due to the negative sign problem.⁵⁰ In most practical applications, an approximated implementation that does not update the nodal positions is used, known as fixed node DMC (FNDMC),⁵⁰ which is employed in this study. The fixed node error due to the approximation is known as “uncontrollable;” hence, any predictions made by FNDMC should carefully consider the limitations of the approximation. In addition to this error, other sources of bias such as time-step bias, the finite size error, and population control bias^60,61 exist in FNDMC, some of which are discussed in Sec. IV.

Regarding ab initio bandgap predictions, as well as the bias due to incomplete self-interaction cancellation described above, another bias is commonly observed in DFT and DMC, which is caused by the choice of pseudopotentials (PP) (i.e., large-core or small-core approximation). This is the topic of this study. In group III nitride semiconductors, d electrons of Ga or In are regarded as semi-core electrons, which may be assigned to either the core or valence regions somewhat arbitrarily.⁶² The large-core treatment (whereby semi-core electrons are assigned to the core) results in an underestimation of cohesion because of the lack of contributions by semi-core electrons to the binding.^47,48 Moreover, shallow cohesion leads to the overestimation of the lattice constants.⁶³ Although the large-core treatment leads to less accurate predictions of cohesion, large-core DFT predictions result in better agreement between bandgap predictions and experimental bandgap values for InAs,⁶⁴ GaN,^47,65 and InN.⁴⁵ Additional research into the dependence of bandgap predictions on core size has involved DFT + U, SIC, and GW methods.^47,48,63,65 Furthermore, the small-core/LDA + U method applied to GaN⁴¹ resulted in a closer agreement between predicted and experimental bandgaps; however, some large-core/SIC studies for InN^44,45 also reproduced reasonable bandgap values (0.58 eV). These previous studies have suggested that possible hybridization⁴⁷ between (N-2s or 2p) and (Ga-3d or In-4d) orbitals could explain the dependence of bandgap estimates on the semi-core treatment. Specifically, two different factors (the self-interaction cancellation problem and the semi-core treatment problem) have gradually been recognized but have not previously been clearly distinguished. Although it is agreed that large-core PP produce larger bandgap predictions,⁶⁶ the large-core approximation does not necessarily produce better bandgap predictions. In DFT, the bandgap predictions are originally underestimated; hence, a larger prediction may appear to exhibit better agreement with experimental values; however, this problem is specific to DFT due to the self-interaction cancellation problem. In methods that are unaffected by the cancellation problem or some other factors leading to underestimation, the small-core approximation might result in reasonable agreement between predicted and experimental values; therefore, the larger bandgap predicted by the large-core approximation will exhibit less agreement with experimental values.

Thus, to investigate the bandgap dependence on semi-core treatment, a method should be employed that does not suffer from the self-interaction cancellation problem. Specifically, this study uses the DMC method,^{24,50,59,67–75} which has rarely been used to investigate the large-core/small-core problem.⁷⁶ This study confirms that a larger bandgap estimate is given by the large-core PP approximation, which is consistent with previous DFT research. However, the small-core PP approximation in DMC produces an overestimated bandgap; hence, the large-core PP approximation produces an even larger overestimate. A comparison of the difference between small-core and large-core overestimates indicates that the bias cannot be explained by pd hybridization,^45,47,65 which is recognized as a possible cause of bias, but is instead due to the lack of screening by semi-core electrons.

We applied FNDMC calculations⁵⁰ to GaN crystals with a wurtzite structure, which are optimized by DFT to provide relatively good agreement with experimental values,⁴⁶ with only slight overestimation of the c-axis lattice constant by 0.3%. All DMC calculations were performed using CASINO.⁷⁷ Guiding functions for the DMC were prepared in the standard Slater–Jastrow form,⁷³ in which the determinant part was composed of Kohn–Sham orbital functions generated by generalized gradient approximation (GGA) within the DFT package QuantumEspresso.⁷⁸ 

Unlike DFT, only the norm-conserving PP⁷⁹ are applicable to the quantum Monte Carlo (QMC) electronic structure calculation as it is based on the many-body wavefunction theory.⁵⁰ To achieve better consistency with the many-body wavefunction theory, the PP are typically generated by accurate atomic calculations within the quantum chemistry framework.^80–84 Several efforts have been made to construct and provide the PP in the form of a library, considering bias due to the Fock exchange,^85,86 the performance comparison between the energy consistency and shape consistency,⁷⁰ the energy level reproducibility at several different electronic configurations,⁸² and the reproducibility of the density matrices.^84,87 We used trail-needs PP for large-core Ga (valence:4s², 4p¹)⁸¹ and opium PP⁸⁸ for small-core Ga (valence:3d¹⁰ 4s² 4p¹). N is described by the trail-needs PP (valence:2s² 2p³).

A supercell size of 4 × 4 × 2 was employed in DMC, which corresponded to 512 (1152) electrons with the large-core (small-core) PP according to available computational resources. Correspondingly, the Kohn–Sham orbitals were generated by the DFT with a k-mesh size of 4 × 4 × 2, which can achieve coincidence within the 0.0046 Ry/cell with the converged value of DFT ground state energy at 8 × 8 × 4 by using a plane-wave cutoff energy, E_cut, of 100 Ry. A larger cutoff energy can achieve coincidence within the 0.014 Ry/cell with the converged value of DFT ground state energy. The generated orbitals were re-expanded using the blip basis set⁸⁹ for DMC. Wurtzite-GaN is a direct transition semiconductor with both a valence bond maximum (VBM) and a conduction band minimum (CBM) at Γ.¹⁰ The gap is evaluated as the optical gap using an N-electron many-body wavefunction to describe an excited state.⁵³ More specifically, the determinant for the excited state is composed of all occupied orbitals, except one orbital at VBM promoted to that at CBM.^59,90,91 Although this promotion method has been used in several previous studies,^53,59,90 the treatment is a sensitive approximation in DMC that does not relax the newly occupied CBM orbital.

The Jastrow part of the guiding functions was optimized using the variational Monte Carlo (VMC) method so that VMC energies were minimized (energy minimization⁹²). A total of 38 variational parameters comprised of the Jastrow functions, including those in the electron–electron part (ee-terms), the electron–nucleus part (eN-terms), and the plane-wave expansion part (p-terms) in CASINO implementation.^77,93 To reduce the time-step bias in DMC calculations, we took the extrapolation scheme toward δt → 0 using the results at δt = 0.01 and 0.0025 a.u.⁻¹, which achieves the best computational efficiency.⁹⁴ Approximately 2560 walker population were employed to equilibrate the sampling distributions so that the total equilibration steps could achieve their diffusion lengths with the Wigner–Seitz radius as a typical characteristic length of the system. Note that fewer or more walkers would affect the result, as the number of walkers affects the population control bias.^60,61 In general, the bias depends on the quality of the trial wave function adopted. In our study, the nodal structure for the excited states is constructed by promoting the VBM orbital without further relaxations applied to the excited states. This implies that the quality of excited-state trial node is worse. Therefore, the population control bias would get more crucial in excited-state calculations for a gap evaluation than in conventional ground-state calculations. The equilibrations were also confirmed by plotting the behavior of sampled local energies (Fig. 1), whereby the statistical samplings were accumulated until the resultant statistical error bars of the bandgaps for small-core and large-core PP approximations did not overlap each other.

Table I presents a comparison between previous theoretical studies and experiments and the results of this study. DFT values, including those obtained in this study (used to generate Slater determinants), are underestimated compared to the experimental result (3.39 eV).¹⁰ Comparing the small-core DFT results between the present generalized gradient approximation (GGA) method and a previous LDA method, we observe that GGA leads to a slightly larger underestimation. This trend is also observed in a previous DFT bandgap study over several insulators.⁹⁵ Exhibiting relatively good agreement with previous DFT results^47,48 at the initial stage (generating the trial nodes), the present DMC method (the imaginary evolution of the many-body wavefunction toward the exact evolution) increases the bandgap beyond the experimental values. This behavior is expected because a similar trend has been reported in several previous DMC studies.^59,91 The bandgap obtained after the large-core treatment in DMC, $Δg(L)$⁠, is 5.28(6) eV, whereas the bandgap obtained after the small-core treatment in DMC, $Δg(S)$⁠, is 4.1(5) eV. The difference in the predictions, $δΔg(DMC)$⁠, between small-core, $Δg(S)$⁠, and large-core treatments, $Δg(L)$⁠, is ∼1.2 eV for DMC, which is higher than the difference for DFT (0.61 eV). The DMC prediction corresponds to the optical gap;⁵³ therefore, the exciton correction should be taken into account^91,96 when we consider the quasi-particle gap. The correction is estimated as the binding energy of excitons, which is observed at ∼0.023 eV⁹⁷ and is smaller than the statistical errors in the present DMC. However, other DMC research⁹¹ has reported very small differences (∼0.01 eV) between the estimates of optical and quasiparticle gaps. Hence, the exciton correction is not considered in the following discussion.

As mentioned previously, pd hybridization is a possible reason for the observed trend, whereby $Δg(S)<Δg(L)$⁠.^66,98 Hybridization occurs between the p-orbitals at VBM from the d-levels of anions and cations, which pushes the VBM edge upward and makes the bandgap narrower than that without hybridization. As hybridization cannot be described by the large-core treatment (lack of Ga-3d), the bandgap is overestimated.⁶⁶ To determine the extent of hybridization, we depict the partial density of states (pDOS) for large-core and small-core calculations in Fig. 2. Calculation conditions for pDOS are the same as that for generating trial wave functions of DMC. Ga-3d and N-2p are not equal in any energy range for small-core calculations [Fig. 2(a)], implying minimal hybridization. Moreover, the shape of the DOS around the Fermi level is relatively similar in small-core and large-core calculations [Figs. 2(a) and 2(b)], which confirms this finding. Furthermore, hybridization is reported to decrease the bandgap by ∼0.1 eV–0.3 eV,⁹⁸ which does not account for the observed difference, i.e., $δΔg(DFT)=Δg(L_DFT)−Δg(S_DFT)$ = 0.61 eV, where $Δg(L_DFT)$ and $Δg(S_DFT)$ represent large-core and small-core DFT bandgaps, respectively. The interaction between semi-core and valence electrons can be partly explained by core polarization.⁹⁹ However, core polarization results in a minimal change in the bandwidth,⁹⁰ which is likely also the case for the bandgap.

Although the above factors (pd hybridization and core polarization) would be more emphasized by DMC due to its vividly captured electronic interactions, the observed $δΔg(DMC)$ value of 1.2 eV is too large to be explained by these mechanisms alone, implying that another factor is responsible for the observed bias. In a previous GW study,³⁸ the bandgap was overestimated (underestimated) when the electronic screening captured by random phase approximation (RPA) was underestimated (overestimated). As neglecting the semi-core electrons leads to underestimation of screening because of the lack of contributions by semi-core electrons, the large-core treatment is expected to result in the overestimation of the bandgap, which is consistent with the observed trend of $Δg(S)<Δg(L)$⁠. Although the neglected semi-core electrons partially contribute to the underestimation of the screening, the GW estimation could change by ∼0.5 eV.³⁸ 

The finding that $δΔg(DFT)<δΔg(DMC)$ is difficult to plausibly explain because the DMC result is composed of two sources obtained from the subtraction of two quantities, i.e., E_g (ground-state energy) and E_e (excited state energy). The values are $Eg(S)$ = −142.280 7(4), $Ee(S)$ = −142.276 0(4), $Eg(L)$ = −24.306 48(4), and $Ee(L)$ = −24.300 42(6) in hartree, which do not exhibit any particular feature that can be used to identify the dominant source of the observed bias. Although the absolute value of δΔ_g appears to increase from DFT to DMC, it is unchanged when represented in terms of the relative differences, $δΔg/Δg(S)$⁠, with values of 27.9% for DFT and 29.3% for DMC, respectively.

As mentioned previously, there are several uncontrollable biases caused by the fixed-node approximation and the finite-size error in QMC. A full clarification of these biases is not feasible in this study because of limited computational resources; however, the authors refer to the following points. As the FNDMC results essentially depend on the choice of fixed nodal structure,^100–102 we should consider whether the key finding of this study, $Δg(S)<Δg(L)$⁠, would differ with a different choice. The ratio between the exchange and correlation of XC in DFT used to generate the trial nodes also affects the FNDMC results, as observed when using the DFT + U method^101,103 and hybrid XC functionals.^100,102 Although it could be reasonable to adjust the nodes according to this ratio (e.g., by using U) so that the ground state energy would be further stabilized,^100,102,103 the adjusted value could not justifiably be used to evaluate the bandgap for excited states.¹⁰¹ Moreover, for NiO, the bandgap predicted by FNDMC is overestimated when the +U value is applied, which is optimized to stabilize the ground state energy.¹⁰¹ The +U for NiO adjusts the nodal surface so that the FNDMC bandgap estimate is increased by ∼1 eV.¹⁰¹ Thus, our bandgap estimates could be affected in the same way by introducing U. Although it is interesting to investigate whether potential modifications of the results could alter the observed trend of $Δg(S)<Δg(L)$⁠, such evaluations are rather costly because the statistical errors should be converged during bandgap estimation to ensure that the difference in the energy difference, δΔ_g, is distinguished. Nodal structures may also depend on the choice of PP. Note that the nodes mentioned here are those of many-body wave functions. The nodal structures of different PP core sizes are, therefore, not comparable as they have different dimensions (3N dim.). A different choice of PP with a fixed core size might modify the nodal topology within a fixed dimension of the space. However, this topic is not well investigated and not expected to be a critical factor because the variety of PP (within the norm-conserving form, which is only applicable to QMC) does not exhibit as diverse a design strategy as that in XC.

Regarding the finite size error, several schemes have been established to correct this bias.^74,104 The analysis of the performance of such schemes applied to excited states has revealed a fairly significant bias in bandgap predictions.^91,102 The finite size error bias would be critical to the final result if the convergence of the bias according to the system size differed significantly between the large-core and small-core PP results. In this case, the trend of $Δg(S)<Δg(L)$ would be altered by the choice of simulation size. Recent advances have enabled application of the twisting average scheme to bandgap evaluation.¹⁰⁵ Averaging is the optimal method for analyzing the bias; however, it is computationally expensive. Therefore, we instead evaluated the bias using the Kwee–Zhang–Krakauer (KZK) scheme within the DFT evaluation.^49,74 The scheme is applicable to isotropic systems,¹⁰⁴ so is somewhat valid even for hexagonal crystals with no specific long axis. According to the results in Table I, we approximately estimated the bias using the energy difference, (KZK)−(LDA), which only increased the bandgap by several meV for both small-core and large-core approximations. This amount is too small to affect the trend of $Δg(S)<Δg(L)$⁠. An analysis using the twisting average¹⁰⁵ resulted in a bias for Si crystals of ∼0.01 eV when evaluated by a finite size of 3 × 3 × 3. Adding to the Madelung correction, the finite size correction is 0.22 eV, which is not negligible; however, the $Δg(S)<Δg(L)$ trend is maintained with differences on the order of ∼1 eV.

This study evaluated the bias in bandgap prediction according to the choice of valence range in pseudopotentials, taking GaN with Ga-d semi-core electrons as an example. To avoid the DFT bandgap problem (i.e., underestimation due to a lack of self-interaction cancellation), we evaluated the bias using DMC methods and compared it with the trend evaluated by DFT. Although the DMC small-core prediction was slightly overestimated when compared with the experimental value, this trend was consistent with previous DMC research on bandgaps. When the semi-core contributions were omitted by using large-core pseudopotentials, the bandgap was overestimated by ∼30% in both DFT and DMC. Therefore, we considered the possible causes of this overestimation, including pd-hybridization, core-polarization, and the screening effects between semi-core d and valence electrons. Based on previous reports, hybridization and polarization cannot explain the magnitude of the observed bias. However, the trend of over-/underestimation of the bandgap depending on over-/under-screening reported by previous GW research seemed to agree with the fact that the large-core approximation overestimated the bandgap. Nevertheless, the finding that the screening is the cause of the observed bias requires further confirmation because the contribution from semi-core d electrons to all electronic screening has not been thoroughly determined.

All computations were performed using the facilities of the Research Center for Advanced Computing Infrastructure at the Japan Advanced Institute of Science and Technology (JAIST). K.H. acknowledges the financial support from the HPCI System Research Project (Project IDs hp190169 and hp200040) and MEXT-KAKENHI (Grant Nos. JP16H06439, JP17K17762, JP19K05029, and JP19H05169). R.M. acknowledges the financial support from MEXT-KAKENHI (Grant Nos. JP19H04692 and JP16KK0097), FLAGSHIP2020 (Project Nos. hp190169 and hp190167 at K-computer), the Air Force Office of Scientific Research (Grant Nos. AFOSR-AOARD/FA2386-17-1-4049 and FA2386-19-1-4015), and JSPS Bilateral Joint Research Projects (with DST, India).

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