First-principles study of bandgap bowing in BGaN alloys

Ternary alloys of the group-III nitrides have been instrumental in the advancement of electronic and optoelectronic devices. Alloys of wurtzite-phase InN, GaN, and AlN enable light-emitting¹ and laser diodes² with wavelengths ranging from red to ultraviolet. Nitride alloys are also used in high-power and high-frequency electronics.^3,4

Boron-containing III-nitride alloys are now being considered as a potential addition to the nitride family,^5,6 with the goal of engineering BGaN-based heterostructures that would exploit differences in bandgap or polarization. Predicting the properties of boron-containing alloys is challenging because the ground state of BN is in the hexagonal phase, while the conventional III-nitrides (AlN, GaN, and InN) are stable in the wurtzite phase; the addition of small concentrations of boron is expected to result in wurtzite-phase structures. In addition, the small size of boron (relative to Al, Ga, and In) is expected to have a major impact on the atomic and electronic structure. The size mismatch is expected to lead to a large miscibility gap when boron is alloyed to the wurtzite III-nitrides.⁷ Nevertheless, a number of experimental growth efforts have demonstrated single-crystal wurtzite BGaN alloys,^6,8–10 albeit at low B concentrations.

Little or no information is available about the properties of BN in the wurtzite phase (wz-BN). First-principles calculations for wz-BN have predicted a large indirect bandgap and substantial spontaneous polarization.¹¹ In this paper, we present a comprehensive computational investigation of the electronic structure of BGaN alloys as a function of B content.

Prior density functional calculations of boron incorporation in GaN^12,13 lack predictive power because of the use of functionals, such as the local density approximation¹⁴ or the generalized gradient approximation (GGA)¹⁵ that significantly underestimate the bandgap and do not allow for a quantitative description of the alloy electronic structure. Akiyama et al.¹⁶ attempted to overcome this issue through the use of a single-shot hybrid functional calculation after atomic positions were relaxed within the GGA. However, they did not have a methodology for disentangling the supercell band structures and also neglected the role of disorder by considering only the energetically most favorable configuration for each concentration.

In this paper, we employ density functional theory with a hybrid functional to accurately describe the atomic and electronic structure of BGaN alloys. The bandgap of wz-BN is indirect, while the bandgap of wz-GaN is direct. A transition from a direct to an indirect bandgap is, therefore, expected as the boron content is increased. The alloys are described using large supercells, and the band structure in those supercells exhibits zone folding; simple inspection thus does not allow distinguishing between direct and indirect bandgaps. We, therefore, utilize a projection scheme for identifying extended states associated with the direct and indirect bandgaps. In addition, we compute dipole matrix elements and verify that the states identified as direct by the projection scheme lead to strong optical transitions.

Our calculations are based on density functional theory as implemented in the VASP code^17,18 with the hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE).¹⁹ These calculations use projector augmented wave pseudopotentials²⁰ with a plane-wave energy cutoff of 420 eV; the Ga $d$ states are included in the core. The HSE mixing parameter $α$ is set to 0.287 and is chosen to produce the best agreement in the structural and electronic properties with the available experimental data. This chosen mixing parameter is smaller than previously used values^21,22 due to a change in the Ga and N pseudopotentials with the release of VASP 5.2.2. For the 4-atom wurtzite primitive cells, the Brillouin zone is sampled with a $7×7×4$  $Γ$-centered Monkhorst-Pack grid; for the 72-atom supercells used for alloy calculations, a $2×2×2$  $Γ$-centered Monkhorst-Pack grid is employed. Atomic relaxation is performed until forces are below 5 meV/Å.

To simulate BGaN alloys, we construct a $3×3×2$ supercell of the 4-atom primitive cell, giving rise to 36 cation sites to populate with Ga or B atoms. For each B concentration, 10 distinct atomic configurations are generated by random substitution of the B atoms on the cation sites, thus mimicking disordered alloys.⁸ We obtained data for B concentrations up to 50%; while we recognize that this value is well above the miscibility limit, the results provide interesting insights into the trends and underlying physics. The resulting properties are averaged over the configurations with equal weights. The lattice parameters of each supercell are varied according to Vegard’s law,

where $x$ is the fractional B concentration and $a i Y$ is the $i$th lattice parameter of crystal $Y$⁠. We verified adherence to Vegard’s law for a subset of alloy supercells by allowing for a full relaxation of the lattice parameters and finding a linear variation with B concentration. Input file generation and subsequent analysis was facilitated by the pymatgen library.²³ 

