Bandgap narrowing in zincblende III–V semiconductors: Finite-temperature full random-phase approximation and general analytical model

The bandgap narrowing (BGN) in zincblende III–V semiconductors is calculated in a finite-temperature full Random-Phase Approximation (RPA) formalism based on an isotropic dispersion model. The cases of n-type and p-type quasi-neutral regions and the case of a neutral electron–hole plasma are elaborated for the technologically important materials GaAs, AlAs, InAs, GaP, InP, GaSb, InSb, zb-GaN, zb-InN, Al 0.3 Ga 0.7 As GaAs 0.5 Sb 0.5 , InP 0.69 Sb 0.31 , InAs 0.4 P 0.6 , InAs 0.4 Sb 0.6 , In 0.52 Al 0.48 As In 0.49 Ga 0.51 P, In 0.53 Ga 0.47 As In 0.5 Ga 0.5 Sb, and zb-Ga 0.5 In 0.5 N (60 cases). In quasi-neutral regions, the correlation energy of the interaction between carriers and ionized dopants adds two terms to the total BGN. At low temperatures, inefficient screening makes the hole term dominant in n-type materials with a large ratio of the valence band to the conduction band (CB) density-of-states. The inclusion of the CB nonparabolicity is decisive here, as it prevents a diverging BGN at high concentrations. For all 60 cases, the BGN is evaluated in the temperature range from 0 to 500 K. A strong temperature dependence over the whole density range is observed


I. INTRODUCTION
Bandgap narrowing (BGN) in semiconductor devices is a many-body effect that occurs in heavily doped regions and under strong optical or electrical excitation. In the first case, as in bipolar transistors (BTs) and field-effect transistors (FETs), the activated doping is compensated by a one-component plasma either in the emitter of a BT or in the source/drain of a FET. In the second case, a two-component plasma is generated, as in optically and electrically pumped laser diodes. Device characteristics can be altered by BGN, examples being the reduced current gain in a BT due to the increased emitter minority concentration or the density dependence of the laser wavelength. Modeling of BGN is facilitated by the condition of charge neutrality, which holds in quasi-neutral regions (negligible field strength) and in a neutral plasma (background doping is small compared to the plasma density). The physical origin of BGN is the many-body Coulomb interaction between free charge carriers as well as between charge carriers and ionized dopants. The resulting total self-energy shifts the band edges such that the bandgap shrinks. Four contributions have to be considered in a neutral electron-hole (eh) plasma-the exchange and correlation energies of electrons and holes, whereas five components contribute in quasineutral regions-the exchange energy of the majority carriers, the free-carrier correlation energies of the majority and minority carriers, and the ionic correlation energies of the majority and minority carriers. In this paper, all self-energies are calculated in the finitetemperature full Random-Phase Approximation (RPA). 1-6 RPA is increasingly accurate at higher densities where BGN becomes significant for device operation, and it is the only method to cover the whole temperature range. In n-type III-V zincblende semiconductors with small electron effective masses, the temperature dependence becomes extreme even at high densities. This is because the ionic correlation energy of holes is the largest contributor to the total BGN here. 7 The shift of the minority carrier band is a consequence of the lowering of the energy of the generated/annihilated minority carrier due to Coulomb attraction by the mobile majority carriers and Coulomb repulsion by the sub-system of immobile ionized dopants. 8 The basic RPA theory as well as results for the eh-plasma in various materials can be found in Ref. 1. The extrinsic semiconductor with a random distribution of dopants was treated in Ref. 4 for T = 0 K. Earlier calculations of the self-energy of the carrier-dopant interaction in the T = 0 limit were performed with the Hartree-Fock variational method for donors distributed on a regular sub-lattice 9 applying the single plasmon pole (SPP) approximation 10 and with second-order perturbation theory for randomly distributed dopants based on RPA screening. 11 The multi-valley case had been discussed for both dopant arrangements in Refs. 12-14. Measurements of BGN are accompanied by an intricate analysis. Three types of experiments have been developed and mostly applied to silicon due to its outstanding role in technology: Absorption experiments, [15][16][17][18][19][20] photoluminescence (PL) experiments, [21][22][23][24] and electrical measurements of BTs. [25][26][27][28][29][30] In absorption experiments, phonon-assisted transitions are present besides the band-band transitions; therefore, the phonon energies must be known as well as the exact position of the Fermi energy. The evaluation of PL experiments is difficult because of the actual initial states and the weak intensity of the spectrum. The electrical method requires the measurement of collector current in the BT as a function of the emitter-base voltage, sheet resistance below the emitter, and minority carrier mobility in the base. An analytical model is then used to calculate the BGN from these quantities. BGN data for III-V materials are rare; for most of the 20 materials considered in this paper, they do not exist. Olego and Cardona 31 performed PL measurements of the BGN in p-type GaAs at 2.1 K for concentrations of (1.6, 4, 9) × 10 19 cm −3 . The same material was analyzed by Titkov et al. 32 by PL spectroscopy at 4.2 K in the density range 3 × 10 18 -4 × 10 19 cm −3 . Tiwari and Wright 33 used hetero-structure BTs to measure the BGN in p-type GaAs at room temperature in the density range (10 18 -2 × 10 19 ) cm −3 . Semikolenova et al. performed absorption measurements of the BGN of n-type InAs at room temperature for tellurium concentrations between 2 × 10 17 and 3 × 10 18 cm −3 . PL spectra of p-type GaSb with densities between 3 × 10 18 and 4 × 10 19 cm −3 at 4.2 K were published by Titkov et al. 32 The BGN in n-type InP at 300 K was studied by absorption and PL experiments by Bugajski and Lewandowski, 35 covering the concentration range 10 16 -10 19 cm −3 . In their careful analysis, they considered the nonparabolicity of the conduction band (CB), band filling, and band tailing.
In view of the sparse experimental data and the complications in extracting the actual BGN, reliable calculations and analytical models derived from them are an important means to provide the electronics and opto-electronics communities with information on the density-and temperature-dependent BGN in zincblende III-V materials.
The paper is organized as follows. Sec. II gives a short outline of the RPA theory of BGN and of the band structure model that allows for both the numerical computation and the derivation of a general analytical model. All material and band-structure related parameters of the 20 zincblende III-V semiconductors are listed. The importance of nonparabolicity is demonstrated for the three materials with the smallest gap (and smallest electron effective mass). It is shown how the hole masses are calculated from Luttinger parameters and how the limitations of the band structure model lead to the allowed density maxima for each material. Available experimental upper bounds of the densities from measured limits of the activated doping are cited and found to always be smaller than the theoretical limits. Section III provides a summary of the results of the numerically calculated BGN as a function of material at various densities and temperatures. The trends in comparison to bandgap and effective masses and the differences between neutral eh-plasma, quasi-neutral n-type, and quasi-neutral p-type materials are discussed. The strong temperature dependence of BGN in n-type materials with small electron effective masses (like InSb) is contrasted with the weak one in n-type materials with large electron effective masses (like AlAs). Section IV presents the formulas of the general analytical model of BGN in quasi-neutral regions and in the neutral eh-plasma as functions of density and temperature, avoiding any material-dependent free fit parameter. The way the multiple integrals over complicated functions are solved is sketched, and the origin of the strong temperature dependence of the BGN in n-type materials is explored. It is shown that nonparabolicity increases the electron-ion correlation energy but decreases the hole-ion correlation energy. It is demonstrated how the latter effect weakens the BGN in n-type materials with small electron effective masses, preventing the bandgap from shrinking to zero in the relevant density range. Conclusions are given in Sec. V. Appendix A contains the complete presentation of the full-RPA BGN in all 20 materials as a function of density in the temperature range from 0 to 500 K. Each figure comprises quasi-neutral n-type, quasi-neutral p-type, and neutral eh-plasma and also displays the corresponding analytical curves. The density limits can be inferred from the endpoint of the density axis. Appendix B provides the derivation of the analytical expressions of the free-carrier correlation energies in quasi-neutral regions and in the neutral eh-plasma, whereas Appendix C presents the derivation of the analytical formulas of the ion-carrier correlation energies in quasi-neutral regions including the effect of nonparabolicity in n-type material, and the temperature dependence of the donor-hole correlation energy including the effect of nonparabolicity on this temperature dependence.

