Localization landscape theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes

TL;DR

This paper presents a new computational approach using localization landscape theory to incorporate quantum disorder effects into classical models, enabling efficient and accurate simulation of carrier transport and recombination in LEDs with disordered heterostructures.

Contribution

The study introduces a self-consistent, faster method to include quantum effects in drift-diffusion models for disordered semiconductors, validated against experimental LED data.

Findings

Accurately models current-voltage characteristics of nitride-based LEDs

Provides a good approximation to the density of states in disordered systems

Enables detailed analysis of quantum effects in carrier transport

Abstract

This paper introduces a novel method to account for quantum disorder effects into the classical drift-diffusion model of semiconductor transport through the localization landscape theory. Quantum confinement and quantum tunneling in the disordered system change dramatically the energy barriers acting on the perpendicular transport of heterostructures. In addition they lead to percolative transport through paths of minimal energy in the 2D landscape of disordered energies of multiple 2D quantum wells. This model solves the carrier dynamics with quantum effects self-consistently and provides a computationally much faster solver when compared with the Schr\"odinger equation resolution. The theory also provides a good approximation to the density of states for the disordered system over the full range of energies required to account for transport at room-temperature. The current-voltage characteristics modeled by 3-D simulation of a full nitride-based light-emitting diode (LED) structure with compositional material fluctuations closely match the experimental behavior of high quality blue LEDs. The model allows also a fine analysis of the quantum effects involved in carrier transport through such complex heterostructures. Finally, details of carrier population and recombination in the different quantum wells are given.

Tables

Table 1. Table 1: Band structure parameters for wurtzite nitride alloys: Vurgaftman  et al. ( 2001 ); Piprek ( 2007 ) bandgap, relative permittivity, and effective masses.

	$E_g$	${\epsilon }_r$	$m_e^{\parallel }$	$m_e^{⟂}$	$m_{hh}$	$m_{lh}$
units	(eV)		($m_0$)	($m_0$)	($m_0$)	($m_0$)
GaN	3.437	10.4	0.21	0.20	1.87	0.14
InN	0.61	15.3	0.07	0.07	1.61	0.11
AlN	6.0	10.31	0.32	0.30	2.68	0.26
Bandgap alloy				InGaN: 1.4
bowing parameter				AlGaN: 0.8

Table 2. Table 2: Piezoelectric coefficients for wurtzite III-N materials. Romanov  et al. ( 2006 )

	$e_{33}$	$e_{31}$	$e_{15}$
units	(C/cm2)	(C/cm2)	(C/cm2)
GaN	0.73	-0.49	-0.40
InN	0.73	-0.49	-0.40
AlN	1.55	-0.58	-0.48

Table 3. Table 3: Parameters of polarization values. Fiorentini  et al. ( 2002 ); Ambacher  et al. ( 2016 )

	$a$	$b$	$c$
units	(C/cm2)	(C/cm2)	(C/cm2)
InxGa1-xN	-0.042	-0.034	0.037
AlxGa1-xN	-0.090	-0.034	0.021

Table 4. Table 4: Simulation parameters of each epi-layer. David and Grundmann ( 2010a ); Li  et al. ( 2014 ); Kumakura  et al. ( 2005 )

	n-GaN	i-InGaN	p-AlGaN	p-GaN
thickness (nm)	200	67	20	100
${\mu }_e$ (cm2/Vs)	200	300	100	32
${\mu }_h$ (cm2/Vs)	23	10	5	5
doping (cm-3)	$5\times {10}^{18}$	${10}^{17}$	$3\times {10}^{19}$	$2\times {10}^{19}$
$E_a$ (meV)	25	NA	215	170
${\tau }_n^{nonrad}$(s)	10	${10}^{-7}$	${10}^{-7}$	$6\times {10}^{-10}$
${\tau }_p^{nonrad}$(s)	$7\times {10}^{-10}$	${10}^{-7}$	10	10
$B_0$ (cm3/s)	$3\times {10}^{-11}$	$3\times {10}^{-11}$	$3\times {10}^{-11}$	$3\times {10}^{-11}$
$C_0$ (cm6/s)	$2\times {10}^{-31}$	$2\times {10}^{-31}$	$2\times {10}^{-31}$	$2\times {10}^{-31}$

Table 5. Table 5: Computation time required for each iteration step when solving the Poisson, DD, landscape, and Schrödinger equations for a given number of nodes, as tested by our home-built software and compared with other disorder models in nitrides. Watson-Parris  et al. ( 2011 ); Schulz  et al. ( 2015 ); Auf der Maur  et al. ( 2016 ); Watson-Parris ( 2011 )

	Node number	Computation time (s)
	(matrix size)
Poisson	428,655	25
Drift-diffusion	428,655	50
Localization landscape	428,655	50
Schrödinger	428,655	63,650
Ref. 14, 62	1,500,000	60,000
Ref. 20	328,000	7,500
Ref. 21	100,000	24,000

Table 6. Table 6: Total computation time of a complete I-V curve ranging from 1.8 V to 4.0 V (0.05 V step) for the Poisson-DD solver, Poisson-DD-Schrödinger solver, and Poisson-DD- 1 / u 1 𝑢 1/u solver, as tested by our home-built software. (*: This value has been obtained from the Table  5 assuming 16 iterations to achieve convergence.)

	Total computation time (s)
Poisson-DD	54,000
Poisson-DD-Schrödinger	*45,882,000
Poisson-DD-$1/u$	90,000

Equations

