Modeling of temperature and excitation dependences of efficiency in an InGaN light-emitting diode

The changes in excitation dependence of efficiency with temperature is modeled for a wurtzite InGaN light-emitting diode. The model incorporates bandstructure changes with carrier density arising from screening of quantum-confined Stark effect. Bandstructure is computed by solving Poisson and k.p equations in the envelop approximation. The information is used in a dynamical model for populations in momentum-resolved electron and hole states. Application of the approach shows the interplay of quantum-well and barrier emissions giving rise to shape changes in efficiency versus current density with changing temperature, as observed in some experiments.


Introduction
Considerable resources are being devoted to advancing InGaN light-emitting diodes (LEDs), largely because of solid-state lighting. Of particular interest is improving device efficiency [1], especially the mitigation of efficiency degradation with increasing excitation, i.e., efficiency droop [2,3,4]. To arrive at a solution, it is important to understand the principal mechanism of emission together with the effects of carrier leakage [5], Auger recombination [2], junction heating [6], carrier and defect delocalizations [7,8,9]. Equally important is the change in efficiency with temperature, not only for application but also for better understanding of the physics contributing to efficiency droop. 1 This paper investigates LED efficiency as functions of current density and lattice temperature. The analysis uses a model that allows direct input of bandstructure properties [10,11]. Band-structure details are important, because underlying emission properties in wurtzite quantum-well (QW) structures are the changes in band energy dispersions, confinement energies and optical transition matrix elements with excitation. These changes arise from screening of piezoelectric and spontaneous polarization fields [12,13,14].
The present approach differs substantially from the total-carrier-density rateequation description used in most discussions involving InGaN LED efficiency. In terms of reproducing experimental efficiency versus injection current data, a particularly successful rate-equation model is the ABC model [2,15]. The model's name derives from the phenomenological constants A, B and C, introduced to account for Shockley-Read-Hall (SRH), radiative-recombination and Auger-scattering carrier losses, respectively. Band-structure properties enter indirectly via these coefficients.
Section 2 summarizes the model and presents the working equations. Section 3 demonstrates the application of the model by calculating internal quantum efficiency (IQE) as a function of injection current for different lattice temperatures. Following Ref. [16], the calculations are performed for a LED with an In 0.37 Ga 0.63 N QW active region. The choices of input parameters are explained and comparison to recent experimental results is discussed. This section also describes the contributions from Auger carrier loss, plasma heating and carrier leakage. Then, the role of the band structure is illustrated by considering a LED with an In 0.20 Ga 0.80 N active region. Section 4 summarizes the paper.