II. RPA THEORY AND BAND STRUCTURE MODEL
The full-RPA BGN is calculated by the quasi-particle shift (QPS) 1 Δa(k) (a = e for electrons, a = h for holes), which is given by the difference between interacting and free dispersion, 82 and equals the real part of the self-energy Σa, It consists of the three terms the unscreened exchange energy Δ x a (k), the correlation energy of free carriers Δ c a (k), and the correlation energy of the interaction between the carriers and ionized dopants Δ i a (k). A random arrangement on regular lattice sites is assumed for the dopants. Since the REVIEW pubs.aip.org/aip/adv dispersion of the QPS in the energy interval between band edge and Fermi energy is rather flat for the relevant densities, 1,41 the dispersive QPS Δa(k) can be replaced by a rigid shift Δa. Then a selfconsistent solution of the problem is not necessary because the rigid shift drops out in the energy difference Ea(k + q) − Ea(k), which becomes E 0 a (k + q) − E 0 a (k), i.e., fully determined by the free dispersion. In the distribution functions, Δa is fixed by the given density in the quasi-neutral region or by the plasma density in the case of the neutral eh-plasma. The condition that the QPS density should not change in first order with respect to Δa(k) − Δa results in 1 , with rigidly shifted bands Ea(k) = E 0 a (k) + Δa and Fermi-Dirac functions fa(k) depending on shifted chemical potentials μ a , The first derivative of the Fermi-Dirac function in (4) plays the role of a weight that filters out energies near the Fermi energy in the low-T/high-density limit and energies close to the band edge in the high-T/low-density limit. The dispersive QPS (3) In the correlation energies (7) and (8), v(q) = e 2 /(ϵ 0 ϵsq 2 ) is the Fourier transform of the bare Coulomb potential with the static permittivity ϵs. Ων = 2πi ν ̵ hβ denotes the Matsubara frequency (ν integer) and n i the concentration of electrically active dopants. The correlation energies are governed by the RPA dielectric function, where the energy difference Ea(k + q) − Ea(k) in the denominator is replaced by E 0 a (k + q) − E 0 a (k) (rigid shift). In the parabolic band model, the QP density na in the rigid shift approximation is defined by with the Fermi integral F 1/2 , the thermal wavelength Λa = (2π ̵ h 2 β/ma) 1/2 , the inverse thermal energy β = 1/k B T, and the Fermi energy ζ a = μ a − Δa. The simple form of the exchange energy (6)-just a single integral (the Fermi integral)-is bound to parabolic bands E 0 a (k) = ̵ h 2 k 2 /2ma. 1,42 In contrast, the rigid-shift expressions of the correlation energies (7) and (8) are obtained by employing (4) regardless of the specific form of the band dispersion. The exchange-correlation energy of the electron and hole plasmas is rather insensitive to band structure details as a consequence of a compensation effect. 1,43 Therefore, parabolic bands are used to evaluate Eqs. (6) and (7) in this paper. However, the nonparabolicity of the CB is crucial for the ionic correlation energy Δ i h,n−type of certain n-type III-V materials in the low-T/high-density regime 7 since, due to the small n-density-of-states (DOS), the Fermi level moves deep into the CB. The electron density ne and the dielectric function ϵ(q, Ων) in (8) are then calculated with the nonparabolic free dispersion, where γ e denotes the nonparabolicity parameter (the common symbol α is used for the mass ratios α {e,h} in this paper). In addition, for the shift of the majority band edge in n-type material, Δ i e,n−type , and for the shift of the minority band edge in p-type material, Δ i e,p−type , the nonparabolic dispersion (11) is applied for electrons. Figure 1 demonstrates the nonparabolicity effect on Δ i h,n−type , which dominates the total BGN in all n-type materials with me ≪ m h (α h ≪ 1). It is particularly strong for the small-gap materials with strong CB-valence band (VB) coupling (Nos. 1-5 in Table I) but is also important for all other direct materials. The only exception is the indirect materials AlAs and GaP, where me ≈ m h (α {e,h} ≈ 0.5).
As the evaluation of (7) comes down to a sevenfold integration and that of (8) to a sixfold integration, an isotropic dispersion model is inevitable. A full numerical treatment based on a more elaborate band structure model is out of reach even in the rigidshift approximation of BGN. Therefore, the real DOS has to be approximated by an isotropic DOS, and its range of validity must be carefully determined, which sets upper bounds on the allowed densities. The isotropic band structure model is obtained as follows: for all direct zincblende materials, the Γ-valley is assumed isotropic with  45,[69][70][71][72]74 The symbols have the following meaning: Eg(0 K) is the fundamental bandgap at zero temperature; me/m 0 is the density-of-states electron mass in units of the free electron mass; m h /m 0 is the density-of-states hole mass in units of the free electron mass; g e,h are the band multiplicities including the spin degree of freedom; μ * /m 0 is the reduced effective mass in units of the free electron mass; α e,h are the mass ratios μ * /m e,h ; γ e is the nonparabolicity parameter for the conduction band; Ry ex is the excitonic Rydberg energy, i.e., the binding energy of the exciton (used for energy scaling); aex is the excitonic Bohr radius (used for length scaling); and ϵs is the static permittivity of the semiconductor in units of the vacuum dielectric constant. The excitonic quantities are very suitable reference units, even at densities where excitons do not exist. 1 Detailed explanations and corresponding formulas are given in Sec. II of the main text.

