Correlating electroluminescence characterization and physics-based models of InGaN/GaN LEDs: Pitfalls and open issues

The ongoing debate on efficiency droop in GaN-based light-emitting diodes (LEDs)^1–6 attests the present limitations of both characterization techniques and numerical models of III-nitride active optoelectronic devices. The inadequacies of the standard drift-diffusion (DD) simulation framework when used to describe carrier transport across a LED active region^7,8 have elicited the inclusion of semiempirical corrections for mechanisms such as tunneling,^9,10 carrier overflow,^11,12 ballistic overshoot,^13–15 “nonlocal” transport between quantum wells (QWs),¹⁶ and Auger-induced leakage.¹⁷ Unfortunately, the effects of those mechanisms are often hard to discriminate experimentally, which results in a problematic validation of the corresponding models. These intrinsic difficulties are compounded by at least three factors:

• the enduring lack of reliable first-principles models for nonradiative and radiative recombination mechanisms and other critical physical parameters (band offsets at heterojunctions,^18,19 interface polarization charges,^20,21 etc.) in realistic InGaN QWs,

• the inconsistencies in the definitions of some optical and transport parameters when used in experimental and in simulation contexts (e.g., “bulk” Auger coefficients vs. Auger rates in two-dimensional systems²²),

• the large (but often unacknowledged) sensitivity of several key parameters to technological/structural details (interface quality, layer homogeneity, doping and composition profiles, etc.) often unaccessible to an accurate experimental determination,

to which one could add the bias towards simpler, single-mechanism explanations against more comprehensive descriptions.^23–25 

We set out to illustrate the main issues faced when performing a combined experimental and simulation-based analysis of InGaN/GaN blue LEDs with an example carried out on a set of devices grown, characterized and simulated by the authors. Section II describes the test devices and presents the results of a systematic electroluminescence (EL) characterization for different temperatures, essential in the determination of the relative importance of competing loss mechanisms.²⁶ After several critical structure-related issues are addressed in Section III (contact geometry, interplay between current crowding, current leakage and vertical nonuniformity of carrier densities, etc.), Section IV applies two- and three-dimensional (2D, 3D) numerical device simulation to the determination of “optimal” recombination parameters allowing to reproduce the experimental L/I ratio as a function of current density and temperature. The results of this study, leading to an explanation of efficiency droop in terms of Auger-like nonradiative losses characterized by a C coefficient in the 10⁻³⁰ cm⁶/s range at room temperature, are critically discussed in Section V, where the uncertainties and limitations of the modeling framework are shown to correspond to a possible C overestimation by about an order of magnitude and, more generally, underscore the problematic identification between the Auger-like C coefficient considered by DD analysis and the actual Auger rates. As a conclusion, Section VI makes briefly the case for an improved microscopic model based on a more realistic description of the dynamics and distribution functions of carriers.

In order to help discriminate the recombination coefficients, our study was carried out on two sets of nominally identical LED structures grown on GaN templates differing for their threading dislocation densities (TDDs) N_disl. Low TDD (LDD) samples have N_disl ≈ 3 × 10⁸ cm⁻², while high TDD (HDD) devices have N_disl ≈ 8 × 10⁹ cm⁻². On top of the substrate, having a 100 μm nominal thickness, is a 2.6 μm-thick GaN buffer layer, with a n-type doping density 4 · 10¹⁸ cm⁻³. The active region, designed for 470 nm emission wavelength and non-intentionally-doped, has a total thickness of 64 nm and includes five InGaN QWs separated by GaN quantum barriers (QBs), with no electron or hole blocking layer (EBL/HBL). The absence of an EBL seems to have little impact on the onset of droop in the devices under study; a discussion is presented in Section III. The epitaxial stack is completed by a 130 nm-thick GaN cap layer with p-type doping level 3·10¹⁹ cm⁻³. From high resolution x-ray diffraction (HRXRD) spectra, the combined thickness of a QW/QB pair has been determined as 11.5 nm in both the LDD and HDD structures, while the In molar fraction in the QWs is in the range 0.15–0.20. As a limiting case, all simulations presented in Sections III and IV will assume in the active region an ideally piecewise constant composition profile with undoped 5 nm-thick In_0.15Ga_0.85N QWs and 6.5 nm-thick GaN QBs; the effects of deviations from these values will be briefly discussed in Section V C.

