Finite-temperature full random-phase approximation of bandgap narrowing in quasi-neutral regions: Theory including nonparabolicity and application to silicon and InGaAs

The bandgap narrowing (BGN) in quasi-neutral regions of semiconductors is calculated in a finite-temperature full random-phase approximation formalism based on a simple isotropic dispersion model including band nonparabolicity. The total quasi-particle shift (QPS) is determined by the exchange-correlation self-energy of the free carriers and the correlation energy of the interaction between carriers and ionized dopants. At cryogenic temperatures, the latter part results in giant shifts of the minority band edge in n-type semiconductors with a large ratio of valence to conduction band density of states, as often present in III – V materials. However, at room temperature, the BGN does not exceed common values. The reason for this behavior is explained analytically. Whereas the exchange-correlation energy of free carriers is known to be insensitive to band structure details, the nonparabolicity of the conduction band (CB) has a strong effect on the ionic QPS of the minority carriers in n-type III – V materials. It strongly reduces the BGN at cryogenic temperatures compared to the case of a parabolic CB. This is demonstrated numerically and also analytically for n-type InGaAs. The BGN in n-type silicon becomes independent of temperature at high concentrations, but in p-type silicon, a weak temperature dependence re-emerges above the Mott density, which also can be attributed to the ionic QPS of the minority electrons. The calculated BGN for quasi-neutral regions in silicon is in good agreement with earlier photoluminescence and more recent photo-conductance measurements.