InGaN LED model
The derivation of spontaneous emission from QW and barrier transitions starts with a Hamiltonian adapted from quantum optics [17]. With this Hamiltonian, Hiesenberg operator equations of motion are derived for the carrier populations and polarizations. These equations contained operator products describing coupling to higher order correlations in light-matter interaction. Truncating at the Hartree-Fock level, as well as adding phenomenologically SRH carrier loss and relaxation contributions from carrier-carrier and carrier-phonon scattering, lead to the working equations for the present investigation [11].
For the QW population n σ,ασ,k ⊥ , where each QW state is denoted by the subscripts σ, α σ , k ⊥ , representing charge (σ = e or h), subband index (α σ ) and in-plane momentum (k ⊥ ). In the above equation, A is the SRH coefficient, γ c−c and γ c−p are the effective carrier-carrier and carrier-phonon collision rates, respectively. Carrier loss by Auger scattering is modeled with a term containing an effective Auger scattering rate γ ag = CN 2 , where C is the Auger coefficient and N is the total (QW and barrier) carrier density. In the summation, σ ′ = h or e, when σ = e or h. For the barrier carrier population n b σ,k , where each barrier state is denoted by the subscripts σ, k, representing charge (σ) and 3-dimensional momentum (k). In the above equation, J is the current density, e is the electron charge and N p σ = k f ε b σ,k , µ p σ , T p . The injected carrier distribution f ε b σ,k , µ p σ , T p is a Fermi-Dirac function with chemical potential µ p σ , temperature T p . In the above equations, ε σ,k ⊥ and ε b σ,k are the carrier energies in the QW and barrier, respectively. From the derivation, the spontaneous emission coefficients are given by where ℘ ασ ,α σ ′ ,k ⊥ and Ω ασ ,α σ ′ ,k ⊥ are the QW dipole matrix element and transition energy, ℘ k and Ω k are the barrier dipole matrix element and transition energy, and ǫ b is the host permittivity. For the asymptotic Fermi-Dirac distributions approached via carrier-carrier collisions f (ε σ,k ⊥ , µ σ , T ) and f ε b σ,k , µ σ , T , the chemical potential µ σ and plasma temperature T are determined by conservation of carrier density and energy. For f ε σ,k ⊥ , µ L σ , T L and f ε b σ,k , µ L σ , T L , which are reached via carrierphonon collisions, the chemical potential µ L σ is determined by conservation of carrier density and the lattice temperature T L is an input quantity. For f (ε σ,k ⊥ , µ ag , T ag ) and f ε b σ,k , µ ag , T ag , which are reached via Auger scattering, there is a loss of carrier density that eventually leads to a chemical potential µ ag of the intrinsic semiconductor, i.e., approximately at mid gap, with its exact location depending on electron and hole dispersions. T ag is determined by noting that energy is conserved in an Auger scattering event.
When summed over all states, Eqs. (1) plus (2) reduce to total carrier density rate equations used in the ABC model, if we assume that the total emission contribution equals BN 2 , where B is an effective bimolecular radiative coefficient. In the summation over states, quantities such as the total carrier density and energy are computed by converting the sum over states to integrals, i.e., where S and h are the surface area and thickness of the active active region consisting of the QW and barrier layers. The conversion to bulk (3-d) density is via division by the QW layer width h in the case of the QW and by the total barrier width h b in the case of the barrier. Further discussions involving implementation and comparison with results from quantum-kinetic calculations may be found in earlier reports [10,11,18]. The dynamical solution gives the carrier densities in QW and barrier states. Band-structure information enters directly into Eqs (1) and (2) via the dipole matrix elements ℘ ασ,α σ ′ ,k ⊥ , ℘ k and carrier energies ε σ,k ⊥ , ε b σ,k . Iterative solution of the k · p and Poisson equations [19] gives the QW energies and dipole matrix elements. For this calculation, the necessary bulk wurtzite material parameters are given in Refs. [20,21,22,23]. When performing the band-structure calculation, quasiequilibrium condition is assumed to determine the QW and barrier bulk densities used in the solution of Poisson equation. This is an inconsistency that is acceptable provided the dynamical solution does not produce carrier distributions deviating too far from quasiequilibrium distributions. Application of the model involves first calculating the band structure for a range of carrier densities by solving the k · p and Poisson equations. This gives the band-structure information needed for the solving the population dynamical equations. To facilitate the numerics, carrier states are grouped into two categories: those belonging to the QWs and those belonging to the barriers. With a high internal electric field, the distinction between QW and barrier states may be ambiguous. In this paper, the choice is made by calculating QW dz |u mσ,βσ (z)| 2 , where integral is performed over the QWs and u mσ,βσ is the wavefunction. The states where the integral is greater than 1/2 are grouped as QW states and the rest as barrier states. For the problem being addressed, which are the excitation and temperature dependences of IQE, the distinction is only important because only QW transitions are affected by the quantum-confined Stark effect (QCSE). For the barrier transitions, the dipole matrix element in the presence of an internal electric field is approximated by an average, where each transition is weighted according to the occupations of the participating states. Grouping the barrier states appreciably reduces numerical demand, which remains substantial because one is still keeping track of a large number of k-states.
In the numerical solution of Eqs. (1) and (2), the band-structure quantities are updated at each time step according to the instantaneous carrier density. When steady state is reached, IQE is determined from dividing the rate of carrier (electron or hole) loss via spontaneous emission by the rate of carrier injection, i.e.
3 Results The above model was applied to compute the temperature dependence of IQE versus current density. A motivation is that adding temperature dependence may produce further physical insight to efficiency droop, as well as provide more stringent testing of the model. In order to relate to experiments [24,25,16], simulations were performed for two LED configurations. The experimental devices had single-QW active regions, which circumvented complications arising from nonuniform population in multi-QW structures. One device consisted a 2nm In 0.37 Ga 0.63 N QW between GaN barriers [24], while the other device consisted a 3nm In 0.20 Ga 0.80 N QW between GaN barriers [25]. The experiments [24,25,9] show features of IQE dependences on excitation and temperature, that if understood, may shed light on LED physics in general and the efficiency droop in particular. They involve changes in shape of IQE versus current density curves with temperature and the relationships among these curves after the onset of droop. In the former, for lattice temperatures 300K and above, both In 0.37 Ga 0.63 N and In 0.20 Ga 0.80 N LEDs exhibit IQE versus current density behavior showing the familiar efficiency droop, as described relatively well by the ABC model. Interestingly, a second IQE bump appears at lower temperatures for the In 0.37 Ga 0.63 N sample and becomes very pronounced at 100K. The ABC model is unable to describe this shape change, which does not occur in the In 0.20 Ga 0.80 N sample at any temperature. A goal of this study is to use our model to reproduce and understand the origin of the second bump and determine whether the different behavior between the two LEDs is intrinsic, i.e., involving only differences in band structure. Additionally, the model will be used to examine the detailed relationship among IQE curves at high current densities. In the In 0.37 Ga 0.63 N case, the curves for different temperatures are clearly separated for the entire measured current density range. On the other hand, for In 0.20 Ga 0.80 N, the IQE versus current curves for different temperatures either cross or merge with each other at high current densities. Previous modeling of these behaviors has lead to puzzling results, such as inferring a decrease in C coefficient with increasing temperature, contrary to what is expected of Auger scattering [9]. Figure 1 shows results from the present model for the In 0.37 Ga 0.63 N LED. The curves describe IQE versus current density at different lattice temperatures. Similar to experiment, there is a double bump excitation dependence in IQE. The lower excitation bump decreases with increasing temperature and vanishes for T L > 250K. Input parameters that are assumed temperature dependent are the SRH and Auger coefficients, A and C, respectively and the carrier-phonon scattering rate γ c−p . They are adjusted to produce IQE behavior resembling those found in Refs. [24] and [16]. The SRH coefficient is decreased with decreasing temperature to give the onset of emission as depicted in the figure. The Auger coefficient is adjusted to reproduce approximately the experimental maximum IQE values. Since Auger carrier loss conserves energy, significant plasma heating may occur. The role of carrier-phonon scattering is to relax the hot carrier distributions back to the lattice temperature. Both C and γ c−p affect the peak IQE. However, the current density where IQE peaks is relatively sensitive to γ c−p , which is chosen to produce an onset of droop in the tens of A/cm 2 range. Figure 2 plots the values of A, C and γ c−p versus temperature used in computing the curves in Fig. 1. A decrease in A with decreasing temperature is expected for defect related loss. Figures 2 (b) and 2 (c) show increases in Auger coefficient and carrier-phonon scattering rate with increasing temperature, which are consistent with microscopic calculations. The values of Auger coefficient are within the range predicted for phonon-assisted Auger scattering [4] and smaller than most values obtained from experimental curve fitting with the ABC model [15]. The values used for carrier-phonon scattering are higher than predicted by quantum kinetic calculations for typical III-N structures. However, it should be noted that those calculations are for intraband scattering among nearby states. The present effective rates represent the relaxation of very energetic states populated by Auger scattering to the QW ground state. The energy differences are in the neighborhood of the bandgap energy (2.7eV ).
To give some insight into the origin of the double bump, Fig. 3 (a) shows the spontaneous emission contributions from QW and barrier states versus current density for the T L = 200K curve in Fig. 1. Clearly, the spontaneous emission from barrier states dominates the emission from the QW at low current density, coinciding with the location of first IQE bump [see left arrows in Figs. 1 and 3 (a)]. Here, optical emission from barrier transitions occurs via the contribution k b k n b e,k n b h,k , as soon as the product of electron and hole populations, n b e,k n b h,k becomes nonzero. However, there may be no optical emission from QW transitions because of absence of localized QW states from the strong piezoelectric and spontaneous polarization fields (quantum-confined Stark effect). There is a transition region where the screening of internal electric field results in the appearance of localized QW states. Here, the small energy separation between QW and barrier states gives comparable QW and barrier populations, especially for the holes. However, the QW emission contribution αe,α h ,k ⊥ b αe,α h ,k ⊥ n e,αe,k ⊥ n h,α h ,k ⊥ can still be negligible, even though the product n e,αe,k⊥ n h,α h ,k⊥ may be appreciable, because QCSE spatially separates electrons and holes in the QWs, resulting in very small dipole matrix elements for QW transitions. At high current density [see right arrows in Fig. 1  and 3 (a)], plasma screening restores the QW localized states and optical-dipole matrix element, leading to QW emission overtaking that from the barrier because of the advantage in 2-dimensional versus 3-dimensional carrier density of states.
The above explanation is supported by Figs. 3 (b) and 3 (c). The solid curves are absolute square of envelope functions for electrons and holes at current densities indicated by arrows in Fig. 3 (a). Each curve is displaced according to its energy for clarity. Envelope functions belonging to QW and barrier states are indicated by black and grey curves, respectively. The former is clearly missing in Fig. 3 (b). The black dashed lines plot the confinement potentials and they indicate a significantly stronger deformation in Fig. 3 (b) compared to Fig. 3 (c) due to the higher internal electrical fields in the lower current density situation. Figure 4 plots the contributions to the T L = 200K IQE curve for various recombination and scattering processes. The solid red curve is the sum of QW and barrier emission, with a slight change in slope at low current density indicating the transition from predominately barrier to predominately QW emission. A second slope change at higher current density arises from the onset of Auger carrier loss. SRH and leakage (or failure to capture) losses from the barrier is depicted by the blue dashed curve. They are primarily responsible for the low IQE at low current densities. The green dotted curve shows that defect (SRH) loss from QW states are essentially negligible. Auger carrier loss accounts for the difference between the solid black curve and the sum of all the other curves.
To illustrate the Auger contribution, Fig. 5 (a) shows the IQE versus current density for two Auger coefficients: C = 2.3 × 10 −31 cm 6 s −1 , which is used in Fig. 1, and C = 10 −34 cm 6 s −1 , which is in the neighborhood of what is typically expected for materials with bandgap energy of approximately 2.7eV. The dashed curve clearly shows that for C less than 10 −34 cm 6 s −1 , Auger carrier loss has negligible effect on IQE. The difference between the two curves indicates the Auger contribution to the IQE curve in Fig. 1 for T L = 200K. During each Auger event, an electron and a hole recombine, with the loss in energy transferred to a second electron or hole. In a wide bandgap material, such as InGaN, the transferred energy is significant and can lead to highly-energetic carrier distributions, depending on the carrier-phonon scattering rate. Because of rapid carrier-carrier scattering, this change in carrier distributions may be described in terms of the plasma temperature T . Figure 5 (b) shows the plasma temperature versus current density for the two Auger coefficients. The solid curve indicates significant rise in temperature for the T L = 200K curve in Fig.  1. The dashed curve deoicts the rise in temperature in the absence of Auger loss, so that heating comes mainly from the relaxation of carriers from the barrier states, where they are injected, to the ground states, where emission occurs.
The role of band structure is investigated by performing simulations for the In 0.20 Ga 0.80 N LED, using the same parameter values as plotted in Fig.  2. Figure 6 shows the IQE versus current density for different temperatures. Clearly indicated is a missing or negligible second IQE bump. Band-structure calculations traced the reason to a smaller internal electric field with lower In concentration in the QW. Consequently, localized QW states exist for the entire temperature and current-density range considered, and their emission exceeds that of the barrier states at all current densities. The absence of a transition from barrier to QW emission dominance gives rise to the lack of a double bump shape.
Experimental evidence supporting the above argument may be found in [25].
There, the IQE versus current density curves shows single-bump behavior for T L between 100K and 300K, as in Fig. 6. However, agreement is not complete for the simulations fail to exactly reproduce experimental observation at high current density, where the experiment shows crossing of the IQE curves for different temperatures. At high excitation, the fitting parameters used in this study which influence LED behavior are C and γ c−p . We had assumed that their values do not change considerably with a change in Indium concentration from Ref. [25], we again consider A, C and γ c−p as free-parameters. Figure 7 shows the IQE results for input parameters given in Fig. 8. Figure 7 (a) indicates better agreement with with Ref. 14. However, agreement between theory and experiment for the In 0.37 Ga 0.64 N LED with the new set of parameters is degraded. On the other hand the single-versus double-bump behavior remains consistent with experiment and appears robust to changes in the fitting parameters.