No.
Material Thanks to the isotropic model, the angular integrations in the correlation energies (7) and (8) and in the dielectric function (9) can be performed exactly. The numerical expense then reduces to a double integral plus the Matsubara sum for the free-carrier correlation energy and to a twofold integration for the ionic correlation energy. The Matsubara sum is transformed into an integral in the case of T = 0 K. For zero temperature, the multiple integrals are solved approximately in Sec. IV to obtain a general analytical BGN model in terms of the material parameters listed in Table I. An important aspect in the context of BGN is the possible maximum density. Applying the theory to unreasonably high densities can readily shrink the gap of certain n-type materials to zero. In quasi-neutral regions, the carrier density is equal to the electrically active doping. The intrinsic limitations of the latter have been subject to intensive research. 47,48 It has been recognized that this limitation is a property of the material and not a chemical or electronic feature of the dopant 47 (e.g., formation of inactive complexes or charge compensation due to the amphoteric nature of shallow dopants). Process-related localized defects stabilize the Fermi energy at the Fermi-level stabilization energy ∼4.9 eV below the vacuum level in many III-V materials. 47,49 The amphoteric nature of these defects leads to a Fermi level pinning in the gap or in the CB, independent of the type or doping level, which explains the saturation of the electrically active doping. Experimental data for the limits of activated doping are only available for some of the 20 materials treated in this paper: InSb, 50,51 InAs, 52,53 GaSb, [54][55][56] InP, 57,58 GaAs, 59 GaP, 60-64 and zb-GaN. 48,[65][66][67][68] These limits were adopted for the density limits since they are, in any case, smaller than the limits imposed by the band structure model. In all other cases and, in general, for the neutral eh-plasma, the maximum densities were calculated from the condition that the Fermi level must not move into the next higher CB valley (electrons) or into the split-off band (holes). The valley separation energies E L − E Γ or E X − E Γ , and the spin-orbit energies Eso were adopted from Refs. 44, 45, and 69. The nonparabolic free dis-  Table I). The point m h = 1.369 m 0 for zb-AlN was omitted for better scalability.
persion (11) was used for the relation between electron density and Fermi level in n-type materials. The upper bounds of the density can be inferred from the end points of the density axis in Figs. 10-29. If available, the nonparabolicity parameters γ e were either directly taken from full-band calculations 72 or from the data collection of Ref. 44. In most cases, they had to be estimated by 73 γ e = 1/Eg (0 K), where the size of the bandgap at zero temperature was adopted from Refs. 44 and 45, except for InN, where a more recent value was used. 69,74,75 The static permittivities were taken over from Ref. 44. If they were not available for certain ternary compounds, they were calculated as the composition average of the permittivities of the binary constituents.
All material-related parameters needed for the BGN calculation are listed in Table I, where the size of the energy gap at T = 0 K was chosen for the ordering of the materials. The degeneracy factors g a account for the multi-valley conduction band (AlAs, GaP), heavy and light hole bands, and spin summation. It is convenient to use normalized quantities. All energies are normalized by the excitonic Rydberg energy Ryex = ̵ h 2 /(2 μ * a 2 ex ) and all lengths by the excitonic Bohr radius aex = h 2 ϵs/(e 2 μ * ), where μ * denotes the reduced effective mass μ * = (1/me + 1/m h ) −1 . The parameters αa = μ * /ma serve to transform excitonic quantities back to electron/hole quantities. Note that αe + α h = 1. Figure 2 shows permittivity (left ordinate) and effective DOS masses (right ordinate) as functions of material, where their ordering is the same as in Table I, i.e., according to the size of the energy gap at T = 0 K. The general trends are that the permittivity decreases from ∼17 to ∼9, and both DOS masses increase. In the case of the CB DOS mass, the indirect materials AlAs and GaP are outliers, which is related to their multi-valley nature. The nitrides are strong outliers with respect to the VB DOS mass (see the distinct peaks of the blue curve). The general trends of permittivity and DOS masses roughly explain the general trend of BGN (see Sec. III).