I. INTRODUCTION
Bandgap narrowing (BGN) is an important phenomenon in semiconductor devices. It is a consequence of either heavy doping or strong optical excitation and high electrical injection. High excitation/injection levels lead to plasma-induced BGN, as in photo-conductive switches, concentrator solar cells, and power devices operated in the on-state. The plasma can be neutral, e.g., in the case of laser excitation of intrinsic regions, or not, e.g., in the case of optical excitation of doped regions. Quasi-neutral regions of electronic devices (emitters of bipolar transistors, source/drain-regions of field-effect transistors) host a one-component plasma that compensates the electrically active dopants. "Quasi" refers to the possible presence of a very small electric field when the device is biased or a small concentration of minority carriers. The manybody Coulomb interaction between free carriers as well as between carriers and ionized dopants results in self-energy and a corresponding shift of the band edges such that the bandgap shrinks. In quasi-neutral regions, one has five components: the exchange energy of the majority carriers, the free-carrier correlation energies of majority and minority carriers, and the ionic correlation energies of majority and minority carriers due to their interaction with the electrically active dopants. These contributions can be treated equally in the theoretical framework of the finite-temperature full Random-Phase Approximation (RPA). [1][2][3][4][5] All measurements of BGN involve electron-hole pairs. In quasi-neutral regions, the minority carrier band shifts because the generated/annihilated minority carrier lowers its energy due to Coulomb attraction by the mobile majority carriers and Coulomb repulsion by the sub-system of immobile ionized dopants that are the source of the majority carriers. 6 As will be shown in this paper, the ionic correlation energy of holes in n-type In 0:53 Ga 0:47 As is even the largest contributor to the total BGN in this material.
Most experimental data on BGN exist for silicon but with a rather large spread, which can be attributed to the involved evaluation procedure. Absorption experiments [7][8][9][10][11][12] require the admixture of phonon-assisted transitions and the knowledge of phonon energies and Fermi level. Problems in photoluminescence (PL) experiments [13][14][15][16] are the nature of the initial states and the weak intensity of the spectrum. The electrical method [17][18][19][20][21][22] uses bipolar transistors where the collector current as a function of the emitterbase voltage, the sheet resistance underneath the emitter, and the minority carrier mobility in the base must be measured. The BGN is then found from an analytical model that connects these quantities and that itself was derived with a number of approximations. For materials other than silicon, BGN data are rare, and for most materials, they are not available at all. Therefore, reliable calculations and fit models derived from them play an important role to fill this gap.
Earlier theoretical methods to derive the self-energy of the carrier-dopant interaction in the T ¼ 0-limit comprise Hartree-Fock variational calculations for a donor system distributed on a regular sub-lattice 23 using the single plasmon pole approximation 24 and application of second-order perturbation theory for a random system of dopants 25 using RPA screening. The role of multi-valley scattering in both dopant arrangements had been the topic of discussions. [26][27][28] RPA is the method of choice to compute BGN for the whole temperature range, and it gains accuracy with an increase in the concentration above the Mott density. The theoretical basics can be found in Ref. 1, which also presents the temperaturedependent QPS in a neutral electron-hole plasma. The case of an extrinsic semiconductor with plasma excitation was treated in Ref. 3 where numerical results based on a random distribution of dopants were given for T ¼ 0 K. In Ref. 29, an RPA-based fit model for silicon device simulation as a function of carrier densities, doping concentration, and temperature was developed. To cover all cases (neutral plasma, quasi-neutral region, depletion region), compromises had to be made in the fitting procedure. As will be shown below, this fit model underestimates the BGN in quasi-neutral regions at very high doping concentrations. Besides the revaluation of BGN for quasi-neutral regions in Si including nonparabolicity of the conduction band (CB), BGN in quasi-neutral In 0:53 Ga 0:47 As will be analyzed in detail. The main finding is that at a temperature of 300 K and as a consequence of strong CB nonparabolicity, the BGN in the n-type material is limited to %10% of the bandgap at the highest measured concentration of electrically active doping. In the p-type material, this amount is %6%.
The paper is organized as follows: Section II gives a brief summary of the RPA theory and shows how nonparabolicity is included. Section III provides results for silicon: numerical curves for the different ionic QPSs, the temperature-dependent total BGN in comparison with absorption and PL data, and the small influence of CB nonparabolicity. To better understand the ionic QPS of the minority and majority bands, they are calculated analytically as a function of doping concentration in the parabolic, T ¼ 0 K-limit. Section IV presents the results for In 0:53 Ga 0:47 As: numerical freecarrier and ionic correlation energies over the whole temperature range, the comparison of parabolic and nonparabolic curves, and the total BGN in n-type and p-type materials as a function of temperature and doping. To better understand the impact of CB nonparabolicity, the ionic correlation energy of minority holes is derived analytically in the T ¼ 0 K-limit. Conclusions are given in Sec. V. Appendix A contains the derivation of the full-RPA dielectric function for nonparabolic bands, its static limit, as well as the high-T/low-density and low-T/high-density limits of the static form, respectively. Appendix B demonstrates the derivation of the full-RPA correlation energy of ion-carrier interaction for nonparabolic bands in the high-T/low-density and low-T/high-density limits, respectively.