Summary
This paper uses an approach to modeling InGaN LEDs that involves the selfconsistent solution of band structure and carrier dynamics. One motivation is systematic incorporation into a model, at a microscopic level and on equal footing, the effects of spontaneous emission, carrier capture and leakage, as well as nonequilibrium effects such as plasma heating. Another motivation is to provide direct input of band-structure properties, in particular, their carrier-density dependences arising from screening of piezoelectric and spontaneous polarization fields. Doing so allows comparison of different devices at the heterostructuredesign level.
The end result is a non-equilibrium microscopic model that provides a more precise estimation of relative strengths of possible physical contributions compared to the commonly used ABC model [2] and is easier to implement than a first-principles, many-body approach [3]. The model reproduces the main experimental features of the observed droop behavior versus temperature relatively successfully. Shape changes of IQE versus current density with temperature are described with the optical emission from quantum-well and barrier transitions treated on equal footing. The IQE droop is described by Auger scattering and resulting plasma heating, with material parameters and their temperature dependences that are consistent with microscopic calculations.
However, there are remaining discrepancies. For example, while the experimental and theoretical maximum IQEs can be made to roughly agree, the current densities where they occur are overestimated by the model. Several mechanisms may be responsible, e.g., effects arising from doping profile, presence of carrier blocking layers and interface irregularities [5,27]. These effects are ignored in the present study, but may be readily incorporated into the model. Also excluded is a more detailed description of defect recombination, e.g. as proposed for describing the experimental data [9,16]. More serious may is the neglect of many-body effects, such as treated in Refs [3] and [14].
There is also room for refinement in treating the connection between the bandstructure and population-dynamics aspects. Carrying out these improvements increases computational demands, and their implementation can be greatly facilitated by having more experimental data and better knowledge details on the experimental configurations.    Fig. 1. Solid curves in (b) and (c) are absolute square of envelope functions at zone center (k ⊥ = k = 0) for electrons and holes at current densities indicated by arrows in Fig. 3 (a). The x-axis is along the growth direction. Each curve is displaced according to its energy for clarity. Envelope functions belonging to QW and barrier states are indicated by black and grey curves, respectively. The black dashed lines plot the confinement potentials.