Through an etching process, the LEDs are patterned into square mesas, having width W = 0.5 mm and total nominal height H = 320 nm (about 130 nm of the n-type buffer have been etched). The p-type top layer is completely covered with a uniform, non-roughened semitransparent square contact of area S = W² = 0.25 mm². Finally, a 60 μm-wide cathode contact metallization is deposited on the partially exposed n-type GaN buffer, surrounding the entire mesa periphery at a 20 μm distance from the mesa sidewalls. A schematic cross-section of the LEDs under study and a detail of the mesh used in the 3D simulations are shown in Figs. 1 and 2.

The I(V) and L(I) characteristics (where V is the applied voltage, I is the current flowing through the anode and cathode contacts, L = η_extrP_op is the emitted optical power, P_op is the total optical power generated in the active region, η_extr is the extraction efficiency) and the EL spectra were measured under pulsed excitation over a controlled temperature range between 15 °C and 120 °C. In Fig. 3 are shown the experimental EL spectra for different temperatures in the LDD and HDD LEDs, respectively. Experimental results indicate that, predictably, LDD LEDs exhibit a significantly higher efficiency: the peak L/I ratio at room temperature in LDD LEDs is about twice as large as in HDD ones, see Fig. 4. LDD devices are also affected by a much smaller temperature sensitivity of the quantum efficiency; for an injection current I = 100 mA, corresponding to a current density J = 40 A/cm², a temperature increase ΔT = 105 °C leads to a reduction in the emitted optical power equal to 50% and just 10% in HDD and LDD LEDs, respectively (see Fig. 5). LDD LEDs show also a significant improvement in the stability of the peak emission wavelength λ_peak with respect to both temperature and injection current (see Fig. 6).

Our simulation analysis has been based on the APSYS modeling software from Crosslight Software Inc. and, to a lesser extent, on the TCAD Sentaurus suite from Synopsys Inc. APSYS and TCAD Sentaurus are well representative of the state of the art in optoelectronic computer-aided design (CAD) tools, whose limitations are the focus of our investigation. Both implement 2D and 3D DD models coupled with relevant quantum corrections, including the self-consistent solution of Poisson and Schrödinger equations in the active region. The valence band is described with a

$k·p$ Hamiltonian, allowing the determination of spontaneous emission as a function of wavelength and carrier density in strained QWs; many-body effects are also accounted for. Specific features of non-cubic semiconductors are taken into account, such as the spontaneous and piezoelectric polarization effects, and the description of carrier transport across the active region may be corrected with semiempirical models for e.g. tunneling and capture/escape processes between QWs and barriers. The nominal values of the material parameters used in the present simulations are reported in Appendix A 1.

Before addressing, in Section IV, the use of 3D simulation to extract “optimal” recombination coefficients allowing to reproduce the experimental L/I ratios as functions of current density, temperature and TDD, some preliminary considerations on our test devices are in order, relevant also to most LEDs discussed in the literature.

The standard ABC model used for correlating current and carrier density in the analysis of experimental data is:²⁷ 