Since one parent material has a direct bandgap and the other an indirect bandgap, we need to track the band-edge eigenstates as a function of composition in order to identify the transition from a direct to an indirect bandgap. Within the alloy supercells, band folding occurs due to the reduced size of the Brillouin zone: eigenstates within the Brillouin zone of the primitive cell are mapped back into the supercell Brillouin zone by a reciprocal space lattice vector, as depicted in Fig. 1(a). To disentangle the electronic structure, the projection scheme of Shen et al. is utilized.²⁴ In this scheme, eigenstates are identified by their projection weight with a band-edge eigenstate from a reference calculation of a parent material strained to the lattice parameters of the alloy. The projection weight is

where $|α⟩$ is a band-edge eigenstate extracted from GaN or BN and $| SC , i ⟩$ is an eigenstate with band index $i$ from the alloy supercell calculation. We identify the supercell eigenstates $| SC , i ⟩$ with the largest projection weights on the direct or indirect band-edge eigenstates $|α⟩$⁠; their eigenvalues are then used to evaluate the direct and indirect bandgaps.

Once the direct and indirect bandgaps have been identified as a function of B concentration, it is customary to fit the resulting bandgap values to a quadratic equation,

where $E g BN$ is the bandgap of wz-BN and $E g GaN$ is the bandgap of wz-GaN. The bowing parameter $b$ describes the quadratic deviation from a linear interpolation of the parent materials’ bandgap.

While the projection scheme allows us to identify the bandgap of each alloy supercell, it does not provide any information about the strength of the optical transitions. Optical absorption measurements are commonly used to experimentally determine the bandgap of a material. The strength of transitions is governed by the dipole matrix element. We compute the dipole matrix elements for all possible direct valence- to conduction-band transitions, allowing us to compare the transitions with a strong dipole matrix element with the direct bandgaps selected by the projection scheme.

Lattice parameters for wz-GaN are found to be $a=3.19Å$ and $c=5.18Å$⁠, in good agreement with experimental values of $a=3.19Å$ and $c=5.19Å$⁠.²⁵ For wz-BN, our calculations predict $a=2.53Å$ and $c=4.18Å$⁠. Nagakubo et al. have shown that hexagonal BN can be transformed into the wurtzite phase under high pressure and temperature and measured the lattice parameters to be $a=2.55Å$ and $c=4.22Å$⁠,²⁶ within 1% of our calculated values.

The band structures of the parent materials are shown in Figs. 1(b) and 1(c). The calculated direct bandgap of wz-GaN is 3.48 eV, close to the experimental value of 3.51 eV.²⁷ In Fig. 1(b), we label the valence-band maximum (VBM) of wz-GaN as $Γ v$ and the conduction-band minimum (CBM) as $Γ c$⁠. The indirect bandgap (⁠ $Γ v → K c$⁠) of wz-GaN is 6.60 eV. For wz-BN, our calculations predict an indirect bandgap of 6.70 eV, slightly differing from previous theoretical predictions due to our choice of the HSE mixing parameter.^11,24 In Fig. 1(c), the VBM of wz-BN is labeled by $Γ v$⁠, and the CBM is labeled by K $c$⁠. The eigenstate in wz-BN that most resembles the CBM of wz-GaN is labeled by $Γ c$⁠, and the corresponding $Γ v → Γ c$ gap is 13.20 eV. Other conduction-band states at $Γ$ that could lead to a direct bandgap would be relevant only at B concentrations exceeding 50%, outside the range examined in the present study. The states labeled $Γ v$⁠, $Γ c$⁠, and K $c$ will be used in the projection scheme to identify the direct and indirect bandgaps in the alloy supercell.

When computing the projection weights, the band-edge eigenstates are taken from the bulk wz-GaN calculations. However, our results indicate that there is no difference in the selected alloy eigenstates when wz-BN states are used as the reference. This is due to the fact that the spatial distribution of the reference eigenstates in wz-GaN and wz-BN is nearly identical once the pristine cell has been strained to the lattice parameters of the alloy composition.

In the alloy supercell, we define the direct bandgap as the difference in the energy eigenvalues of the states with the largest projection weight onto $Γ v$ and $Γ c$⁠, and similarly, the indirect bandgap is defined based on the states with the largest projection weight onto $Γ v$ and K $c$⁠. The extracted bandgaps as a function of alloy composition are shown in Fig. 2 and a fit to Eq. (3) is included for both bandgaps.

The indirect gaps of GaN and BN are very similar and the alloy bandgap exhibits only small deviations from linear interpolation. The bowing parameter for the indirect bandgap is found to be 0.55 eV. Deviations from linearity are much larger for the direct gap, and it also turns out that a quadratic fit with a single bowing parameter is inadequate to describe the direct bandgap over the entire range of B concentrations. To quantify the quality of the quadratic fit, we compute the mean squared error (MSE)