II. THEORETICAL BACKGROUND
BGN in RPA quality is calculated by the quasi-particle shift (QPS) 1 Δ a (k) (a ¼ e for electrons, a ¼ h for holes), which is the difference between free and interacting dispersion (the vector character of momenta and coordinates is only written explicitly where needed, and the normalization volume is set to unity), It represents the real part of the self-energy Δ a (k) ¼ Re Σ a k, E a (k) þ i0 þ ½ and consists of three terms where Δ x a (k) is the unscreened exchange term, Δ c a (k) is the correlation term of free carriers, and Δ i a (k) is the correlation term of the interaction between carriers and ionized dopants. The latter are assumed to be randomly arranged on regular lattice sites. The dispersion of the QPS is rather flat in the relevant energy interval between band edge and Fermi energy for all densities and temperatures. 1,29 Therefore, the dispersive QPS Δ a (k) can be replaced by a rigid shift Δ a . There is no need to solve the self-consistency problem then because the only occurrence of Δ a is in the distribution functions, where it is fixed by the given density in the quasineutral region. An elegant way to compute the rigid shift follows from the requirement that the QPS density should not change in the first order with respect to Δ a (k) À Δ a , which leads to 1 with rigidly shifted bands E a (k) ¼ E 0 a (k) þ Δ a and Fermi functions f a (k) depending on shifted chemical potentials μ a , The weight function in (3) (the derivative of the Fermi function) 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 explicit expressions for the three components of the dispersive QPS in (2) can be found in Refs. 1, 3, and 29 and are not repeated here. After insertion into Eq. (3), the rigid shifts take the compact form 1,3 In the case of parabolic bands, the QP density n a in the rigid shift approximation is given by with the Fermi integral F 1=2 , the thermal wavelength Λ a ¼ (2π h 2 β=m a ) 1=2 , the inverse thermal energy β ¼ 1=k B T, and the Fermi energy ζ a ¼ μ a À Δ a . The multi-valley conduction band, heavy-and light-hole bands, and spin summation are condensed in degeneracy factors g a . The effective masses m a are understood as density of states (DOS) effective masses m e ¼ (m 2 t m l ) 1=3 for aniso- for holes. This yields a tractable isotropic dispersion model needed for analytical BGN modeling that also allows one to include nonparabolicity. More comments on the band structure model will be given in Secs. III and IV. The remaining ingredients in Eqs. (5)-(7) are the bare Coulomb potential v(q) ¼ e 2 =(ϵ 0 ϵ s q 2 ) with the static permittivity ϵ s , the Matsubara frequency Ω ν ¼ 2πi ν hβ (ν integer), and the concentration of electrically active dopants n i . The correlation energies are determined by the RPA dielectric function As the rigid QPS drops out in the energy differences, the dispersion can be replaced by the free dispersion E 0 a . It should be noted that the single-integral form of the exchange term (5) is bound to parabolic bands E 0 a (k) ¼ h 2 k 2 =2m a , 1,30 whereas (6) and (7) are obtained by (3), regardless of the specific form of the dispersion. It is well known that band structure details have only a minor impact on the exchange-correlation energy of the electron and hole plasmas due to a compensation effect. 1,31 Therefore, Eqs. (5) and (6) are evaluated with parabolic bands in this paper. However, as will be shown below, in the low-T/high-density limit, the term Δ i h can take large values in n-type III-V materials. In this case, the Fermi level is located high in the CB, where the dispersion strongly deviates from the parabolic shape. To study how much nonparabolicity changes the BGN, electron density n e and dielectric function ϵ(q, Ω ν ) in (7) are calculated with the nonparabolic free dispersion where γ e denotes the nonparabolicity parameter. Both dielectric function and ionic QPS can be expressed by single integrals. This is demonstrated in Appendixes A and B, respectively. The T ¼ 0 K-limits are also given explicitly there as they serve to validate the tricky numerical solutions at very low but finite temperatures and to explain the high values of Δ i h in n-type materials with a large ratio of valence to conduction band density of states. It is convenient to use normalized quantities. All energies are normalized by the excitonic Rydberg energy Ry ex ¼ h 2 =(2μ * a 2 ex ), and all lengths by the excitonic Bohr radius a ex ¼ h 2 ϵ s =(e 2 μ * ), where μ * denotes the reduced effective mass μ * ¼ (1=m e þ 1=m h ) À1 . Useful parameters are α a ¼ μ * =m a . All values for Si and In 0:53 Ga 0:47 As are given in Table I.