where A, B, C are the “effective” Shockley-Read-Hall (SRH), radiative, and Auger recombination coefficients in the QWs, q is the electron charge, V_QW is the volume of all QWs, and I_leak accounts for all other loss mechanisms. This model is based on two main implicit assumptions about electron and hole densities n(x, y), p(x, y) in the QWs comprised into the active region: both n and p are essentially constant along x and y, and n ≈ p. These conditions should be satisfied also when more detailed models for the ABC coefficients are used, e.g. when their possible saturation at high carrier densities is taken into account^27,28 (as a general rule, the Auger recombination rate is proportional to n³ only if the Maxwell-Boltzmann carrier statistics can be used, but this is not usually the case in LEDs and lasers) or density-activated defect recombination (DADR) is considered.²⁹ Moreover, all leakage mechanisms, including carrier overflow (thermionic emission),^11,12 non-capture by the QWs or ballistic overshoot,^13–15 tunneling,^9,10 and Auger-induced leakage¹⁷ are usually either neglected or empirically approximated with n-dependent models. Some groups have also proposed ABCD models to describe leakage, and numerical estimates of a third (or higher) order coefficient corresponding to its effects have been presented e.g. in Refs. 30–33.

Contrary to the one-dimensional (1D) charge-control models based on Eq. (1), a 3D DD simulation framework can take into account selfconsistently carrier nonuniformity and leakage. Therefore, it may be of general interest to discuss in detail the impact of the structural features of the devices under examination on the calibration of the comparison between measured and simulated characteristics.³⁴ 

First, our LEDs, as well as most non-vertical structures, are affected by nonnegligible current crowding (see e.g. Refs. 35–41). The simulated lateral profile of the vertical component of the electron current density J_{n, y}(x) in the present devices at a current I = 100 mA (J = 40 A/cm²) is reported in Fig. 7, where a slightly stronger crowding may be observed for increasing temperature both in the LDD and HDD case. A map of the modulus of the total current density in the mesa at the same bias level is presented in Fig. 8, where the aspect ratio of the mesa W/H > 10³ is mirrored by the extensions of the considered x and y intervals. Not only the current in the vicinity of the mesa walls is significantly higher than the laterally averaged current reported in the L(J) and efficiency plots, but also carrier density levels in photogenerating regions may be significantly higher near the sidewalls than in the case of uniform lateral distribution at the same current levels, as shown in Fig. 9. Hence, Auger-like effects could be strongly enhanced^37,42–44 and the corresponding C coefficient could be overestimated from an analysis of experimental L/I data neglecting crowding effects.

Second, the active region of the devices under study is not bounded by an EBL, which is usually included in order to minimize electron leakage from the active region at high current levels. However, in the present structure, the onset of efficiency droop may be observed at current densities where leakage is probably still negligible. Fig. 10 shows that, with the selected values of material parameters, numerical simulations predict essentially no electron leakage from the active region into the p-cap at a current I = 100 mA where droop is already significant in low-TDD structures, over the entire temperature range 15–120 °C. (As discussed in Ref. 45, electron leakage may become less relevant for increasing temperature.) From the same figure, the importance of current crowding may be appreciated by comparing the J_{n, y}(y) profiles at different lateral coordinates.

Third, in addition to suppressing leakage, an EBL could help achieving a more uniform distribution of carriers and EL intensity among QWs, and hence a higher efficiency. An uneven vertical carrier distribution among QWs may have an effect similar to current crowding in enhancing Auger-like recombination mechanisms. In fact, lateral and vertical carrier nonuniformities translate collectively in a “reduced effective volume” of the active region, as discussed e.g. in Ref. 46. In order to check the possible effects of an EBL in the LEDs considered in the present study, a set of simulations has been performed on a structure including a 50 nm-thick Al_0.2Ga_0.8N layer with a 10¹⁸ cm⁻³ p doping density. (See e.g. Refs. 24,47 for discussions on the optimal choice of EBL thickness, composition and doping level.) Little impact from the EBL insertion has been observed in the simulations, up to a current I = 100 mA (see Fig. 11), either on leakage current and vertical carrier distribution or on current crowding.

As a final note in this discussion of carrier density distributions, we recall that realistic device modeling requires the inclusion of temperature-dependent incomplete dopant ionization. In the present study, activation energies E_D = 20 meV and E_A = 200 meV have been assumed in GaN for Si and Mg, respectively. These energies correspond to ionized donor and acceptor densities in the n-buffer and p-cap