where $n$ is the number of data points, $Y i$ is the observed value, and $Y i ^$ is the fitted value. The MSE for the quadratic fit up to a given B concentration is shown in Fig. 3. Up to 20% B, the MSE for the direct gap is negligibly small; this indicates that if we only care about boron concentrations up to 20%, a very good fit can be obtained by including only data points up to this concentration in the fit (see Fig. 2); the corresponding bowing parameter is 8.68 eV. For comparison, the bowing parameter obtained from fitting all direct-bandgap data is 10.65 eV.

Experimentally, the bandgap of BGaN has been measured using optical techniques.^9,10 However, in the analysis of the experimental data to extract a bowing parameter by fitting to Eq. (3), the indirect gap of wz-BN was used as the value for $E g BN$⁠. Instead, the value of the direct gap should be used. A reanalysis of the data in Ref. 10, using the value of the direct bandgap of wz-BN determined in our study (13.20 eV), results in a bowing parameter of 11.48 eV. This value is in reasonable agreement with our calculated value. The remaining discrepancy may be attributed to alloy fluctuations (apparent from the increase in the full-width half maximum of the experimental data). The more pronounced decrease in bandgap observed in Ref. 9 may be due to larger fluctuations.

A simple linear interpolation of bandgaps would predict that the crossover from a direct to an indirect bandgap would occur at a B concentration of 32%. Because of the large bowing of the direct gap, the crossover actually occurs at a concentration well above 50% (see Fig. 2). Another consequence of the large bowing is that the bandgap does not increase much when small amounts of boron are added to GaN. For instance, at a boron concentration of 3% (which is the highest concentration at which high-quality random alloys could be grown by MBE⁸), a linear interpolation of the gaps would indicate that the alloy bandgap is 0.30 eV higher than the GaN gap; but when bowing is taken into account, the actual increase in the gap is only 0.05 eV. We note that these conclusions apply to strain-free alloys. Strain in epitaxially grown films may lead to additional changes in the gap. The boron incorporation may also be useful for strain engineering.

The standard deviations in Fig. 2 reflect the spread in bandgap values obtained from calculations with different atomic configurations; the deviations are larger for the indirect bandgap than for the direct gap. This difference can be attributed to the fact that the eigenstate labeled by $K c$ has mostly $p$-orbital character, as opposed to the $s$-orbital character of the eigenstate $Γ c$ [see Figs. 1(b) and 1(c)]. $p$-like states are less symmetric than $s$-like states and, therefore, are more sensitive to the atomic relaxations, explaining why the K $c$ energy eigenvalue displays a larger variance.

In Fig. 4, finally, we show the magnitude of the computed dipole matrix elements, with the quadratic fit of the direct-bandgap data from the projection scheme (for B concentrations less than 20%) overlaid for comparison. Very good agreement between the two methods is found for B concentrations up to 20%. The agreement indicates that the direct bandgap identified by the projection scheme is an allowed transition in the alloy and would be identified as the bandgap in an optical experiment. The strong matrix elements parallel to the $c$ axis [Fig. 4(b)] systematically occur at slightly higher energy transitions than those perpendicular to $c$ [Fig. 4(a)]; this difference can be attributed to the crystal-field splitting of the VBM, where the $p z$ orbital is lower in energy than the $p x$ and $p y$ orbitals.

Using first-principles computational techniques, we elucidated the electronic structure of BGaN alloys. We used a projection scheme to identify the direct and indirect conduction-band minima; a transition from direct to indirect is predicted to occur at B concentrations above 50%. A bowing parameter of 0.55 eV was found for the indirect bandgap. The direct bandgap cannot be described by a quadratic fit, but for B concentrations less than 20%, a bowing parameter of 8.68 eV is sufficiently accurate. We also calculated dipole matrix elements which confirmed that the direct-bandgap states selected by the projection scheme correspond to allowed optical transitions. The large bowing of the direct bandgap means that for low B concentrations, the bandgap increases only very slowly with the B content. The incorporation of small amounts of B into wz-GaN is, therefore, not an effective means of increasing the bandgap for strain-free alloys but may be useful for strain engineering.

The research reported here was primarily supported by the National Science Foundation (NSF) through the Materials Research Science and Engineering Center at UC Santa Barbara (No. DMR-1720256) (Seed). D.W. was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DE-SC0010689. Computational resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF (No. ACI-1548562).