III. RESULTS FOR SILICON
Application of the empirical nonparabolicity model (10) to silicon necessitates an isotropic DOS. The use of DOS effective masses m e ¼ (m 2 t m l ) 1=3 to take account of the anisotropy of CB valleys is already an approximation in the correlation energies (6) and (7) with the RPA dielectric function (9). Starting with ellipsoids and transforming k and q to elliptical coordinates would not give the same expressions because of the occurrence of the Fourier transformed Coulomb potential v(q). The light and heavy hole masses m lh , m hh are angle-averaged parameters, the light-hole band is nonparabolic, and the split-off band with its maximum at 44 meV below the VB edge is neglected. The latter simplification sets a limit to the hole density, i.e., calculated BGN values become rather uncertain above 5 Â 10 19 cm À3 . Heavy doping also induces band tails, an effect completely ignored in this paper. Figure 1 shows the ionic part of the average QPS Δ i e for neutral n-and p-type silicon, and Fig. 2 shows the ionic part of the average QPS Δ i h for neutral n-and p-type silicon as a function of activated doping concentration and temperature. Note that the Journal of Applied Physics While in a neutral electron-hole plasma the contribution of Δ i a to the total BGN becomes small under conditions of high excitation/injection, this is not the case in heavily doped quasi-neutral regions. The reason is that only majority carriers, and only with a fixed density n i are available for screening there. If their DOS is small, the screening effect is weak, and Δ i a of the minority carrier band can become significant. Despite the relatively large DOS of Si, this is clearly seen for p-type silicon in Fig. 1(b) and, on a smaller scale, for n-type Si in Fig. 2(a).
To understand the different behavior of majority and minority carriers it is instructive to derive the ionic correlation energy in a fully analytical form. This is only possible for T ¼ 0 K.