$ND+≈1.7·1018,2·1018$ cm⁻³ and

$N_A^- \approx 2.7 \cdot 10^{17}, 10^{18}$

$NA−≈2.7·1017,1018$ cm⁻³ at T = 15 °C and T = 120 °C, respectively. An example of the importance of incomplete ionization is presented in Fig. 12, where the significant increase of the fraction of ionized acceptors at T = 120 °C contributes to more uniform distributions of carriers and radiative emission among the QWs of a LDD LED,⁴⁸ which partly compensate the adverse effects of temperature on IQE through increased nonradiative recombination coefficients and decreased radiative processes.

A seemingly trivial simulation detail is the appropriate assignement of the boundary conditions on the mesa walls. In the present 2D simulations, in order to halve the computational box, the device symmetry along the lateral direction (x) has been exploited by applying Neumann boundary conditions along a vertical axis at the center of the mesa, where the reference x = 0 has been taken (see Figs. 1 and 2). While the effect of this halving on the current flow across the device is obvious, the impact on the corresponding emitted optical power must be checked in order to avoid possible double counting, depending on the definition used by the CAD software for the power generated/emitted from spontaneous emission.⁴⁹ 

A less trivial symmetry issue is related to the square-ring shape of the cathode contact. Since the n-side metallization surrounds the entire mesa periphery, in the presence of current crowding a fully 3D simulation should be required.^7,38,50–55 Thanks to the device symmetries and lateral dimensions, in the following Sections a quasi-3D approximation has been adopted by determining temperature-dependent correction factors f_{J, ring}(T; V) = J_3D(T; V)/J_2D(T; V) and f_{L, ring}(T; V) = L_3D(T; V)/L_2D(T; V), used to weight the 2D-simulated current density and emitted optical power when compared with the corresponding experimental values. This approximation has been validated against fully-3D simulations performed with TCAD Sentaurus on increasingly realistic test structures (see Figs. 2 and 13). The severe current crowding predicted by both APSYS and TCAD Sentaurus for the LED structures under consideration lead to f_{J, ring}(T; V) ≈ f_{L, ring}(T; V) ≈ 2 (see Fig. 14). The use of uncorrected 2D models to fit experimental results in LEDs affected by current crowding may result in significant underestimation of the Auger-like recombination coefficients, as examined in the next Section. It is worth emphasizing that, however important, this modeling error can be easily prevented, and should not be confused with the possible sources of C overestimation mentioned in Section III A (mainly concerning 1D approximations) and the most critical ones discussed in Sections V B and V C.

Either with a charge-control model based on Eq. (1) or with a 3D numerical device simulator, it is straightforward to match a single experimental L/I curve of any given LED by adjusting the recombination parameters ABC. Unfortunately, such an exercise is of modest practical relevance, because of non insignificant concerns about the uniqueness of the “optimal” coefficients determined by the fitting process.

In the present study, the availability of experimental information not only on a wide range of injection currents and temperatures, but also on samples differing nominally only for their TDD (and hence for their SRH lifetimes) has allowed the definition of a large set of constraints, hopefully removing much of the arbitrariness from the optimization procedure leading to the appraisal of the recombination coefficients. As additional self-imposed restrictions, the coefficients C_n, C_p, used in the simulations to describe the Auger-like recombination rate as R_Auger = n²pC_n + np²C_p, were assumed to exhibit the approximately linear dependence on T observed in Refs. 67,68, while the A coefficient was assumed to be roughly proportional to N_disl⁶⁹ (see also the recent studies on the impact of TDD on LED efficiency presented in Refs. 63 and 70–76). Auger coefficients were not forced to be the same in LDD and HDD LEDs because, although nominally identical, secondary ion mass spectrometry (SIMS) measures suggest that the QWs of the LDD devices have probably a smoother profile, which could imply a lower Auger rate.⁷⁷ The radiative recombination coefficient B was not used in the present calculations because the radiative recombination rate R_rad is self-consistently determined by the simulator in the active region. The optimization space was effectively limited also by the absence of an EBL and by the seeming negligible contribution of leakage in the devices under study, discussed in Section III, which remove possible effects from e.g. uncertainties in the valence band offset at the barrier/EBL interface.⁴⁵ Some effects of model and structural parameter variability will be addressed in Section V.