$$\{\begin{array}{l}\nabla \cdot \left(\epsilon \nabla \phi \right)=e\left(n-p+N_A^{-}-N_D^{+}\pm {\rho }_{pol}\right)\\ J_n=n{\mu }_n\nabla E_{Fn}\\ J_p=n{\mu }_p\nabla E_{Fp}\\ \nabla \cdot J_{n,p}=\pm e\left(A_0+B_{0}np+C_0\left(n^{2}p+n{p}^2\right)\right)\\ A_0=\frac{np-n_i^2}{{\tau }_n\left(p+n_i{e}^{\frac{\left(E_i-E_t\right)}{k_{B}T}}\right)+{\tau }_p\left(n+n_i{e}^{\frac{\left(E_t-E_i\right)}{k_{B}T}}\right)}\\ n={\int }_{E_c}^{+\infty }{\mathrm{DOS}}_{n,bulk}(E)\cdot f_n(E)dE\\ p={\int }_{-\infty }^{E_v}{\mathrm{DOS}}_{p,bulk}(E)\cdot f_p(E)dE\end{array}$$

$$\hat{H}=-\frac{{\hslash }^2}{2m_{e,h}^{\ast }}\Delta +E_{c,v}$$

$$\hat{H}{u}_{e,h}\left(\vec{r}\right)=1,$$

$$\begin{array}{r}n={\int }_{1/u_e}^{+\infty }{\mathrm{LDOS}}_{3D}(E)\cdot \frac{1}{1+\exp \left(\frac{E-E_{Fn}}{k_{B}T}\right)}dE\\ p={\int }_{-\infty }^{1/u_h}{\mathrm{LDOS}}_{3D}(E)\cdot \frac{1}{1+\exp \left(\frac{E_{Fp}-E}{k_{B}T}\right)}dE\end{array}$$

$${\mathrm{LDOS}}_{3D}(E)=\frac{\sqrt{2}{m}_{e,h}^{\ast }}{{\pi }^2{\hslash }^3}\sqrt{|E-1/u_{e,h}|}$$

$$x\left(r_i\right)=\frac{\sum _j\text{mbox}atom(j)\times e^{-\frac{{\left(r_j-r_i\right)}^2}{2{\sigma }^2}}}{\sum _j{e}^{-\frac{{\left(r_j-r_i\right)}^2}{2{\sigma }^2}}},$$

$$E_g^{\text{mbox}{In}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$$-1.4x(1-x),$$

$$E_g^{\text{mbox}{Al}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$$-0.8x(1-x),$$

$${\epsilon }_r^{\text{mbox}{In}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$${\epsilon }_r^{\text{mbox}{Al}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$$m^{\ast ,\text{mbox}{In}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$$m^{\ast ,\text{mbox}{Al}_x\text{mbox}{Ga}_{1-x}\text{mbox}N}=$$

$$\displaystyle{\bf P}^{pz}=\left[{\bf e}\right]\cdot\left[{\bf\epsilon}\right]=\left[\begin{array}[]{c}e_{15}\epsilon_{xz}\\ e_{15}\epsilon_{yz}\\ e_{31}\left(\epsilon_{xx}+\epsilon_{yy}\right)+e_{33}\epsilon_{zz}\end{array}\right]$$

$${\mathbf{P}}^{sp}=ax+b(1-x)+cx(1-x),$$

$$\nabla \cdot {\mathbf{P}}^{\mathrm{total}}={\rho }_{pol}\left(\vec{r}\right)$$

$$\mathrm{IF}=\frac{e}{k_{B}T}\frac{\partial V}{\partial \ln (J)}$$

$$n$$

$$p$$

$$J\propto N_c{N}_v\exp \left(-\frac{E_g}{k_{B}T}\right)\exp \left(\frac{eV}{k_{B}T}\right)\approx n_i^2\exp \left(\frac{eV}{k_{B}T}\right),$$

Full text

I Introduction

In a previous paper (referred to here as LL1, Ref. 1) we have introduced the new localization landscape (LL) theory to describe energy levels, localized states and density of states in disordered materials, with applications to the nitride materials case. A second paper  (refered to here as LL2, Ref. 2) showed the capability of the theory to compute materials system for time independent, frozen states entering optical absorption below the bandgap in alloy quantum wells (QWs). This phenomenon called the Urbach tail was shown to be well described by the density of states and the wavefunction overlap given by the landscape theory. The present paper will address the physical effects of disorder on carrier transport and recombination in a much more complex situation, that of the full multilayer heterostructures which make the light emitting diodes (LEDs), devices of utmost importance for energy savings and so far very poorly modeled due to the lack of a proper description of the effects of disorder.

Applying the LL theory to nitride-based alloys aims at solving long standing issues in these important materials. In recent years, nitride-based materials indeed play an increasing role in semiconductor markets including high power devices and light emitters due to the large bandgap range and full visible spectrum achieved by nitride ternary alloys.Weisbuch et al. (2015); Nakamura et al. (1991); Pengelly et al. (2012) However, a complete fundamental analysis of the intrinsic material properties is still missing to reconcile experiments and theoretical calculations. For instance, the explanations for the droop effect in nitride-based LEDs include electron overflow,Kim et al. (2007); Özgür et al. (2011) Auger recombination,Iveland et al. (2013); Shen et al. (2007); Binder et al. (2013a); David and Gardner (2010) poor hole injection,Ni et al. (2008) and carrier localization-delocalization.Mukai et al. (1999); Watson-Parris et al. (2011) While various direct measurements favor the Auger recombination, however requiring a large Auger coefficient when considering direct Auger recombination processes,Kioupakis et al. (2015) a number of other indirect experiments analyzed with simple or ad-hoc models form the basis for claiming the other mechanisms. All these above mechanisms are probably influenced by material compositional disorder effects which are neglected in device simulations or are only included at an elementary level. The disorder is best described from experiments including atom probe tomography (APT) with indium atoms distributed randomly in InGaN/GaN quantum wells (QW).Watson-Parris et al. (2011); Bennett et al. (2011); Shivaraman et al. (2013); Mazumder et al. (2013); Riley et al. (2014) Such compositional disorder due to local indium fluctuations should have an important role in determining electrical and optical properties in nitride-based LEDs or laser diodes. Recent atomistic calculations show that the localization effect induced by indium fluctuations will cause strong inhomogeneous broadening of the lowest transition energies.Schulz et al. (2015) Radiative recombination coefficients were also found to decrease with increasing indium concentration due to increased carrier localization in the random alloy fluctuations.Auf der Maur et al. (2016) However, atomistic simulations are impractical for calculating a full LED structure with contacts, multiple quantum wells (MQWs), electron blocking layer (EBL), and indium fluctuations. They do not really allow computing basic device characteristics, such as current-voltage curves and electro-luminescence spectra, where the determination of many quantum levels and the description of transport among such states is required.

In fully ordered materials, e.g., pure compounds, because of the small computation time and of the mature development of the technique, the classical drift-diffusion (DD) equations coupled with the Poisson equation are widely used to describe transport and optical properties although they treat carriers as semi-particles with a renormalized effective mass. Many results of quantum theory are however implicitly introduced as energy levels, density of states, quantum Fermi-Dirac statistics, and transport parameters, such as carrier mobilities and diffusion coefficients. In the case of LED simulations, radiative and non-radiative recombination mechanisms are described by a Shockley-Read-Hall coefficient $A$, a radiative coefficient $B$, and a non-radiative Auger recombination coefficient $C$. However, when comparing computational results using typical material parameters for ideal quantum wells or barriers without compositional disorder to experiments, the classical transport model leads to turn-on voltages either much too large in LEDs Yang et al. (2014); Karpov (2011); Chen et al. (2009) or too small in electron barriers.Browne et al. (2015) It also does not satisfactorily model the droop of the internal quantum efficiency (IQE).Piprek (2010) In general, researchers then use reduced piezo-electric polarization charge and larger Auger coefficient compared with experimental measurements to fit the experimental quantum efficiency Piprek (2010) and the turn-on voltage.Yang et al. (2014); Wu et al. (2012)

Quantum models such as the non-equilibrium Green’s function formalism (NEGF) have been proposed to model quantum transport properties  .Shedbalkar et al. (2016); Datta (1995) They provide a good description of quantum effects such as tunneling which are absent from DD models. However, to describe indium fluctuations Wu et al. (2015) multi-dimensional models are needed and the large burden of computation time ($>$ thousands of CPU hours) makes the NEGF approach impractical.

We previously incorporated indium fluctuations into the classical Poisson and DD model, by taking into account the fluctuating conduction and valence band potentials. We found that inclusion of these random potentials led to the enhancement of Auger recombination due to the higher local carrier density. In addition, the many percolation paths through local minima of potential energy for carrier diffusion did reduce the LED turn-on voltage.Yang et al. (2014); Wu et al. (2012) It was also computed that carrier localization induced by indium fluctuations has a strong influence on the broadening of the light emission spectrum.Yang et al. (2014) However, we still did not take into account the wave nature of electrons leading to localization and delocalization in a disordered potential.

In this paper, we will implement the LL theory Filoche and Mayboroda (2012); Arnold et al. (2016) into semiclassical Poisson and drift-diffusion equations. We model carrier dynamics including transport and recombination by using known parameters (mobilities, $A$, $B$, and $C$ coefficients) from experiments. The quantum effects affecting in-plane and perpendicular transport which arise from indium fluctuations are taken into account by effective electron and hole energies (the landscape energies), electron-hole overlap, disordered densities of states. Thanks to the efficiency of landscape computations, we carry out all calculations self-consistently for the disordered LED system.

II Simulation Methods

In this section, we describe the simulation framework, including how to apply the LL theory into Poisson and DD equations. In the standard classical picture, Poisson and drift-diffusion equations are solved self-consistently to obtain the conduction and valence band edges which are the potential energies $E_{c,v}$ for electrons and holes. The set of equations is

[TABLE]

where $\varphi$ is the electrostatic potential, $E_{Fn}$ and $E_{Fp}$ are the quasi Fermi levels of electrons and holes, DOS${}_{n,p,bulk}(E)=\sqrt{|E-E_{c,v}|}\cdot\sqrt{2}m^{*3/2}/\left(\pi^{2}\hbar^{3}\right)$ is the bulk density of states (DOS), $m^{*}$ is the effective mass of the electron or hole, and $f_{n}$ and $f_{p}$ are the Fermi distribution functions for electrons and holes, respectively. Note that if all the layers of the simulated semiconductor structure are considered as homogeneous, then the local free electron and hole carrier densities, $n$ and $p$, will only depend on the growth direction $z$. $E_{c}$ and $E_{v}$ are the local conduction and valence band potential energies, respectively. $N_{A}^{-}$ and $N_{D}^{+}$ are the activated doping densities of acceptors and donors, respectively. $\rho_{pol}$ is the density of polarization charges, which can be computed by taking the divergence of the total polarization in the space $(\nabla\cdot\bf{P}^{\rm total})$ including spontaneous and piezoelectric polarization fields. $J_{n}$ and $J_{p}$ are the electron and hole current densities, respectively. This paper uses the Shockley-Read-Hall (SRH) model to account for the defect-related non-radiative (NR) recombination through a rate $A_{0}$ [Eq. (1)], where $\tau_{n}$ and $\tau_{p}$ are the NR carrier lifetimes dependent on the growth condition, which are taken in this paper as $10^{−7}$ s (300 K), a value typically associated with low-defect density nitride LEDs.De Santi et al. (2016); Zhang et al. (2009) $k_{B}$ is the Boltzmann constant and $T$ is taken as room temperature. $E_{t}$ is the trapping energy level assumed to be located at the mid-gap, and $E_{i}$ and $n_{i}$ are the intrinsic energy level and intrinsic carrier density, respectively. $B_{0}$ ($3.0\times 10^{−11}$ cm3s*-1*) is the intrinsic radiative recombination coefficient. Our $B_{0}np$ rate represents the usual radiative recombination rate with $B$ including the effect of overlap of the wave functions across the active region.David and Grundmann (2010a) In contrast, $B_{0}$ is a bulk coefficient in our model and the separations of the electron and hole distributions in the QW both along z due to the quantum-confined Stark effect (QCSE) and in-plane due to indium fluctuations are considered in the $np$ term. $C_{0}$ ($2.0\times 10^{−31}$ cm6s*-1*, see Ref. 15) is the Auger recombination coefficient, where the influence of the electron and hole overlap is included in the $n^{2}p$ and $np^{2}$ terms. $B_{0}$ and $C_{0}$ are considered temperature-independent in our simulations.

At this point we should emphasize that the choices of $B_{0}$ and $C_{0}$ are somewhat arbitrary if even not inconsistent: as the computation takes into account disorder induced localization effects and QW confinement Stark effects, we should take values for the $B_{0}$ and $C_{0}$ bulk parameters without disorder and electric field, but any experimental value will incorporate such effects of disorder and QW confinement. Only those values for bulk binary compounds such as GaN would be disorder effect free, but then the effect of the bandgap change with QW alloying might be significant. For $B_{0}$, Kioupakis et al. Kioupakis et al. (2013) calculated that the change from GaN to In0.25Ga0.75N modifies $B_{0}$ by less than 10% without taking disorder into account, with $B_{0}$=$6.4\times 10^{-11}$ cm3s*-1* for 25% indium. On the experimental side, the bulk $B_{0}$ value ($7.0\times 10^{−11}$ cm3s*-1*) extracted from InGaN/GaN double heterostructure (DH) LED experiments would give a value without QCSE, although with disorder, but it assumes that $n$ (carrier density) is constant across the structure.David and Grundmann (2010b) However, in real DHs, the carrier density is not a constant as a large polarization electric field is present in the DH. A very high carrier density still exists locally to screen the polarization field in the DH which locally increases the radiative recombination rate compared to the space-averaged one. Therefore, the $B_{0}$ value obtained in a DH while assuming uniform carrier concentration is not correct and cannot be taken as the $B_{0}$ bulk input parameter without disorder and electric field as it is an overestimation. Thus, we rather take an experimental value of $3.0\times 10^{-11}$ cm3s*-1* for $B_{0}$ David and Grundmann (2010a) in QWs as a better approximation in this simulation, although a parametric evaluation should be done at a later point, both experimentally and theoretically.

The $C_{0}$ issue is even more complex: the Auger recombination term $Cn^{3}$ (more precisely $C_{eeh}n^{2}p+C_{ehh}np^{2}$, where $C_{eeh}$ and $C_{ehh}$ are the electron-electron-hole and electron-hole-hole Auger coefficients, respectively) incorporates disorder impact through two effects: increase in local carrier densities through carrier localization; increase in the Auger coefficient through wave function localization.Vaxenburg et al. (2013) The starting $C_{0}$ value for the computation should then be free of both effects, and cannot come from experiment. Theory Kioupakis et al. (2015) gives a value of $\sim 10^{-33}$ cm6s*-1* for the direct Auger process in InGaN and $\sim 10^{-31}$ cm6s*-1* for the indirect phonon-assisted Auger process. By taking alloy effects through a supercell model, Kioupakis et al. Kioupakis et al. (2015) calculate an alloy-assisted Auger coefficient of a few $10^{-31}$ cm6s*-1*, and an overall Auger coefficient of $\sim 10^{-30}$ cm6s*-1*. As we have not yet calculated the effect of alloy disorder on the Auger coefficient in the landscape model,Filoche et al. (2017); Lentali et al. (2017) we consider as a starting value for $C_{0}$ a low indium content experimental value from David and Grundmann David and Grundmann (2010a) of $\sim 2.0\times 10^{−31}$ cm6s*-1*. Doing this, we might underestimate the effect of the alloying on the Auger coefficient, but we will still capture the major effect of local carrier concentration increase due to alloying, as well as the effects of the internal electric field.

If we want to solve the Schrödinger equation to take the disorder-induced quantum effects into account, we should use the disordered potential energy in a Schrödinger solver to calculate the wave functions and eigen-energies. Then the carrier density distribution can be obtained by the wave function distribution and relative eigen-energy levels. When the carrier density is obtained, this should be plugged into the Poisson-DD solver and the corresponding equations should be solved iteratively until convergence. In addition, a self-consistent 3D Poisson-Schrödinger solver is highly time consuming, as will be discussed later in Appendix.

As a result, we rather apply the theoretical landscape model proposed by Filoche et al. Filoche and Mayboroda (2012); Arnold et al. (2016); Filoche et al. (2017) to obtain the equivalent semiclassical confining potential seen by the carriers, as described in LL1 (Ref. 1). According to this theory, an effective potential can be found which captures the complex interference pattern created by the carrier wave functions in the original disordered potential and transforms it into a semiclassical confining potential which localizes the carriers in different regions. Additionally, the long-range exponential decay, characteristic of Anderson localization,Anderson (1958) is explained as the consequence of multiple tunneling in the dense network of barriers created by this effective potential.Arnold et al. (2016) Therefore, both quantum localization/confinement and tunneling effects are described in the LL theory.

The Hamiltonian entering the Schrödinger equations for electrons and holes reads

[TABLE]

The landscapes $u_{e,h}(\vec{r})$ for electrons and holes are defined as the solutions of:

[TABLE]

and $1/u_{e,h}$ are the effective potentials incorporating the localization properties of the solutions of the Schrödinger equation.Filoche et al. (2017); Arnold et al. (2016) The boundary conditions for Eq. (3) can be either Dirichlet, Neumann or periodic. After obtaining $1/u_{e}$ and $1/u_{h}$ for electrons and holes, respectively, we use these as the input potential energies for the Poisson and DD equations to replace the original terms $E_{c}$ and $E_{v}$. In the carrier density calculation, $1/u_{e}$ and $1/u_{h}$ determine locally the bottom energy for the local DOS (LDOS) of the disordered system (Fig. 1):

[TABLE]

where $E_{Fn}$ and $E_{Fp}$ are the quasi Fermi energies for electrons and holes, respectively.

The LDOS in landscape theory can be computed from Weyl’s law in 3D:Arnold et al. (2016); Filoche et al. (2017)

[TABLE]

The fact that the DOS based on the LL is simply obtained by replacing the original potentials $E_{c,v}$ with the effective potentials $1/u_{e,h}$ makes the LL theory easily implementable into the classical Poisson and DD model. It should be emphasized at this point that Eq. (4) and (5) will give a good description of optical and transport properties, much better than through the use of a model with disorderless QWs, only because the LDOS spectrum is well described in the landscape theory as was shown in LL1 (Ref. 1), section IVB.2.

The schematic flowchart of the entire simulation process is shown in Fig. 2. First, the spontaneous polarization charges and piezoelectric fields are computed. After the Poisson equation is solved to obtain $E_{c}$ and $E_{v}$, the landscape equations for the conduction and valence band are solved giving the effective potentials $1/u_{e}$ and $1/u_{h}$ [Eq. (3)]. Then, the carrier densities of electrons and holes are calculated from Eq. (4) using $1/u_{e}$ and $1/u_{h}$, and fed back to the Poisson-DD equations to be solved in a self-consistent manner. When the potential energy difference between two consecutive iterations is smaller than $10^{-5}$ eV, we consider the simulation loop as having converged and the iterations stop.

In the following sections, we model realistic GaN-based LED structures with indium fluctuations to study the impact of disorder as calculated through the LL theory.

III The fluctuating potential in InGaN QWs and AlGaN EBL

We adopt In0.14Ga0.86N and Al0.15Ga0.85N as the average alloy composition of the QWs and EBL, whose dimensions are illustrated in Fig. 3(a). The indium and aluminum atoms are randomly distributed in the QWs and EBL and the composition maps shown in Fig. 3(b) and 3(c) are obtained via the Gaussian averaging method detailed below.

Note that other alloy systems, for instance InxGa1-xAs, are well described by models that do not capture any effect of disorder-induced localization, such as the virtual crystal approximation (VCA) in which each potentially disordered site is substituted by an artificial atom interpolating between the properties of the actual components. In the VCA the maps of Fig. 3(b) and 3(c) would be substituted by an homogeneous atom distribution with the mean composition of the alloy. However, nitride-based materials are characterized by composition fluctuations which induce large polarization-related local electric fields, large band offsets between GaN and InN, and heavy carrier effective masses, both larger than in the arsenide alloy system. Simulations based on the VCA in this context fail to provide a correct description of the local variations of the physical observables (density of states, carrier distribution, etc.).Nguyen et al. (2004)

Then one has to define the maximum length scale over which the rapidly oscillating distribution of atoms can be averaged to obtain a continuous fluctuating potential while preserving the effects of disorder on the electronic properties of the system. Such length scale can be predicted by the LL theory, which is able to provide the effective potential fluctuations “seen” by the carriers. However, to avoid using a circular argument, let us estimate this length scale from an independent general theory of disorder: Baranovskii et al. Baranovskii and Efros (1978) showed that the spatial scale of the fluctuations affecting the transport of the electrons and holes is given by the de Broglie wavelength $\lambda=\hbar/\sqrt{m^{*}E_{0}}$, where $m^{*}$ is the carrier effective mass and $E_{0}$ is the energy scale of the band edge broadening due to disorder. In LL2 (Ref. 2) a value of $E_{0}\approx 50$ meV or $\approx 30$ meV (fluctuations of $E_{c}$ and $E_{v}$ respectively) was calculated for InxGa1-xN QW layers over a large range of indium concentrations ($x$=10%-30%). Taking these values of $E_{0}$ and the carrier effective masses of GaN (Table 1, note that $m^{*}_{\rm GaN}>m^{*}_{\rm InGaN}$) gives a lower bound on the spatial size of fluctuations in InGaN alloys of 2.8 nm for electrons and 1.1 nm for heavy holes, of the same order as the fluctuations of the effective confining potentials obtained from the LL theory, as it will be shown in Fig. 5(b) and 5(d). In the computational framework presented here the atomic disorder is smoothed by a Gaussian averaging over a length scale of $2\sigma\approx 1.2$ nm, that is smaller than the spatial size of fluctuations seen by the carriers calculated from Ref. 42. Therefore the averaged atomic distribution we use still incorporates disorder on the relevant scale for carrier transport. The details of the algorithm generating the electric potential map from the atomic distribution are described in the following.

At first, a cubic grid is constructed with a spacing corresponding to the average distance between cation atoms in GaN ($a=2.833$ Å). Then we randomly assign at each cation site either an indium (aluminum) or gallium atom for the InGaN (AlGaN) alloy, as shown in Fig. 4. For each atom site $i$ the local averaged alloy composition $x(r_{i})$ is determined from the Gaussian averaging method as

[TABLE]

where the sum goes over all atom sites $j$ of the domain, $atom(j)$ is zero or unity as decided by the random generator, and $\sigma$ is the half width of the Gaussian broadening parameter.

In this paper and in LL2 (Ref. 2) we fixed $\sigma=2a$ ($\approx$ 0.6 nm), which gives an average alloy composition along the growth direction of an InGaN QW that matches APT data.Wu et al. (2012) In addition, we observe that when $\sigma=a$ is used, the indium composition map exhibits very large fluctuations, from 0% to more than 60%, that are strongly localized in space and behave like single-atom fluctuations, which are beyond the applicability of the effective-mass Schrödinger equation. Such choice of the Gaussian parameter impacts the results of the Poisson-drift-diffusion solver, which takes as an input the strongly fluctuating real potentials $E_{c}$ and $E_{v}$. However the calculated effective potentials used in the Poisson-drift-diffusion-landscape model are observed to be substantially unchanged for $\sigma$ values ranging from $a$ to $2a$, as the rapid fluctuations of the real potentials are smoothed by the LL theory.

In our computations, we separate the process of atom grid generation and computation mesh construction to make the random alloy generator independent from the mesh elements. At each mesh node all the material parameters used in the simulation (e.g., bandgap, dielectric constant, and effective mass) are assigned according to the local alloy composition map $x(r_{i})$. If the mesh node does not coincide with the atom grid position, the linear interpolation of grid map will be used to determine the composition. The III-nitrides material parameters we used are shown in Table 1. All parameters of InGaN and AlGaN alloys are obtained by an interpolation method as

[TABLE]

The band offsets between GaN/InGaN and GaN/AlGaN conduction bands are assumed to be 63% of the bandgap difference.

To model the 3D strain distribution in the disordered system, we adopt the 3D continuum strain-stress model solved by the finite element method (FEM) to calculate the strain distribution over the entire LED before solving the Poisson and DD equations.Hsu et al. (2015) The calculated strain is transformed into the piezoelectric polarization field as

[TABLE]

where $\epsilon_{xx}$, $\epsilon_{yy}$, $\epsilon_{zz}$ are normal strains and $\epsilon_{yz}$, $\epsilon_{zx}$, $\epsilon_{xy}$ are shear strains. ${\bf P}^{pz}$ is the strain-induced piezoelectric polarization. $e_{15}$, $e_{31}$, and $e_{33}$ are the piezoelectric coefficients (Table 2), and other terms are zero due to the symmetry of wurtzite crystal structures. On the other hand, the spontaneous polarization values related to the GaN buffer layer is obtained by the following equation:

[TABLE]

where the $a$, $b$, and $c$ coefficients can be found in Table 3. By taking the divergence of the total polarization ${\bf P}^{\rm total}$, which includes the spontaneous and piezoelectric polarization field (${\bf P}^{\rm total}={\bf P}^{sp}+{\bf P}^{pz}$) over the entire domain, we calculate $\rho_{pol}(\vec{r})$ as

[TABLE]

This induced fixed polarization charge $\rho_{pol}$ at different locations is finally implemented into the Poisson equation as the initial condition, as shown in the flowchart of Fig. 2.

IV The localization landscape in nitride LEDs

In this paper, the size of the simulated domain is 35 nm $\times$ 35 nm $\times$ 387 nm with a full LED structure including six-pairs MQW, an EBL, and $p$ and $n$ transport layers, as shown in Fig. 3(a). The simulation will model the current injection, transport and carrier screening of polarization fields, by solving the equations self-consistently. The geometric structure was meshed by the Gmsh program,Geuzaine and Remacle (2009) where the mesh has 1,265,291 nodes and 7,662,428 tetrahedral elements. The mesh grid size is 0.5 nm $\times$ 0.5 nm in the $x$-$y$ plane and a gradual mesh technique was used for the grid size in the $z$-direction ranging from 0.12 nm to 20 nm. The Schrödinger equation is of course an eigenvalue problem, whereas the LL model solves a much simpler linear equation. This significantly reduces the computation time in each iteration step by a factor of $\sim$ 1000 as compared to a Schrödinger solver. In addition, the computation time of the landscape equation is approximately the same as that for the DD equation (with both electrons and holes). The detailed computation time required to solve each equation and a comparison with other models are given in the Appendix.

We now implement the LL theory into the Poisson-DD model and solve these equations self-consistently (Fig. 2) to account for disorders. This solves the carrier density and transport including the quantum effects of disorder inasmuch the landscape model results represent those of a Schrödinger solver.Filoche et al. (2017) At a given applied bias to the LED structure the 3D LL is computed in a self-consistent manner starting from the original electron and hole potentials. As an illustration, the 2D energy potential maps corresponding to the mid-plane of a QW and the 1D band diagram of the structure along the $z$-direction are shown in Fig. 5 and Fig. 6, both for the original and the effective potentials $1/u_{e,h}$ at a bias of 2.8 V.

Figures 5(a) and 5(b) show the conduction band potential and the corresponding landscape potential $1/u_{e}$ computed self-consistently in presence of QW disorder. $1/u_{e}$ appears to be smoother compared to $E_{c}$ because the landscape theory flattens the rapid fluctuations not “seen” by the quantum states of the disordered system.Filoche et al. (2017); Arnold et al. (2016) The local peak potentials are lowered and smoothed due to quantum tunneling effects. Besides in-plane quantum effects, the energy reference of $1/u_{e}$ in Fig. 5(b) is also raised with respect to $E_{c}$ by quantum confinement along $z$ [see Fig. 6(b)]. Note that the difference between the $E_{v}$ and $1/u_{h}$ maps is smaller, as shown in Fig. 5(c) and 5(d). This is due to the heavy hole effective mass (1.829 $m_{0}$), which is much heavier than the electron effective mass (0.159 $m_{0}$) in the In0.14Ga0.86N QW. As a result, the quantum confinement and tunneling effects experienced by heavy holes are much reduced.

Similar effects can be clearly observed in the band diagram of the structure shown in Fig. 6 obtained when averaging in-plane $1/u_{e}$ and $1/u_{h}$: $1/u_{h}$ remains fairly close to the valence band edge, while in $1/u_{e}$ the minima of the conduction band edge are raised considerably and the barriers are appreciably lowered. We emphasize that the LL theory captures the quantum effects in a disordered layer over a wide range of energies, far from only the ground states, and thus describes the effective band diagram of a quantum semiconductor structure (see LL1, Ref. 1), therefore influencing both in-plane and perpendicular carrier transport through the LED. This will prove essential in the threshold voltage for carrier transport through the heterostructure, and is at the root of the proper description of the perpendicular transport, while using layers with homogeneous averaged composition lead to vastly overestimated threshold voltages.

V Electrical and optical properties calculated from the localization landscape theory

After understanding how the LL theory predicts the effective local potentials for electrons and holes used to calculate the carrier concentrations in a self-consistent loop with the Poisson-DD equations, we move to study current-voltage characteristics, carrier densities and the quantum efficiency of the whole GaN-based LED structure. The simulation parameters can be found in Table 4 including the doping and recombination coefficients. For the sake of comparison we compute different optoelectronic properties of the LED using the landscape-based Poisson-DD solver with random alloy fluctuations (“$1/u$-Poisson-DD”) and compare it to the classical Poisson-DD method with random alloy fluctuations (“Poisson-DD”) or with uniform QWs and EBL (“Poisson-DD (uniform)”), developed in our previous work.Wu et al. (2012) Note that only the first modeling approach takes into account quantum effects due to disorder because of the implementation of the LL theory.

In Fig. 7 we compare the carrier distribution of electrons and holes computed in the mid-plane of the third QW of the LED structure. The classical Poisson and DD solver treats carriers as particles and, as a consequence, the local carrier density fluctuates strongly reflecting the rapid spatial oscillations of the alloy composition [Fig. 7(a),(c)].

On the other hand, in the landscape model quantum effects are included and the smoother effective potentials [Fig. 5(b),(d)], via Eq. (4), produce more uniform carrier distributions which are better representative of the standing wave nature of the localized quantum states [Fig. 7(b),(d)].

The I-V characteristics of the LED calculated by different modeling approaches are shown in Fig. 8. Usually,Yang et al. (2014) the classical Poisson-DD model without disorder leads to a very large turn-on voltage when assuming 100% theoretical polarization charge, and in several works Piprek (2010); Wu et al. (2012) an internal charge reduced by 50% was used to realize a turn-on voltage and IQE more in line with experiment. Both calculations with disorder, “Poisson-DD” and “$1/u$-Poisson-DD”, use the same input random indium distribution and 100% theoretical polarization charge. Due to the higher effective bandgap of the $1/u$ potentials [Fig. 6(b)], the current density computed using the landscape is slightly lower before or near turn-on voltage compared to the classical Poisson-DD solver with adjusted 50% polarization charge. When the applied voltage increases above the threshold, the current density computed using the landscape becomes larger due to the lower effective barriers of the $1/u$ potentials [Fig. 6(b)], which effectively reduce the internal resistance for carrier transport. The forward voltage ($V_{f}$) computed by $1/u$ corresponding to 20 A.cm*-2* is around 3.0 V at 300 K, which matches quite well commercial blue LED data (2.8-3.1 V). For instance, Nichia Co. reported in Ref. 52 two blue LEDs, a high-efficiency one with $V_{f}$=2.89 V at 10 A.cm*-2* and a high-power one with $V_{f}$=3.10 V at 35 A.cm*-2*, both values being in good agreement with the I-V characteristic predicted by the $1/u$-Poisson-DD model. The remaining difference could be attributed to leakage paths via V-pit structures or dislocation lines,Li et al. (2016) to the absence in our modeling of tunneling in perpendicular transport, or to the internal temperature of real LED devices being higher than 300 K, leading to a lower $V_{f}$.

The average carrier densities along each QW of the active region are shown in Fig. 9. The black curves correspond to a structure with uniform In0.14Ga0.86N QWs and polarization charge reduced to 50% solved by the classical Poisson and DD model, whereas the blue curves are the QW calculated by the Poisson and drift-diffusion implementing the landscape.

Since at a given LED bias voltage the injected current in the two structures, with and without disorder, can be very different [see Fig. 8], it is more relevant to compare the carrier densities realized from the different models at a given current density, therefore in comparable conditions of band filling and carrier-induced electric field screening. We also plot in red in Fig. 9 the calculation of the carrier injection when the disorder is taken into account only through the changes in conduction and valence band levels, without the use of the landscape model to account for localization effects. The results lie somewhat in between those of the uniform material model and the landscape model. As we can see, at the same injected current density of 20 A.cm*-2*, the $1/u$-Poisson-DD model predicts a smaller carrier density as a consequence of the larger electron-hole spatial overlap and higher recombination rates induced by disorder. Moreover the modeling based on the LL shows that carriers are still quite inhomogeneously injected [see Fig. 16(b)], leaving room to improvement through better designs of the active region as one wishes homogeneous carrier injection to diminish the highly nonlinear Auger recombination.

Figure 10 shows the simulated effective $B$ coefficient, $B_{\rm eff}$, for the last QW (p-side) as a function of current density, defined as the recombination rate divided by the product of the QW averaged electron and hole concentrations. As can be seen, $B_{\rm eff}$ increases with current density, mainly due to electric field screening. The increase is particularly important for the simulated uniform QWs, as disorder will smooth out large potential fluctuations. The size of the calculated change in $B_{\rm eff}$ however points that the frequently used ABC model with constant $A$, $B$, $C$ parameters is of little use to quantitatively analyze recombination phenomena in LEDs.Weisbuch et al. (2015)

Figure 11 shows the $z$-component of current and the corresponding effective quantum potentials in the mid-plane ($x$-$y$ plane) of the third QW. As shown in Fig. 11(d), the current finds percolation paths through lower potential regions, enabling a decrease of the turn-on voltage, where a much smaller voltage is needed to reduce the polarization induced potential barrier. The existence of such percolation paths with high current densities is the basis of the reduced turn-on voltage caused by indium fluctuations induced by disorder, together with the in-plane averaged effects of confinement and tunneling (see Fig. 6). Similar effects of InGaN alloy fluctuations has been observed and verified in unipolar structures.Browne et al. (2015); Fireman et al. (2016)

VI Discussion

The above calculation of the I-V characteristic based on $1/u$ only assuming experimental material parameters (Tables I-IV) represents well the experimental data. We turn now to more precise results to identify the calculated LED internal features. Figure 12 displays the ideality factor (IF) of the I-V characteristics calculated as

[TABLE]

The $1/u$ simulation remarkably reproduces the few experimental data available in high quality LEDs.Auf der Maur et al. (2016); Binder et al. (2013b); Zhu et al. (2009) At low bias and current, IF is near 2, as in the Sah-Noyce-Shockley theory Sah et al. (1957) due to SRH recombination in the depletion region (here mainly in the QWs). Increasing the bias, the IF diminishes close to unity, as modeled in perfect p-n junctions where current is dominated by diffusion in the neutral regions of the junction. This is to be expected: in the bias region where recombination is dominated by a bimolecular radiative recombination process, the current density is approximately $B_{0}~{}np$. Expressing $n$ and $p$ as

[TABLE]

This yields

[TABLE]

which corresponds to IF=1. The minimum IF, near unity, is reached at a current density of the order of 0.1 A.cm*-2*, in the range of the experimental measurements.Auf der Maur et al. (2016); Binder et al. (2013b); Zhu et al. (2009); David et al. (2016) At even larger bias, the IF increases again as the series resistance dominates the device characteristic. Our IF calculation shows that the “$1/u$-Poisson-DD” model including random alloy fluctuations provides an excellent overall description of the transport properties of LEDs.

It must be remarked that an almost ideal IF does not mean at all that carriers are uniformly distributed in the structure (cf. Fig. 9) pointing to the inherent shortcomings of electrical measurements to assess that critical phenomenon. We can also calculate the leakage current exiting the active layer region, as shown in Fig. 13. As can be expected from the low turn-on voltages modeled, quite smaller than the GaN bandgap, very little carrier leakage is expected under usual operating conditions. The calculation indeed shows that leakage is negligible until the bias voltage reaches 3.4 V, value corresponding to the GaN bandgap, with total currents above 500 A.cm*-2*. In this calculation, the sheet resistance is not considered since only vertical transport is calculated. If sheet resistance would be considered, one would even need higher voltage to reach the flat band condition since an extra potential drop will occur in the current spreading layer.

Turning to the light emission efficiency curve, shown in Fig. 14, the peak IQE obtained from $1/u$ is slightly higher than the classical Poisson-DD model on fluctuated QWs because of the more homogeneous in-plane radiative efficiency and carrier injection in each QW. As Fig. 15 shows, the in-plane radiative recombination distribution spreads out in the $1/u$-Poisson-DD model because carriers are localized on larger domains in the QW plane [Fig. 9(b),(d)]. Let us recall that the influence of wave function overlap in our model is included in the $np$ term instead of the radiative recombination coefficient $B_{0}$. Therefore, as shown in Fig. 9, the distribution of electrons and holes along the growth direction is more symmetric with respect to the mid-plane of the QWs due to quantum confined wave properties. They have then better overlap compared to the classical DD model, leading to a higher radiative efficiency.

Figure 16(b) shows the calculation of the integrated recombined current both radiatively and nonradiatively in each well for the uniform QW and disordered QW models at a current density of 20 A.cm*-2*. The inhomogeneity of carrier injection is well displayed by the total recombination current decreasing from the p-side, on the right, to the n-side. Moreover, the IQE of each QW is also calculated to examine the contribution of each QW as shown in Fig. 16(a). The IQE of the QW close to the p-layer is smaller in the classical Poisson-DD model due to a smaller overlap as compared to the $1/u$-Poisson-DD model. The opposite trend between the classical Poisson-DD model (higher IQE in the QWs close to the n-layer) and $1/u$-Poisson-DD model (higher IQE in the QWs close to the p-layer) also reveals a better current injection through MQWs for the $1/u$-Poisson-DD model due to reduced effective barriers for electrons. Besides, the Auger recombination will start to dominate in the last QW (p-side well) at higher carrier densities. Due to the increasing local carrier density in the disorder case, the integrated Auger recombination can be enhanced compared to normal QW (see Fig. 13), which matches well the experimental droop with a $C_{0}$ coefficient in line with theory. On the other hand, the different IQE values among QWs can be attributed to the inhomogeneity of carrier injection, which again shows that adopting a constant carrier density ($n$) in the ABC model is incorrect to represent the IQE behavior in MQW structures. While small, the recombination currents in QWs other than the last one cannot be neglected, as is sometimes done in ABC modeling. In the 50% polarization uniform QWs case, the intentionally reduced polarization field makes electrons and holes overlap better due to the more homogeneous distribution in each well, which does not match experimental observations.David et al. (2008) It also reduces the influence of Auger recombination due to lower local carrier density (without localization), so that the droop effect is less pronounced. A larger $C_{0}$ would then be needed to represent the experimental Auger recombination.

VII Conclusions

In this work, we successfully implemented a novel method, namely the LL theory of disordered systems, to model the carrier transport and optical emission of LED heterostructures including the effects of intrinsic disorder in nitride-based material alloys. According to the LL theory, a function $1/u$ acts as an effective semiclassical confining potential which allows us to account for in-plane quantum confinement and tunneling effects due to the random indium fluctuations. The carrier dynamics can then be modeled through the classical drift-diffusion equations in an efficient self-consistent way. With the landscape model, computations are much faster than the conventional Schrödinger eigensolvers, especially in 3D, typically by a factor $\times$100-1000, allowing self-consistent calculations.

The I-V characterization of LEDs matches very well experimental measurements as a result of reduced energy barriers and of percolative transport. While the ideality factor of the LED has a near perfect behavior as a function of injected current, the carrier distributions are still very inhomogeneous throughout the heterostructures. This shows that measuring electrically the structure does not provide a clear insight on the internal electronic processes. On the optical properties side, the landscape maps for electrons and holes give us a good estimate of electron and hole overlap, leading to more accurate simulations of LED IQE.

In principle, this method is not only restricted to modeling nitride-based devices, but can be expanded to model other disordered semiconductor materials and structures.

Acknowledgements.

This work was supported by the Ministry of Science and Technology in Taiwan for Taiwan-ANR collaboration project under Grant MOST 104-2923-E-002-004-MY3, 105-2221-E-002-098-MY3 and by the French National Research Agency (ANR) under grant ANR-14-CE05-0048-01. Svitlana Mayboroda is partially supported by the Alfred P. Sloan Fellowship, the NSF CAREER Award DMS-1056004, the NSF MRSEC Seed Grant, and the NSF INSPIRE Grant. Additional support for Marco Piccardo, James S. Speck and Claude Weisbuch was provided by the DOE Solid State Lighting Program under Award # DE-EE0007096.

Appendix A Computation time of the landscape theory applied to LEDs compared to other models

Table V shows the approximate computation time required for each iteration step when solving different equations, such as the Poisson, drift-diffusion, LL and Schrödinger equation. The total computation time required to compute a typical LED I-V curve (45 bias values) using different methods is given in Table 6, where an average number of 16 iteration steps until convergence is needed for each bias (Poisson-DD and Poisson-DD-$1/u$). The ARPACK solver Lehoucq et al. (1997) for the generalized eigenvalue problem was used. For the solution of the inverse problem, the PARDISO solver Schenk et al. (2001) was used. Typically, the environment of clusters to which the computation is submitted is 2 Intel Xeon E5-5650V2 8 cores 2.6 GHz CPUs with 396GB memory.

Let us draw a comparison with the computation time deduced from the simulations of InGaN/GaN QWs incorporating random indium fluctuations as reported by other groups (Table 5). D. Watson-Parris et al. used the finite difference method to solve the 3-D effective mass Schrödinger equation,Watson-Parris et al. (2011) where the node size is about 1,500,000. Their computation time for one iteration is about 60,000 seconds, which is quite similar to our Schrödinger simulation model and is extremely time consuming. Besides, the self-consistent loop cannot be done due to this long computation time so that the Poisson-Schrödinger solver is not self-consistent. Concerning atomistic simulations, S. Schulz et al. adopted the empirical tight binding method (TBM) and valence force field model to account for the strain-induced polarization field and band structure of the QW, while the perpendicular carrier transport was overlooked.Schulz et al. (2015) The simulation domain is limited near the single QW region of 10 nm $\times$ 9 nm $\times$ 10 nm containing $\sim$ 82,000 atoms. Hence, the typical node size is about 82,000$\times$4, with an estimated computation time of 7,500 seconds. M. Auf der Maur et al. also applied the empirical TBM to model indium fluctuations in InGaN QWs.Auf der Maur et al. (2016) Although the classical Poisson and drift-diffusion model was used to solve the electrostatic potential, the atomistic calculation is not performed self-consistently with the classical model. The dimension of the atomistic simulation is still limited near the QW region (10 nm $\times$ 10 nm $\times$ 11 nm containing $\sim$ 100,000 atoms), which cannot be used to model full MQW LED structure in view of the computation time. Therefore, the key problem for quantum solvers such as the effective mass approximation Schrödinger equation or TBM, is the huge amount of computation time which makes such models impractical for the full structure LED simulation.

It can be concluded that the landscape model coupled to the Poisson-DD equations is much more computationally efficient with respect to state-of-the-art quantum solvers, while still incorporating quantum effects such as tunneling and quantum confinement.