III. SUMMARY OF FULL-RPA RESULTS
The total full-RPA BGN as a function of material, i.e., in dependence on the increasing bandgap, is shown in Fig. 3 for 300 K and in Fig. 4 for 20 K. Two densities have been chosen for all 60 cases: 10 18 and 10 19 cm −3 . In some cases, the higher density already exceeds the and the cases n-type (red), p-type (blue), and eh-plasma (green). Lacking data points at 10 19 cm −3 are due to the density limitations described in the text. allowed limit, as discussed earlier; hence, the corresponding values were omitted.
Looking at the room-temperature results first, n-type, p-type, and eh-plasma differ only a little at the lower density (filled circles). The average trend is an increase in BGN from ∼15 to ∼50 meV. At the higher density, p-type and eh-plasma still behave similarly with an average increase from ∼30 to ∼100 meV, but the BGN in direct ntype material is significantly larger due to the small electron mass leading to reduced screening. 7 One observes an average increase from ∼40 to ∼120 meV. The indirect materials AlAs and GaP make an exception since their CB-DOS and VB-DOS are comparable in strength. Therefore, the BGN in n-type material is almost identical to the one in p-type material here. The trends closely follow the trends of the DOS effective masses (see Fig. 2). The maxima produced by the nitrides are reflected in the BGN curves, albeit scaled-down a lot.
Turning to the low-temperature case of 20 K in Fig. 4, the behavior of p-type material and eh-plasma is not only almost identical at the lower density of 10 18 cm −3 but also very similar to the one at room temperature. To a lesser extent, this also holds at a  higher density of 10 19 cm −3 . The behavior of n-type material, however, differs drastically from that at 300 K. At lower densities, the BGN exhibits a slightly decreasing tendency up to the indirect materials, followed by an increase for the wide-gap nitrides. Nevertheless, the BGN is relatively independent of material, taking values between ∼50 and ∼80 meV. In direct n-type materials, the BGN becomes dramatically increased as a consequence of the weak screening of the hole-ion interaction by majority carriers (electrons), an effect that gains importance with decreasing temperature (compare Fig. 5 and the detailed BGN representations in Appendix A). The analytical explanation for this temperature dependence will be given in Sec. IV. One can see an increasing tendency of the BGN from ∼130 to ∼200 meV, only interrupted by the indirect materials AlAs and GaP, where again the values are the same as in p-type materials due to the similar DOS size. In Figs. 5 and 6, the BGN of the small-gap direct material InSb is compared with that of the wide-gap indirect material AlAs. The density limit for n-type InSb in the left part of Fig. 5 stems from the measured maximum of activated doping (7.5 × 10 18 cm −3 ) 51 and results in a BGN that already slightly exceeds the value of the zerotemperature gap. Interestingly, for a limit of 4.5 × 10 18 cm −3 , as found in a standard low-pressure growth process, 50 the gap would exactly shrink to zero at T = 0 K. In all other cases in Figs. 5 and 6, the upper bounds of the density are imposed by the band structure model. The striking difference between InSb and AlAs is the behavior of n-type materials. Whereas in AlAs, the BGN is similar to that in p-type and eh-plasma (note the similar BGN scale and the almost identical temperature dependence, which ceases above 10 19 cm −3 ), the BGN in n-type InSb exhibits a strong temperature dependence that becomes even more pronounced at higher densities. At T = 0 K and 10 18 cm −3 , the BGN is 15 times larger than in p-type material. As outlined in Ref. 7, this is due to the small electron effective mass and the resulting weak screening of the hole-ion correlation energy. It will be shown in Sec. IV that the strong temperature dependence arises from the occupation probability of the minority (hole) band.
Appendix A comprises the presentation of BGN for all 20 materials (n-type and p-type) and for the neutral eh-plasma (in the intrinsic form of these materials) as a continuous function of density from 10 15.5 cm −3 to the maximum at the temperatures T = 0, 20, 77, 300, and 500 K. The full-RPA curves for T = 0 K   13), or (iv) the size of the spin-orbit energy Eso that determines n i,max . The latter is rather small [(6-19) meV in the ordering of the nitrides given in Table I], 45 which leaves only a narrow energy interval compliant with the 2-valley VB model. The ratio n h,M /n i,max increases from 5 to 12. Actually, (i) and (ii) are not independent of each other since n h,M ∼ (p M m h ) 3 . As the Mott parameter is related to the lattice constant, it decreases from zb-InN toward zb-AlN, which could partly explain the discrepancy. The analytical T = 0 curves are displayed over the whole density range, neglecting carrier freeze-out. Figure 7 compares the available experimental data with the full-RPA results. The left panel shows PL 31,32 and hetero-structure BT 33 data for p-type GaAs. They fairly match the theoretical curves. In the case of p-type GaSb (right panel), the measured BGN by PL 32 is systematically larger than the RPA one. The opposite is found for n-type InP (absorption data 35 ) in the high-density range. Absorption Interestingly, full-RPA results (based on the assumption of a random impurity distribution) are in reasonable agreement both at low and very high densities.