Journal of Applied Physics
The nonparabolicity effect is minor in Si because of the relatively small value of the nonparabolicity parameter (γ e ¼ 0:5=eV) and the large DOS of the CB that prevents that the Fermi level can move deep into the band with growing density. Therefore, in this section, the analytical treatment is restricted to the parabolic case.
The nonparabolic case will be presented for InGaAs in Sec. IV. The starting point is the parabolic and zero-temperature limit Eq. (B10) in Appendix B with the dielectric function Eq. (A13) in Appendix A. Inserting Eq. (A13) into Eq. (B10) yields by substitut- Here, the Fermi momentum of the majority carriers q F,b is defined by (13) and "b" is the index of the majority carrier band. For majority carriers [first term in (11)], an analytical approximation of the integral is difficult to find. The function tanh À1 ( Á Á Á ) has an integrable singularity at z ¼ 2 (double the amount of the Fermi momentum). The function in angular braces is bell-shaped with a maximum of 2 at z ¼ 0. 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 Therefore, one can assume s b ( 1 for very large n i , which results in a pronounced maximum of the integrand in (11) at some z , 1. The function tanh À1 ( Á Á Á ) can then be linearized and one obtains The last expression is a reasonable analytical approximation of the zero-temperature limit of the ionic QPS in terms of band structure parameters (g a , α a ), reproducing both the slope of the doping dependence at high concentrations and the different amounts of minority and majority carrier band shifts. Of particular interest is the ratio of the minority band shift in the material with the lower DOS Δ i e,pÀtype to the other shifts Δ i h,pÀtype , Δ i e,nÀtype , and Δ i h,nÀtype . Since the DOS is given by a . Therefore, with the parameters of Table I, the DOS of p-type silicon is about three times smaller than that of n-type silicon. Using Eq. (15), one arrives at the following ratios: The numbers at n i ¼ 10 correspond to a doping concentration of 1:94 Â 10 20 cm À3 and the numbers in braces are the numerical values extracted from the black solid curves in Figs. 1 and 2 for this density.
The total BGN in (a) n-type and (b) p-type silicon including the weak nonparabolicity effect is shown in Fig. 3. The temperature dependence starts to vanish above the Mott density, but a slight temperature effect re-emerges in p-type silicon at high densities (n i . 1). It originates from the contribution of Δ i e,pÀtype [see Fig. 1(b)]. Figure 4 shows a comparison of the theory with optical BGN measurements in n-type silicon (a) and p-type silicon (b). Data points are from Aw et al. 11  The black dashed curves are the result of the fit model for device simulation proposed in Ref. 29. In this model, screening had been fitted for the case of a neutral electron-hole plasma (n e ¼ n h ) to also cover the cases of high optical excitation and strong electrical injection in pn-junctions. Application to quasineutral regions yielded an acceptable agreement with the rather scattered experimental data, and no re-fitting was done for this case. As can be seen, the fit model underestimates the BGN in quasi-neutral regions of silicon by 30 meV around a concentration of 10 20 cm À3 . The reason is again the ionic QPS of minority carriers, which in a neutral plasma is significantly suppressed by the screening effect of the species that are minority carriers in quasi-neutral regions.
IV. RESULTS FOR IN 0.53 GA 0.47 AS To provide a more complete picture for InGaAs, first the freecarrier correlation energy Δ c a is depicted in Fig. 5. The majority carrier contribution is strongly temperature-dependent, but its magnitude remains small, in particular, in n-type InGaAs. The free-carrier QPS of minority carriers is in the order of 10 meV at 10 19 cm À3 , the largest possible concentration of electrical active doping. This maximum concentration corresponds to n i ¼ 100 (a À3 ex ¼ 1:2 Â 10 17 cm À3 ).
Due to the relatively small CB DOS of InGaAs, the Fermi level can move deep into the band at high n-doping. The nonparabolic dispersion then has a strong effect on the ionic QPS. This is demonstrated in Figs. 6 and 7. For electrons in the n-type material [ Fig. 6(a)], Δ i e,nÀtype doubles for nonparabolic CB at the highest density. As expected, in the p-type material, the nonparabolicity effect is much weaker [ Fig. 6(b)]. In general, Δ i e is limited to 5 meV, and the temperature dependence is only relevant at intermediate concentrations.
The situation is completely different for holes. If they are minority carriers [ Fig. 7(a)], Δ i h,nÀtype becomes giant in the low-T/high-density range, if a parabolic CB is used. It reaches 500 meV at n i ¼ 100, hundred times more than Δ i h,pÀtype [ Fig. 7  With decreasing doping concentration and rising temperature the value of Δ i h,nÀtype quickly shrinks to "common" values. At room temperature, one finds 8 meV at n i ¼ 10 and 48 meV at n i ¼ 100.   , i.e., the value goes to zero and becomes independent of carrier type and doping species.
As can be seen in Fig. 7(a), nonparabolicity of the CB strongly decreases the value of Δ i h,nÀtype at low temperatures. An analytical explanation of the nonparabolicity effect is now given for T ¼ 0 K. In the important case of minority holes in the n-type material, one obtains instead of Eq. (11), second term, with α eζ e,np q , which follows from Eq. (B9) inserting (A12). Note that the Fermi energy ζ e,np is lower than ζ e,par of the parabolic case ζ e,np ¼ 1 2γ e ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4γ e ζ e,par q À 1 % ζ e,par 1 À γ e ζ e,par : The second term in the denominator of Eq. (18) is treated in the same way as before, i.e., approximating H e (t, z) for small z. This results in where γ eζe,np ( 1 had to be assumed. The ionic correlation energy of minority holes (18)  : The ratio between nonparabolic and parabolic shift is computed using relation (20),

ARTICLE scitation.org/journal/jap
The number in braces is the numerical value extracted from the solid and dashed black curve in Fig. 7(a). Note that this 60%-reduction at n i ¼ 100 is caused by the factor multiplying γ e , which is proportional to (n i =N e ) 2=3 , i.e., by the small DOS of the electrons in InGaAs. The question arises how reliable this result is in view of the simplistic isotropic two-band dispersion model. Compared to silicon, the situation is more favorable here. The direct CB can be considered as isotropic, and the split-off hole band is more than 0.3 eV away from the VB edge. For Δ i h,nÀtype , screening is solely due to electrons. As holes are minorities, the Fermi level is close to or in the CB; hence, the valence bands (VBs) are just probed over an energy interval k B T which becomes narrow at low temperatures. This makes nonparabolicity of the hole bands negligible. The remaining error is then caused by disregarding their warping. But this only affects the pre-factor as can be explicitly seen from the last line in Eq. (B9) showing that Δ i h,nÀtype (β ! 1) is proportional to m h . Therefore, an error estimate can be obtained from the variation of m hh over the warped iso-energy surface. Citing a threeparameter warping model for GaAs (Ref. 44 and references therein), one finds (0:67 . . . 1:73) m hh , where m hh is the average over azimuthal and polar angles. The error bounds for Δ i h,nÀtype are hence similar. However, provided one could perform the fourfold angle integration numerically (together with the two remaining integrations), the angle dependence would be efficiently averaged out. Based on this, effects of the real band structure cannot be expected to nullify the main findings of this chapter: the very large value of Δ i h,nÀtype at low temperatures and its significant reduction due to nonparabolicity of the CB and increasing temperature.
The total BGN in quasi-neutral In 0:53 Ga 0:47 As is shown in Fig. 8. In the n-type material, it is completely dominated by the behavior of Δ i h,nÀtype . However, the nonparabolicity of the CB decreases the zero-temperature value to 237 meV at n i ¼ 100 (1:2 Â 10 19 cm À3 ). The corresponding value in the p-type material is 41 meV. At room temperature, the values are 69 and 45 meV, respectively. The dotted line indicates the published limit of electrical activation in the n-type material 41 (1:4 Â 10 19 cm À3 ). A similar value seems to hold for p-type In 0:53 Ga 0:47 As. 42 Hence, the above BGN values can be roughly considered as upper limits too.

V. CONCLUSION
The analysis of BGN in quasi-neutral regions of Si and In 0:53 Ga 0:47 As based on the finite-temperature RPA formalism including band nonparabolicity led to the following conclusions.
In n-type silicon, the gap shrinks by 119 meV at an active donor concentration of 1 Â 10 20 cm À3 . In p-type silicon, at the same density of active acceptors, the gap shrinks by 144 meV at 20 K and by 131 meV at 300 K, respectively. Both the difference in magnitude between doping types and the temperature dependence in the p-type material have the same origin-the ionic QPS of the minority carrier band (electrons) in the material with the smaller DOS (p-Si). A smaller DOS results in weaker screening.
In In 0:53 Ga 0:47 As, this effect becomes extraordinary at low temperatures. The 26 times smaller DOS in the parabolic model and the five to six times smaller mass of electrons make screening in the n-type material so weak that the ionic QPS of the minority carrier band (holes) takes giant values (450 meV at T ¼ 0 K at a donor density of 1:4 Â 10 19 cm À3 , i.e., more than half of the bandgap). However, this result is due to the assumption of a parabolic CB. The strong degeneracy enforces to include nonparabolicity which lowers the above value to 204 meV. A nonparabolic band effectively implies a larger DOS and a larger average mass to be used in a parabolic model, thus more screening. However, with increasing temperature, the ionic QPS of minority holes quickly shrinks to small values (29 meV at T ¼ 300 K) due to the T À1=2 -dependence in the high-T limit. For the BGN, this means a limit of 72 meV at room temperature. In p-type In 0:53 Ga 0:47 As, the DOS is large and screening is strong. The contribution from the ionic QPS of minority electrons is just 4 meV, and the total BGN amounts to 46 meV. This is accompanied by a very weak temperature dependence (43 meV at 20 K).
The nonparabolicity of the CB must be taken into account in all cases where the DOS/mass of holes is much larger than that of electrons, as often encountered in III-V materials. The paper has demonstrated that this is possible-numerically for all temperatures, and even analytically for the ionic QPS at T ¼ 0 K, if a simple isotropic two-band model is used, where the multi-valley CB and the heavy and light-hole bands are incorporated by DOS effective masses. There is hope that a fully analytical T ¼ 0 K-limit can also be found for the free-carrier correlation energies. Padé interpolation 1,29 then would enable a BGN model as a function of doping and temperature essentially free of material-specific fit parameters.
A serious problem is the impact of all the band structure effects, which are not considered in this work, like anisotropy, warping, contribution of split-off band, DOS tails, and polaronic mass enhancement. Even though a full numerical treatment based on an elaborate band structure model seems to be out of reach even in the rigid-shift approximation of BGN, some aspects could be tested on a smaller scale, e.g., in the screening function.

ACKNOWLEDGMENTS
I am obliged to Dr. Pietro Altermatt (Trina Solar) for making me aware of Refs. 12 and 38.

APPENDIX A: RPA DIELECTRIC FUNCTION FOR NONPARABOLIC BANDS
This appendix provides the derivation of the full-RPA dielectric function for nonparabolic bands, its static limit, as well as the high-T/low-density and low-T/high-density limits of the static form, respectively.

Frequency-dependent form at arbitrary density and temperature
With the Fourier transform of the bare Coulomb potential v(q) ¼ 8πa ex Ry ex =q 2 and the nonparabolic free dispersion (10), the RPA dielectric function (9) takes the form Energies are normalized by Ry ex , wave numbers by a À1 ex . The k z -axis is aligned with q, hence cos Θ ¼ k Á q=(jkkqj). Setting τ(t) ¼ 8γ a α a kqt with t ¼ cos Θ, one obtains The τ-integral can be solved using which gives Inserting into Eq. (A2) yields þ [h a (k, 0) À g ν,a ] ln [h a (k, 0) À h a (k, q) À g ν,a ]g, (A5) with the auxiliary functions Note that h a (Àk, q) ¼ h a (k, Àq) and h a (k, 0) ¼ h a (Àk, 0). In the case of parabolic bands, the k-integral in Eq. (A5) is integrated by parts to remove the logarithmic singularities. 1 This is impossible for nonparabolic bands at T . 0. Therefore, the final frequencydependent form for numerical integration is À g ν,a j jarctg g ν,a j j h a (k, 0) À h a (k, q)

Static form at arbitrary density and temperature
For the ionic part of the correlation energy, only the static limit g ν,a ! 0 of Eq. (A6) is needed. It can be written as where h a (κ, q) ¼ have to be used.
3. Debye limit and T = 0 K-limit of the static form The Debye limit 43 of the static RPA dielectric function follows from Eq. (A7) by q ! 0 and using Boltzmann statistics (hightemperature, long-wavelength limit), This is exactly the same as which is the expression in the case of parabolic bands. Therefore, the form of the Debye limit is independent of the band dispersion as expected.
To obtain the T ¼ 0 K-limit, the κ-integral in Eq. (A7) is integrated by parts where the derivative of the distribution function with respect to κ becomes a delta function Then, Only the majority carriers of the quasi-neutral regions contribute (ζ a . 0), whereas the minority carriers are frozen out. From (A12) the parabolic limit (γ a ! 0) is easily retrieved, (A13)

APPENDIX B: CORRELATION ENERGY OF ION-CARRIER INTERACTION
This appendix provides the derivation of the full-RPA correlation energy of ion-carrier interaction for nonparabolic bands in the high-T/low-density and low-T/high-density limits, respectively.

Arbitrary density and temperature
The correlation energy of ion-carrier interaction is calculated by Eq. (7). After normalization and integration over angles, it reads For the derivative of the dielectric function (A7) with respect to the Fermi energy, one has to apply @f a @ζ a ¼ À @f a @κ under the κ-integral.

Debye limit and T = 0 K-limit
In the high-T/low-density regime, the effect of nonparabolicity is negligible, and one obtains from Eq. (7) with the dielectric function (A10), the well-known result 1,29 In quasi-neutral regions n e þ n h ¼ n i , which simplifies Eq. (B3) to where the abbreviationζ a ¼ ζ a (1 þ γ a ζ a ) was introduced. For the derivative of the static dielectric function with respect to the Fermi energy @ϵ(q, 0)=@ζ a in the limit β ! 1, one needs where the static dielectric function ϵ T¼0 (q, 0) is defined in Eq. (A12) and the function h a in Eq. (A8).

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