The “optimal” SRH lifetimes and Auger-like coefficients, determined with the nominal values of all other material parameters listed in Appendix A 1, are reported in Table I along with previous estimates of the recombination coefficients for InGaN/GaN LEDs (see also Refs. 78,79). Both in the LDD and HDD case, the ratio between the present estimates of the C coefficients at 15 °C and 120 °C is approximately equal to the ratio between the corresponding absolute temperatures. If current crowding effects were neglected when comparing experiments and simulations, i.e. if f_{J, ring} ≈ f_{L, ring} ≈ 1 were assumed, C would be underestimated by a factor ≈1.5÷1.9 (see comment at the end of Section III A).

A comparison between measured and simulated IQE(J; T) with the “optimal” τ_SRH and C values is shown in Fig. 15. Fig. 16 reports the experimental and simulated EL spectra for different temperatures and currents in a LDD and HDD LED, respectively; for recent comparable experimental studies on the dependence of λ_peak on current and temperature, see e.g. Refs. 80–82. It is worth mentioning that, for a realistic description of radiative emission, it has been essential here to take into account the broadening of interband transitions due to inhomogeneities in QW thickness and composition.^83–85 This can be appreciated from Fig. 17, which compares simulated EL spectra at different temperatures with and without broadening in a LDD LED. The spectra are normalized, as the adopted model for broadening does not affect the total generated optical power or the IQE. The inclusion of many-body Coulomb interactions^86–92 as implemented in the present versions of the simulators has not resulted, with the exception of E_g renormalization, in a real improvement in their predictive capabilities, due to the large set of widely-adjustable empirical parameters involved in the corresponding models.

While a systematic uncertainty quantification is beyond the reach of the present analysis of GaN-based LEDs, we will attempt here a critical discussion of the results derived in Section IV. Our focus will be on what we consider the weakest points of the DD modeling framework, whose combined limitations could determine, in particular, a severe overestimation of the Auger-like coefficient also beyond the present case study.

In the previous analysis, a standard DD model for carrier transport across the active region has been adopted, not including carrier flyover (ballistic transport)^13–15,24 or Auger-induced leakage.¹⁷ The dependence of carrier velocity on the electric field has been approximated with a Canali-like description⁹³ using state-of-the-art model parameters.^94–96 We have also experimented modified transferred-electron models⁹⁷ as well as hydrodynamic and energy-balance descriptions,⁹⁸ without noticeable improvements in the predictive capabilities of the simulator.

Among “quantum–corrected” transport models, a well-established approach available in most commercial simulators involves the use of separate continuity equations for carriers in unconfined bulk states (3D) and in bound QW states (2D).⁹⁹ This split-state description, separating the 3D and 2D quasi-Fermi levels in the active region, requires the definition of an appropriate recombination-like coupling term, in turn involving a set of time constants τ_QW associated with the capture/escape processes.^100,101 A recent body of work using split-state models in InGaN-based LEDs^{16,24,102,103} suggests that capture time constants τ_QW ≈ 10⁻⁶ s are required in order to obtain realistic IQE values if other loss mechanisms (Auger recombination in particular) are assumed to be negligible. The need of such long τ_QW has led the authors of Ref. 24 to rule out incomplete QW capture (at least, in its present formulation) as a possible droop mechanism. A similar conclusion seems appropriate also in the case of the present LED structures, since not only low capture rates have to be assumed (τ_QW ≈ 10⁻⁷ s), but also the IQE dependence on current density and temperature cannot be satisfactorily reproduced when long τ_QW are used in the simulations (see Fig. 18).