IV. ANALYTICAL MODEL WITHOUT MATERIAL-DEPENDENT FREE FIT PARAMETERS
In this section, a fully analytical model of BGN as a function of density and temperature in terms of the material and band structure parameters of Table I is presented. The goal is not the highest accuracy but to avoid free fit parameters for each case [as performed in Refs. 77 and 78 for T = 0 K and local-density approximation (LDA)]. This facilitates physical insight and the possibility to quickly calculate BGN for other materials or ternary compositions, as long as they obey (or can be approximated to) the isotropic band structure model consisting of one CB valley and two VB valleys. Wurtzite materials are excluded because of their complicated VB structure caused by the crystal-field splitting.
For p-type material and the neutral eh-plasma, this automatically interpolates between the high-and low-density regimes. In the case of direct n-type materials, the hole-ion correlation energy causes a strong temperature dependence at any density of interest. This is because the main origin is the occupation probability of the minority holes, which are non-degenerate no matter how strong the n-doping is. Therefore, ΔE g,T=0 is replaced by ΔE g,0 (T) with the analytical model for Δ i h,n−type (T) (see below), but the T = 0-limits in the case of all other energy contributions. The total BGN is, therefore, calculated by ΔEg(n, T) = 1 + U(n, T) 1/ΔE g,T→∞ (n, T) + U(n, T)/ΔE g,0 (n, T) .
Here, n equals the concentration n i of activated dopants in quasineutral regions or the plasma density ne = n h of the neutral ehplasma, respectively. Note that BGN is defined as a positive quantity, i.e., the negative sum of all QPS contributions (6)- (8). The Debye limits read 1,5 where the first term is the exchange energy and the second term is the sum of all correlation energies. Nonparabolicity is negligible in the low-density range.
The analytical calculation of ΔE g,0 (n, T) involves the following terms: The expression of the exchange energy Δ x {e,h} is well known, 1 and numerous attempts have been made to fit power laws in terms of density to the RPA free-carrier correlation energies of the ehplasma Δ c {e,h} (reviewed in Ref. 1). In Ref. 7, the ionic correlation energies Δ i {e,h} at T = 0 K were derived for n-and p-doped quasineutral regions of InGaAs, with a first estimate of how strongly nonparabolicity impacts Δ i h,n−type . In this paper, all terms are calculated analytically starting from Eqs. (7) and (8) using the RPA dielectric function (9). A central role in the multiple integrations of complicated functions is played by the length where q F,b is the Fermi momentum of the degenerate majority carriers, and "b" is the index of the majority carrier band. Hence, s b is proportional to the Fermi wavelength λ F,b in units of the excitonic Bohr radius, The assumption s b ≪ 1 is sometimes used in this paper to simplify the integrands. Since 1/q F,b is normalized by the excitonic Bohr radius aex, this assumption works better as the Fermi level moves deeper into the majority band (increasing density, weak DOS) and the excitonic Bohr radius is large, i.e., in the case of direct n-type III-Vs with small electron effective mass me. On the contrary, if g b is large and also me, as in the case of n-type AlAs and GaP, the assumption s b ≪ 1 is rather poor, as s b = 1 already corresponds to a density of ∼ 10 21 cm −3 . This finally shows up in a worse agreement between analytical and numerical ΔEg(n, T)-curves. Therefore, also s b ≫ 1 is tested and utilized if more appropriate. The analytical form of BGN is first given in terms of the length s b , which is more compact, and then as a function of the normalized density for more convenient use.

A. Quasi-neutral regions
In quasi-neutral regions, ΔE g,0 (n, T) reads in terms of the normalized length s b In Eq. (20), "b" is the index of the majority carrier band (b = e for n-type, b = h for p-type material), and "a" is the index of the minority carrier band (a = h for n-type, a = e for p-type material). The ordering of the five terms is the same as in Eq. (17). Note that the third term, the minority free-carrier correlation energy, only contains the index of the majority carrier band (b), since minority carriers do not contribute to screening. Nevertheless, the effective mass of the minority carriers finally comes into play after multiplication with the normalization energy (excitonic Rydberg). Equation (21) is the nonparabolicity correction that only applies to n-type material. Note that it increases the electron-ion correlation energy (fourth term) but decreases the hole-ion correlation energy (fifth term) with different powers. Equation (22) describes the reduction of the hole-ion correlation energy with increasing temperature at intermediate and high densities. It also applies to n-type material only and is strongly affected by nonparabolicity (note the occurrence of φ np in the argument of arctan[⋅ ⋅ ⋅]). The zero-temperature limit (β → ∞) is φ T = 1. Figure 8 shows the ability of this model to reproduce the T-dependence of the hole-ion correlation energy for all densities >10 15.5 cm −3 over the whole temperature range. By way of example, the cases InSb (small gap), InGaAs (intermediate gap), and zb-GaN (wide gap) have been chosen. One observes a decreasing T-dependence and a changing shape of the n-dependence with a rising gap.
From Eq. (8), one can see that the temperature dependence has two origins: (i) the distribution function of minority carriers in the factor (∂n h /∂ζ h ) −1 and (ii) the occupation probability of the majority carriers contained in the dielectric function (9). From a second-order Taylor expansion of the latter by 1/β around the Fermi energy, one obtains a quadratic T-dependence which, however, is weak and vanishes at high densities. It is negligible compared to (i), which results in the term (22). The derivation of this term is given in Appendix C.
Replacing s b by the normalized density n using Eqs. (18) and (19), i.e., with In Eqs. (23)-(25), n represents the density n i of ionized dopants (normalized by a −3 ex ). Note that BGN is normalized by the excitonic Rydberg energy Ry ex , as is the thermal energy 1/β and the inverse of the nonparabolicity parameter 1/γ e .