In the previous Sections, spontaneous emission has been calculated with the present APSYS default parameters of the valence band

$k·p$ Hamiltonian, reported in Table II and derived from Ref. 105 for GaN and Ref. 106 for InN. However widely adopted, those parameter values are affected by a considerable uncertainty, which directly affects the LED characteristics. As a demonstration of the importance for device simulation of the description of the valence bands near Γ, Fig. 19 compares the computed IQE and EL spectra corresponding to the default band parameters and to a more recent

$\bm {k} \! \cdot \! \bm {p}$

$k·p$ model derived from many-body perturbation theory in the G₀W₀ approximation.¹⁰⁴ The newer

$\bm {k} \! \cdot \! \bm {p}$

$k·p$ model, in particular its larger crystal-field splitting of GaN, leads to a nonnegligible decrease in both the EL intensity and the IQE at all injection currents; if used in Section IV, it would have led to significantly lower “optimal” estimates of both the A and C coefficients.

Fig. 20, reporting the simulated IQE as a function of the injection current for different values of the C Auger-like coefficient and of the TDD-SRH-like lifetime, could be used to assert our confidence in the nonradiative recombination coefficients estimated from the analysis presented in Section IV.

The order of magnitude of our A and C values is indeed in general agreement with several simulation studies discussing Auger effects on LED droop^{3,4,23,24,47,107} and recent experimental studies on Auger recombination, focusing on its temperature¹⁰⁸ and composition^65,66 dependence and reaffirming the importance of Auger processes in the blue spectral range. Notably, the Auger-like coefficient determined in the present LED study is increasing with temperature, i.e. it would be compatible with Auger processes in InGaN/GaN QWs rather than with density-activated defect recombination or leakage from the active region, both decreasing with increasing temperature.^26,29,45

However, in our opinion, an accurate determination of the actual Auger rates remains elusive to a combined experimental and simulation analysis based on state-of-the-art DD-based optoelectronic CAD tools, not only because the applicability of the ABC model to QWs is somehow questionable,²⁷ but also on account of possible enhancements of Auger recombination by Coulomb interactions¹⁰⁹ and of the significant uncertainties on the radiative losses (see Section V B and Ref. 92), which are also strongly affected by the microscopic details of the MQW heterostructures. If the relationship between recombination coefficients in bulk and in QWs is not better explored, beyond the attempt presented in Ref. 110, the very “portability” of C between different devices is dubious.²² Direct calculations of Auger recombination in QW structures, usually based on

$k·p$ models, agree in suggesting the importance of QW composition profiles^77,111 and polarization fields¹¹² on the direct process rates. Moreover, fluctuations of the In molar fraction across the active region^113,114 could have significant effects on both leakage and Auger recombination.^115,116

Even assuming an ideally uniform QW profile across the device, as we have done in the previous Sections, the determination of C is significantly influenced by the uncertainties on QW thickness and composition. As an example, Fig. 21 compares the IQE corresponding to the nonradiative recombination coefficients estimated in the previous Section (with the “worst case” 5 nm-thick In_0.15Ga_0.85N QWs) and the results obtained assuming thinner (2 nm) and deeper (up to 0.20 In molar fraction) QWs. Electron leakage and current crowding are minimally affected by these changes, while the use of the worst-case-geometry A and C coefficients with the thinner QWs leads to a significantly larger droop, and the C estimate could be reduced approximately by a factor 3 with 2 nm-thick In_0.2Ga_0.8N QWs.

In addition to these uncertainties, recent modeling work¹⁷ suggests that experimentally-derived C values could overestimate the real ones up to a factor 2 since they do not take into account Auger-induced leakage effects. Finally, possible non-Auger contributions to the Auger-like coefficients could derive from photon quenching mechanisms¹¹⁷ or from strong nolinearities in the n-dependence of SRH-like recombination processes determined by multiple-level defects,⁷⁶ and would further hamper the discrimination between Auger and other nonradiative mechanisms from an analysis of the LED characteristics.