B. Neutral electron-hole plasma
In the case of the neutral eh-plasma, both exchange energies are present, but the carrier-ion correlation energies are absent [see Eq. (17)]. The analytical model for T = 0 K and parabolic bands reads as a function of the normalized lengths s {e,h} , with The first term in (26) is the sum of the exchange energies; the second and third terms are the correlation energies of electrons and holes, respectively. These terms differ from the corresponding majority-   In Fig. 9, the numerical RPA results of the band edge shifts induced by the free-carrier correlation in InGaAs at T = 0 K are compared with the outcomes of the analytical models. The freezeout of carriers was disregarded to display the whole density range. The agreement is excellent for the shift of the CB edge and reasonable for the shift of the VB edge. This difference in accuracy is due to se ≪ s h (αe ≫ α h ). Figure 9 also shows that the correlation energies in a neutral eh-plasma are very similar to the free-carrier correlation energies of the majority carriers in quasi-neutral regions. Size and slope are only slightly larger, which can be attributed to the stronger screening in a two-component plasma.
The comparison of the analytical BGN model with the numerical full-RPA is presented in Appendix A for all 60 cases (Figs. 10-29). In p-type material and in the eh-plasma, the temperature dependence is only significant at intermediate densities. For better visibility, the analytical curves were plotted just for 0 K and room temperature here. The Padé approximation leads to a relatively large deviation in the transition region, 1,41 which is unavoidable. The density and temperature dependences are correctly reproduced in all
60 cases. Very good agreement is found for the direct n-type materials, except for InAsSb at high densities because of its extreme nonparabolicity. Larger deviations occur for the indirect materials AlAs and GaP, as se ≪ 1 is a poor assumption here. For many p-type materials, the accuracy is satisfactory as well; exceptions are  InAsSb, InSb, and the nitrides. For the neutral eh-plasma, the analytical model works best at high densities but deviates markedly at intermediate densities due to the error caused by the Padé interpolation.

V. CONCLUSION
The BGN in quasi-neutral n-type and p-type regions and in a neutral electron-hole plasma of 20 technologically relevant zincblende III-V semiconductors was calculated in full-RPA as a function of density in the temperature range from 0 to 500 K. The sevenfold (plasma) and sixfold (quasi-neutral regions) integration is enabled by a careful isotropization of the band dispersion, including the important CB nonparabolicity in n-type materials. The neglect of higher CB valleys and of the split-off VB imposed by the simplified band structure model results in upper bounds for the allowed density in each case. The available experimental limits of the activated doping are always smaller than these theoretical limits. The gap can never shrink to zero, even at zero temperature. An exception is n-type InSb, where measured doping limits would lead to a vanishing gap at 0 K.
At room temperature, the trend of the BGN as a function of material (gap) follows the trend of the DOS effective masses. At a density of 10 18 cm −3 , the BGN increases from ∼15 meV (InAsSb) to ∼50 meV (AlN) with only a little difference between n-type, p-type, and eh-plasma. At the higher density of 10 19 cm −3 , the BGN REVIEW pubs.aip.org/aip/adv of p-type and eh-plasma is still similar with an average increase from ∼30 meV (InAsSb) to ∼100 meV (AlN), but the BGN in direct n-type material becomes larger due to the small electron mass (reduced screening) and increases from ∼40 meV (InAsSb) to ∼120 meV (AlN). An exception is given by the indirect materials AlAs and GaP due to their similar CB-and VB-DOS. The nitrides are strong outliers with respect to the VB DOS mass, which becomes visible as maxima in the BGN curve. At very low temperatures, p-type material and eh-plasma show almost identical behavior at a density of 10 18 cm −3 , which is not much different from that at 300 K. At the higher density of 10 19 cm −3 , p-type material and eh-plasma are still comparable, but in direct n-type materials, the BGN becomes giant as a consequence of the weak screening of the hole-ion interaction by the majority of electrons. At a density of 10 18 cm −3 , the BGN is relatively independent of material, taking values between ∼50 meV and ∼80 meV at 20 K. At 10 19 cm −3 , it increases from ∼130 meV (InAsSb) to ∼200 meV (AlN), again with the exception of the indirect materials AlAs and GaP.
The strong temperature dependence of the BGN in n-type materials with small CB DOS is traced back to the occupation probability of minority holes. It persists at all electron concentrations. In p-type material and eh-plasma, the temperature dependence of the BGN is much weaker; it is only significant at low and intermediate densities and always ceases at high densities.
The second major result of the paper is the derivation of an analytical BGN model without material-dependent free fit parameters, applicable to all densities and temperatures. Compared to numerical full-RPA, it reduces the central processing unit (CPU) time from days (the worst case is a very low but non-zero temperature and a high density) to less than seconds. The analytical model reveals the dependencies on band-structure and material parameters and allows for a quick calculation of the BGN in materials not explicitly treated in this paper. The only restriction here is compliance with an isotropic dispersion model (which excludes wurtzite materials). Besides other model simplifications (e.g., neglect of disorder-induced DOS tails and polaronic mass enhancement), isotropization is the key not only for the numerical full-RPA calculations but also for the analytical model. After the exact angular integrations, threefold integrals over complicated functions are approximately solved for the free-carrier correlation energies at 0 K and for the non-perturbative temperature-dependent hole-ion correlation energy, including the nonparabolicity of the CB. This is enabled by taking advantage of the ratio between Fermi wavelength and Bohr radius and by applying other approximation methods. As a result, the functional dependence of the BGN(n,p,T) on material and band structure parameters can be preserved. Comparison between analytical and numerical curves shows very good agreement in many cases, in particular for the density and temperature dependence of the BGN in n-type materials. Larger deviations occur in the case of n-type InAsSb at high densities (extreme nonparabolicity), for the indirect materials AlAs and GaP (due to their comparable CB and VB DOSs), and for some p-type materials such as InAsSb, InSb, and the nitrides.
A serious attempt has been made to express the BGN of ternary compounds by a composition average of the BGNs of their binary constituents (without or with bowing), as is common for the bandgap. However, although possible to some extent for a particular density and temperature, it fails for a wider density interval. This is probably due to the involved dependence on band filling and screening.

Conflict of Interest
The author has no conflicts to disclose.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Quasi-neutral regions
The zero-temperature correlation energy of free carriers is calculated by Eq. (7) in the parabolic band approximation (justification given in the main text). The integration of angles can be performed exactly. Furthermore, the limit β → ∞ is utilized. First, in the factor ∂ϵ/∂n, the derivative of the distribution function of majority carriers can be generated by partial integration, which turns into a deltafunction and removes one integral. Second, the Matsubara sum is replaced by ∑ ν → 2β/π ∫ ∞ 0 dt with t = πν/β. In the momentum integration, the variable q is normalized by the Fermi momentum of the majority carriers: q → zq F,b = z √ ζ b /α b . This results in REVIEW pubs.aip.org/aip/adv The RPA dielectric function (9) becomes fully analytical at T = 0 K, with the normalized length s b defined in Eq. (18), Note that only majority carriers (index "b") contribute to free-carrier screening in quasi-neutral regions. Expression (B2) is obtained from Eq. (A6) in Ref. 7 in the parabolic and zero-temperature limits. It must be drastically simplified to enable the double integration of momentum (z) and frequency (t) in (B1). For that purpose, the single plasmon pole (SPP) approximation 3,10,39,40 is a good candidate.
Here, it is used for a one-component plasma, where ωp is the plasmon frequency and κ is the Thomas-Fermi momentum (the static Debye limit is ϵ D = 1 + κ 2 /q 2 ). After normalization, one obtains for the case of a degenerate carrier gas at T = 0 K, and ϵspp(z Note that ϵspp is not used in the factor ∂ϵ/∂n of Eq. (B1) but only in the screening function (1 − 1/ϵ), which reads The correlation energy (B1) takes the form where the integration variable t was changed to τ = t/ζ b . The goal is to approximate the integrand such that the correct dependence on density (contained in s b ) is preserved as well as possible. After testing various ways, the following turned out to be the best.

a. Majority-carrier band edge
For the majority-carrier band edge (first term in curly braces), the frequency dependence in the denominator of the SPP screening function [first line in Eq. (B9)] is skipped (τ → 0). This allows us to perform the remaining τ-integration exactly, leading to The numerator of the integrand is replaced by zΘ and after exact integration in the analytical expression To obtain the density dependence in the form of a simple power law, the limits s b ≪ 1/6 and s b ≫ 1/6 are considered, The first limit (high density, short Fermi wavelength) yields an n 1/6 -dependence, and the second limit (low density, long Fermi wavelength) yields an n 1/4 -dependence. In Fig. 30, se and s h at a density of 10 19 cm −3 are plotted for all materials. Neither of the two limits is actually valid, as one can see by comparison with the green dashed line (1/6). A possible simple fit to the full-RPA curves (compare Fig. 9) is the following: (i) the limit s b ≪ 1/6 is chosen, (ii) a logarithmic weakening 5,36,38 is performed, and (iii) the numerical pre-factor is adjusted. The result is [the second term of Eq. (20)], i.e., a logarithmic weakened n 1/6 -dependence.

b. Minority-carrier band edge
For the minority-carrier band edge [second term in the curly braces of Eq. (B9)], the τ-integration can be performed exactly, and one is left with the single integral, The integrand is replaced by a simpler and integrable function that yields the correct limits z → 0 and z → ∞ and that closely follows the course at intermediate z, Except for the occurrence of α b in the last term of the denominator, this approximate integrand equals the integrand in Eq. (B11). Proceeding as for the majority-carrier band edge, one obtains A feasible power-law fit to the full-RPA curves (compare Fig. 9) is obtained as follows: (i) the limit s b ≫ 1/6 is chosen, and (ii) the numerical pre-factor is adjusted. The result is [the third term of Eq. (20)], i.e., an n 1/4 -dependence.
It should be noted that the power-law fits for majority and minority band edges are somewhat arbitrary, as the applicability of the limits strongly depends on material and density. This might explain the many variants suggested in the past, e.g., n 1/6 in Ref. 37, n 1/4 in Ref. 38, and n 7/30 in Refs. 38 and 41.

Neutral electron-hole plasma
Screening by electrons and holes with density n = ne = n h determines the free-carrier correlation energies in a neutral ehplasma. Analytical modeling proceeds similarly to the case of a one-component plasma. Charge neutrality imposes a simple relation between the quasi-Fermi levels (parabolic dispersion, T = 0 K), In the screening function (1 − 1/ϵ), the SPP approximation is now applied to a two-component plasma, (B20) Using relation (B19), plasmon frequency ωp and Thomas-Fermi momentum κ at T = 0 K are given by The screening function reads with α ba = αa The form of the correlation energy at zero temperature is similar to the majority-carrier term in Eq. (B1), (B25) With the same approximations as for Eqs. (B10) and (B11), one obtains where α ba is defined in Eq. (B24). Exact integration yields the analytical expression REVIEW pubs.aip.org/aip/adv Again, the limits of very small and very large s b are considered in order to simplify the model, As can be seen from Fig. 30 The outcome of this general analytical model is displayed for the free-carrier correlation energy in InGaAs (Fig. 9) and as part of the total BGN of the neutral eh-plasma in all materials in the right panel of Figs. 10-29.

APPENDIX C: ANALYTICAL APPROXIMATION OF THE ION-CARRIER CORRELATION ENERGIES
This appendix provides the derivation of the analytical expressions of the ion-carrier correlation energies in quasi-neutral regions including the effect of nonparabolicity in n-type material, and the temperature dependence of the donor-hole correlation energy including the effect of nonparabolicity on this temperature dependence.
The correlation energy of the ion-carrier interaction is calculated from Eq. (8). After normalization and integration over angles, it is given by The static RPA dielectric function reads [compare Ref. 7, Eqs.
where ha(κ, q) = have to be used. For the derivative of ϵ(q, 0) with respect to the Fermi energy, one needs under the κ-integral. To obtain the T = 0 K-limit, the κ-integral in Eq. (C2) is integrated by parts, where the derivative of the distribution function with respect to κ becomes a delta function, Then For the derivative of the static dielectric function with respect to the Fermi energy ∂ϵ(q, 0)/∂ζ a in the limit β → ∞, one needs which follows from Eq. (C5) neglecting nonparabolicity for the minority carriers. Using this in Eq. (C2) leads to where in the term ∼ Θ(−ζ a ) nonparabolicity was again neglected. The κ-integral can be calculated analytically, For the second line, the relation ln(|x + 1|/|x − 1|) = 2 tanh −1 (1/x) was used. Inserting Eq. (C10) with (C11) and Eq. (C8) into Eq. (C1), the correlation energy of ion-carrier interaction in the T = 0 K-limit becomes where the static dielectric function ϵ T=0 (q, 0) is defined in Eq. (C7) and the function ha in Eq. (C3). Introducing normalized momentum variables of majority carriers by z = q/q F,b = z √ α b /ζ b and t = κ/ √ζ b , this can be written in the form after inserting (C7) into Eq. (C12). Analytical integration necessitates the simplification of the function H b (t, z). The shape of the total integrand strongly depends on the size of the parameter s b . In terms of the normalized doping concentration (normalized by a −3 ex ), the parameter s b becomes [compare Eqs. (18) and (19)], Therefore, one can assume s b ≪ 1 for very large n i , which results in a pronounced maximum of the integrand at some z < 1. The function H b (t, z) can then be linearized in z, which results in With that, the t-integration in the dielectric function can be performed easily, giving The second term is negligible compared to the first one as long as 4γ bζb ≪ 1 can be assumed, meaning that the Fermi energy in the corresponding parabolic band must be much smaller than the bandgap. For InAsSb and InSb, this becomes critical at higher densities, as can be seen in Figs A last simplification is the replacement that closely matches the numerical RPA result. Using the notation φ np (s b ) = h b (1, 0) for the nonparabolicity correction factor finally results in the ionic correlation energies in Eq. (20). Note that γ h = 0. As discussed in Appendix B 1, the strong temperature dependence of the hole-ion correlation energy Δ i h,n−type (T) mainly arises from the occupation probability of holes, whereas the majoritycarrier screening is only weakly T-dependent and vanishes at high concentrations. Therefore, Eq. (C1) becomes Since γ h = 0 is used throughout, nonparabolicity only needs to be taken into account in the factor 1/ϵ 2 T=0 (q, 0) (screening by majority electrons). Furthermore, the T-dependent factors (∂n h /∂ζ h ) −1 and ∂ϵ(q, 0)/∂ζ h can be calculated with Boltzmann statistics. From Eq. (10), one obtains after normalization and from Eq. (C2) after partial integration with H h (κ, q) defined in Eq. (C14). Its parabolic version reads ). Inserting into Eq. (C23) and introducing the new integration variables t = κ √ β and z = q √ αe/ √ζ e, consequently, where τ(z, β) = ∫ za(β) a 2 (β) = α hζe 4αe β = α h g 2 e (2παese) 2 β. (C28) The hole-ion correlation energy can now be written as . (C29) Its zero-temperature limit 7 follows from lim β→∞ τ(z, β) = 2 √ παe/(α hζe z 2 ). In the low-T/high-density range, the denominator in (C27) can be taken at t = 0, and He(t, z) can be linearized in z: He(t, z) ≈ 2zhe(t, 0)/t. The dielectric function then reads ϵ T=0 (z, 0) = 1 + s e z 3 ∫ 1 0 dt t He(t, z) = 1 + 2sehe(1, 0)/z 2 , which results in For an analytical solution of the integral in (C30), the error function is approximated as This results in where he(1, 0) was renamed by φ np (se), which is the symbol for the nonparabolicity correction factor Eq. (21) in the main text. The last term in curly braces in Eq. (C32) is small compared to the second one and will be omitted.
In this approximation, the temperature dependence has the form ∼ [1 − 2 arctan(c √ T/n 1/6 )/π]. Due to the assumption of very large Fermi momenta in the dielectric function, z ≪ 1, it only matches the numerical results at low temperatures and extremely high densities, which are outside the validity range.
The opposite assumption, z ≫ 1, would lead to He(t, z) ≈ 8t/z. The dielectric function then reads ϵ T=0 (z, 0) = 1 + 8se/(3z 4 ), resulting in a temperature and density dependence of the correlation energy given by ∼ √ n/T as in the Debye limit [compare Eq. (16)]. Combining both limits with the static form of the SPP approximation (B7) leads to cumbersome expressions that are not helpful. Therefore, to cover the cases z ≪ 1 and z ≈ 1, a viable solution is to REVIEW pubs.aip.org/aip/adv replace arctan[(⋅ ⋅ ⋅) 1/2 ] by arctan[(⋅ ⋅ ⋅) 3/4 ] in (C33). This not only preserves the zero-temperature limit but also results in excellent agreement with the numerical curves for all n-type materials over the whole density and temperature range, as can be seen from Fig. 8 and the left panel of Figs. 10-29.