After a discussion of some important pitfalls stemming from 1D and 2D simulation of LED structures – foremost the possible incorrect determination of the C coefficient in the presence of current crowding – we have presented a combined experimental and simulative study on InGaN/GaN LEDs where a seeming excellent match between measured and calculated IQE characteristics and EL spectra has been achieved over a wide range of temperatures and operating currents, by using as the only fitting parameters the A and C coefficients conventionally associated with SRH-like and Auger recombination.

However, even if the present study is based on a large and consistent set of experimental data, we have shown that the estimates of A and especially C are affected by large error margins because of the many sources of uncertainty, related not only to device structure and technological details (distribution, activation and ionization of dopant species; composition profiles of QWs and barriers in the active region; uniformity of contacts; parasitic resistances and other interface effects; etc.) and to material parameters (band offsets and polarization charges at heterojunctions, energy gap and electronic structure of GaN and InGaN in presence of strain, defects, nonideal interfaces, etc.), but also to a possibly inadequate description of microscopic processes in the active region (radiative recombination, intra- and inter-subband carrier scattering, capture/escape between bound and unbound states, carrier flyover and tunneling etc.) and of their problematic inclusion in a device-level DD CAD suite.

These uncertainties are compounded by the nontrivial relationship between the bulk C coefficient used in the simulations and the corresponding 2D Auger rate, and by the possible role played by non-Auger recombination mechanisms with a strong nonlinear dependence on the carrier density; a calculation approach for QW structures “beyond

$k·p$” is needed^68,118,119 also in order to take advantage of possible polarization-engineering opportunities.^112,120,121

The prime requisite for a better understanding of experimental information is probably an improved description of carrier transport across the active region. Promising efforts are represented by momentum-resolved models^23,122 and analytic-band Monte Carlo (MC) transport simulations,¹²³ but the need for a full-band MC approach (FBMC, possibly integrated with microscopic models of radiative and nonradiative recombination processes) is apparent, especially after the recent debate following the proposed identification of Auger recombination as the dominant mechanism in droop from electron emission spectroscopy experiments,¹²⁴ where FBMC analysis has suggested that the LED structure under investigation was probably unsuitable to recover an Auger signature.¹²⁵ 

The authors thank Simone Chiaria, Andrew Phillips and Michel Lestrade for useful discussions, Simone Vaccari and Marco Ferretti for their contribution to the experimental characterization. The work at Politecnico di Torino and Boston University was supported in part by the U.S. Army Research Laboratory through the Collaborative Research Alliance (CRA) for MultiScale multidisciplinary Modeling of Electronic materials (MSME).

The nominal values of the material parameters used in the present simulations are listed in Table II (see also Refs. 47,106,126,127), where a and c are the lattice constants, E_g is the energy gap, α and β are the coefficients of the Varshni model for the temperature dependence of E_g,¹⁰⁶ χ is the electron affinity, Δ_so and Δ_cr are the spin-orbit and crystal-field splitting, ɛ₀ is the relative dielectric constant,

$me∥$ and

$m_{e}^\perp$

$me⊥$ are the longitudinal and transverse components of the electron effective mass, A₁…A₆ are the Luttinger parameters of the

$\bm {k} \! \cdot \! \bm {p}$

$k·p$ Hamiltonian, D₁…D₆ are the valence band deformation potentials, c₁₁…c₄₄ are the elastic constants for wurtzite materials, P_SP is the spontaneous polarization, μ_{n, p} are the electron and hole mobilities (assumed to be independent of temperature, ionized impurity densities, and In molar fraction). A linear interpolation between the values of the binary compounds has been used for all In_xGa_1−xN parameters, with the exception of E_g and P_SP, where a quadratic correction x(1 − x)b has been included, with bowing factors b equal to 3.0 eV and -0.038 C m⁻², respectively. A conduction band offset ΔE_c/ΔE_g = 0.67 has been assumed for all heterojunctions, regardless of the In molar fraction.