Localization landscape theory of disorder in semiconductors I: Theory and modeling

TL;DR

This paper introduces a novel modeling approach using localization landscapes to predict quantum localization effects in disordered semiconductor alloys, enabling efficient carrier transport simulations without solving the Schrödinger equation.

Contribution

It develops a new computational model integrating localization landscape theory into drift-diffusion simulations for disordered semiconductors, capturing quantum effects efficiently.

Findings

Accurately predicts localization regions and energies of carriers.

Replicates Schrödinger equation results in one-dimensional structures.

Accounts for quantum tunneling and confinement effects without complex quantum calculations.

Abstract

We present here a model of carrier distribution and transport in semiconductor alloys accounting for quantum localization effects in disordered materials. This model is based on the recent development of a mathematical theory of quantum localization which introduces for each type of carrier a spatial function called \emph{localization landscape}. These landscapes allow us to predict the localization regions of electron and hole quantum states, their corresponding energies, and the local densities of states. We show how the various outputs of these landscapes can be directly implemented into a drift-diffusion model of carrier transport and into the calculation of absorption/emission transitions. This creates a new computational model which accounts for disorder localization effects while also capturing two major effects of quantum mechanics, namely the reduction of barrier height (tunneling effect), and the raising of energy ground states (quantum confinement effect), without having to solve the Schr\"odinger equation. Finally, this model is applied to several one-dimensional structures such as single quantum wells, ordered and disordered superlattices, or multi-quantum wells, where comparisons with exact Schr\"odinger calculations demonstrate the excellent accuracy of the approximation provided by the landscape theory.

Tables

Table 1. Table 1: Band structure parameters for wurtzite GaN and AlN alloys: band gap, relative permittivity, 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
AlN	6.0	10.31	0.32	0.30	2.68	0.26

Table 2. Table 2: Interpolation parameters for polarization in Al x Ga 1-x N.

a	b	c
-0.090	-0.034	0.021

Table 3. Table 3: Piezoelectric coefficients for wurtzite GaN and AlN alloys.

	$e_{33}$	$e_{31}$	$e_{15}$
units	(C.cm-2)	(C.cm-2)	(C.cm-2)
GaN	0.73	-0.49	-0.40
AlN	1.55	-0.58	-0.48

Equations

$$\hat{H}=-\frac{{\hslash }^2}{2m}\Delta +V.$$

$$\hat{H}u=-\frac{{\hslash }^2}{2m}\Delta u+Vu=1,$$

$$|\psi (\vec{r})|\le Eu(\vec{r}).$$

$$-\frac{{\hslash }^2}{2m}\left[\frac{1}{u^2}\mathrm{div}\left(u^2\nabla {\psi }_1\right)\right]+W{\psi }_1=E{\psi }_1.$$

$$⟨\psi ⟩\hat{H}\psi =\frac{{\hslash }^2}{2m}⟨u\vec{\nabla }\left(\frac{\psi }{u}\right)⟩u\vec{\nabla }\left(\frac{\psi }{u}\right)+⟨\psi ⟩\hat{W}\psi$$

$$|{\psi }_{{\vec{r}}_0}\left(\vec{r}\right)|\lesssim e^{-{\rho }_E(\overrightarrow{r_0},\vec{r})}.$$

$${\rho }_E\left(\overrightarrow{r_0},\vec{r}\right)=\underset{\gamma }{\min }\left({\int }_{\gamma }\sqrt{{\left(W\left(\vec{r}\right)-E\right)}_{+}}ds\right),$$

$$u$$

$$\mathrm{with}{\alpha }_i$$

$${\alpha }_i$$

$$=\frac{1}{E_i}⟨1⟩{\psi }_i$$

$$u\approx \frac{⟨1⟩{\psi }_0^{(m)}}{E_0^{(m)}}{\psi }_0^{(m)}$$

$${\psi }_0^{(m)}\approx \frac{u}{\Vert u\Vert }$$

$$E_0^{(m)}$$

$$=\frac{{\iiint }_{{\Omega }_m}u(\vec{r})d^{3}r}{{\iiint }_{{\Omega }_m}{u}^2(\vec{r})d^{3}r}$$

$$H\left(\vec{r},\vec{k}\right)=\frac{{\hslash }^2{k}^2}{2m_e^{\ast }}+E_c\left(\vec{r}\right),$$

$$\mathrm{IDOS}(\mathrm{E})$$

$$=\frac{2}{(2\pi {)}^3}{\iiint }_{\vec{r}}\left({\iiint }_{\frac{{\hslash }^2{k}^2}{2m_e^{\ast }}\le E-E_c(\vec{r})}{d}^{3}k\right)d^{3}r$$

$$\mathrm{IDOS}(\mathrm{E})={\iiint }_{\vec{\mathrm{r}}}\mathrm{LIDOS}(\mathrm{E},\vec{\mathrm{r}}){\mathrm{d}}^3\mathrm{r}$$

$$\mathrm{LIDOS}(\mathrm{E},\vec{\mathrm{r}})$$

$$=\frac{1}{3{\pi }^2}{\left(\frac{2m_e^{\ast }}{{\hslash }^2}\right)}^{\frac{3}{2}}{\left[E-E_c(\vec{r})\right]}^{\frac{3}{2}}$$

$$\mathrm{LDOS}(E,\vec{r})=\frac{1}{2{\pi }^2}{\left(\frac{2m_e^{\ast }}{{\hslash }^2}\right)}^{\frac{3}{2}}\sqrt{E-E_c(\vec{r})}$$

$$\mathrm{LDOS}(E,\vec{r})=\frac{1}{2{\pi }^2}{\left(\frac{2m_e^{\ast }}{{\hslash }^2}\right)}^{\frac{3}{2}}\sqrt{E-W(\vec{r})}$$

$$n(\vec{r})={\int }_{W(\vec{r})}^{+\infty }\frac{1}{1+e^{\frac{E-E_F}{k_{B}T}}}\mathrm{LDOS}(E,\vec{r})dE$$

$$\mathrm{IDOS}(\mathrm{E})$$

$$=\frac{2}{\pi }{\int }_z\sqrt{\frac{2m_e^{\ast }(E-W(z))}{{\hslash }^2}}dz$$

$$\mathrm{LDOS}(E,x)=\frac{1}{\pi }\sqrt{\frac{2m_e^{\ast }}{{\hslash }^2}}\frac{1}{\sqrt{E-W(z)}}$$

$${\mathrm{LDOS}}_{3D}(E,z)={\int }_{W(z)}^E{\mathrm{LDOS}}_{\text{mbox}2D}(E-E_1)$$

$$\times \left(\frac{1}{\pi }\sqrt{\frac{2m_e^{\ast }}{{\hslash }^2}}\frac{1}{\sqrt{E_1-W(z)}}\right)d{E}_1$$

$${\mathrm{LDOS}}_{\text{mbox}3D}(E,z)$$

$$=\frac{1}{2{\pi }^2}{\left(\frac{2m_e^{\ast }}{{\hslash }^2}\right)}^{\frac{3}{2}}\sqrt{E-W(z)}$$

$$-\frac{{\hslash }^2}{2m}\Delta u_K+(V(\vec{r})+K)u_K=1$$

$$|\psi (\vec{r})|\le (E+K)u_K(\vec{r})$$

$$\{\begin{array}{l}\mathrm{div}\left({\epsilon }_r\vec{\nabla }\phi \right)=\frac{q}{{\epsilon }_0}\left(n-p-N_D^{+}-N_A^{-}\pm {\rho }_{pol}\right)\\ \mathrm{div}({\vec{J}}_n)=R+G_n\\ \mathrm{div}({\vec{J}}_p)=-R+G_p\\ {\vec{J}}_n=n{\mu }_n\vec{\nabla }{E}_{Fn}\\ {\vec{J}}_p=p{\mu }_p\vec{\nabla }{E}_{Fp}\\ n(\vec{r})={\int }_{E_c}^{+\infty }\frac{1}{1+e^{\frac{E-E_{Fn}}{k_{B}T}}}{\mathrm{LDOS}}_n(E,\vec{r})dE\\ p(\vec{r})={\int }_{-\infty }^{E_v}\frac{1}{1+e^{\frac{E_{Fp}-E}{k_{B}T}}}{\mathrm{LDOS}}_p(E,\vec{r})dE\end{array}$$

Full text

I Introduction

Alloy semiconductors are ubiquitous in many modern semiconductor devices, where the use of heterostructured materials can drastically improve the device performances (see, e.g., Ref. 1). The main property engineered here is the bandgap difference between the various materials associated in the heterostructures. However, due to the large lattice mismatch between pure compound semiconductors which would lead to highly defective materials if employed as such, one has to associate binary compounds and alloys of binary compounds, leading to ternary alloys, eventually quaternary alloys. The resulting lattice constants are sufficiently close to obtain growth of high-quality materials. One thus retains part of the band-gap discontinuities between the binary compounds which allows one to confine carriers in double heterostructures or quantum wells.Alferov (2001); Kroemer (2001) An additional useful property is the modulation of the refractive index, which proved to provide an additional beneficial effect of crucial importance to achieve room-temperature continuous wave lasers, opening the way to optical telecommunication systems. Many other properties are impacted by alloying, some of them not desirable: for instance, the compositional fluctuations induce an additional scattering mechanism for charge carriers which diminishes the carrier mobilities.

GaN-based compounds are among the semiconductors displaying the largest variety of properties of heterostructures and alloys. In addition to the usual features, they show, due to their large ionic composition and crystalline structure (most often wurtzite along the $c$ axis), spontaneous and piezoelectric fields which strongly impact their electrical and optical properties.Bernardini et al. (1997) They also exhibit significant effects of the intrinsic spatial compositional fluctuations of their alloys. The random indium content in InGaN multiple quantum well (MQW) structures can vary locally from 10% to 23% within a few nanometers, for an average composition value of 17%.Wu et al. (2012) These fluctuations can induce a strong modification of carrier spatial distributions, of recombination rates, and of the overall light emission efficiency of the device. They have long been identified as responsible for the short carrier diffusion lengths which in turn lead to high emission efficiencies in spite of the high density of defects still present in the best grown heterostructures.Nakamura (1998)

Accounting for carrier localization induced by the local material disorder in semiconductor materials and devices is a daunting task. Usually, it requires solving the Schrödinger equation (both for electrons and holes) for random realizations of the disordered potential, and determining the energies and the spatial structure of the localized quantum states. Actually, this problem relates to the famous Anderson localization phenomenon Anderson (1958) which has triggered an enormous literature, and several of whose aspects still remain puzzling after more than 50 years of research.Lee and Ramakrishnan (1985); Evers and Mirlin (2008); Lagendijk et al. (2009) From 1983, early evidence of the role played by localization in semiconductors has been observed in conductivity measurements around the metal-insulator transition of three-dimensional (3D) doped charge-uncompensated silicon.Rosenbaum et al. (1983) Theoretical approaches such as the self-consistent scaling theory of Anderson localization have been able to successfully describe this disorder-induced phase transition.Vollhardt and Wölfle (1982) More recently, Anderson localization has been found to play a key role in the transport properties of low dimensional media such as disordered graphene.Ostrovsky et al. (2006) Yet, most of the theories developed to account for Anderson localization rely on statistical averages (through correlation functions) and scaling hypotheses.Mello et al. (1988); Hofstetter and Schreiber (1993) Moreover, the fermionic nature of the carriers and the electron-electron and electron-hole interactions (through the electric field) increase the difficulty to reproduce the complex and intricate behaviors of semiconductors. The existing methods for computing transport needs to go down to the atomistic level using techniques such as non-equilibrium Green’s functions and coherent potential approximation, at a very high computational cost.Kalitsov et al. (2012) As a consequence, one still lacks a model able to predict efficiently carrier localization and its consequences in realistic disordered semiconductor devices, with the often added complication of multi-layered heterostructures.

One of the most puzzling aspects of Anderson localization is the strong spatial confinement of the one-particle quantum states, attributed to destructive interferences between different propagation pathways in a disordered potential. Recently, a new theory has been proposed, which allows one to accurately predict the localization regions of the carriers, and the density of states (DOS) in the disordered potential created by the fluctuations of material composition, without having to solve the Schrödinger equation.Filoche and Mayboroda (2012); Arnold et al. (2016) This groundbreaking theory is based on a new mathematical tool, the localization landscape (LL), which is the solution to a Schrödinger-type equation with uniform right-hand side.

We present here the implementation of this tool into semiconductor materials, and its use in a semi-classical transport model of semiconductor devices. We show how the formalism of the theory enables us to efficiently predict in semiconductor structures the wave functions and eigen-energies of the confined states, the overlap between electrons and holes, the DOS, and the carrier distribution. This implementation conserves a local formulation, adding only to Poisson and transport equations a different partial differential equation (PDE). It accelerates the computation time by several orders of magnitude compared to the Schrödinger-Poisson-drift-diffusion (DD) type approach.

The results of this localization theory are further applied to the specific case of nitride semiconductors in two companion papers, one showing experiments and theory of the Urbach tail of InGaN quantum wells (Piccardo et al.,Piccardo et al. (2017) hereafter called LL2), the other on the simulation of full light-emitting diode (LED) structures (Li et al.,Li et al. (2017) hereafter called LL3).

II The localization landscape theory

We present first the main features of the LL theory introduced in Refs. 17; 18. According to this theory, the precise spatial location of quantum states in a potential $V(\vec{r})$ can be predicted using the solution $u(\vec{r})$ of a simple associated Dirichlet problem, called the localization landscape. The quantum states and the energies of particles with mass $m$ are, respectively, the eigenfunctions and the eigenvalues of the Hamiltonian of the system defined as

[TABLE]

With this notation, the landscape $u(\vec{r})$ is defined as the solution to

[TABLE]

It is shown in Ref. 17 that the subregions hosting the localized eigenfunctions are delimited by the valley lines of the graph of $u$, see Figs. 1(a) and (b) for the case of two-dimensional (2D) potential (for a 3D potential, these valleys would be surfaces). This property directly derives from a fundamental inequality satisfied by any eigenfunction $\psi$ of $\hat{H}$ with eigenvalue $E$, normalized so that its maximum amplitude is equal to 1 (see Ref. 17 for the proof):

[TABLE]

In other words, the small values of $u(\vec{r})$ along its valley lines Filoche and Mayboroda (2012) constrain the amplitude of $\psi$ to be small along the same lines and, as a consequence, localize low energy eigenfunctions inside the regions enclosed by these lines. The landscape $u$ therefore exhibits a partition of the entire domain into a set of subregions, each of these subregions localizing the carriers. But, as exposed hereafter, much more information can be extracted from the localization landscape $u(\vec{r})$.

II.1 Effective localization potential

Not only $u(\vec{r})$ controls the eigenfunctions, but also the function $W(\vec{r})\equiv 1/u(\vec{r})$ (homogeneous to an energy) can be interpreted as a confining potential that is related, among others, to the exponential decay of the Anderson localized states away from their main localization subregion. This property can be proved by transforming the original Schrödinger equation through the introduction of an auxiliary function $\psi_{1}$ such that $\psi\equiv u\psi_{1}$. A straightforward computation yields:

[TABLE]

One can see from this equation that the auxiliary function $\psi_{1}=\psi/u$ thus obeys a Schrödinger-type equation in which the original potential $V(\vec{r})$ has disappeared. Instead, a new function $W(\vec{r})$ now plays the role of an “effective confining potential”. The valleys of $u$, which are the boundaries of the localization subregions,Filoche and Mayboroda (2012) now correspond to the crest lines (or watershed lines) of $W$ and thus act as barriers to the auxiliary function $\psi_{1}$. This function $\psi_{1}$, as well as the initial eigenfunction $\psi$, is now localized in the basins of $W$ [see Fig. 1(c)].

Actually, it was proved that $W$ plays exactly the role of an effective potential thanks to the following identity satisfied by any quantum state $\ket{\psi}$:

[TABLE]

This identity shows that the energy $E$ of any quantum state $\ket{\psi}$ can never be smaller than the one it would have in a potential $W\left(\vec{r}\right)$. Consequently, according to Agmon’s inequality,Agmon (1982, 1985) the quantity $(W-E)$ controls exponentially the decay of $\psi(\vec{r})$ in the regions where $E<W$. Mostly, the eigenfunction decays exponentially with a rate proportional to $\sqrt{W-E}$ in the barriers where $W>E$. More precisely, the decay at point $\vec{r}$ of the amplitude $\psi_{\vec{r}_{0}}\left(\vec{r}\right)$ of an eigenfunction centered in $\vec{r}_{0}$ and of energy $E$ is expressed through the inequality:

[TABLE]

where $\rho_{E}(\vec{r_{0}},\vec{r})$ is the Agmon distance between $\vec{r}_{0}$ and $\vec{r}$. This Agmon distance (depending on $E$) is defined as:

[TABLE]

where $\gamma$ minimum is the geodesic path (i.e. the path of minimum length) between points $\vec{r}_{0}$ and $\vec{r}$ and $ds$ is the elementary path on that geodesic.

The LL therefore provides an estimate of the decay of the quantum state away from its main existence region. This decay corresponds to the tunneling effect and is more commonly known in quantum mechanics as a result of the Wentzel-Kramers-Brillouin (WKB) approximation. The theory used here is its mathematical generalization, which holds for any potential $W$ satisfying $\braket{\psi}{\hat{H}}{\psi}\geq\braket{\psi}{\hat{W}}{\psi}$ for all quantum states. The fact that $W$ plays the role of an effective potential finely shaping the quantum states is crucial for deriving an accurate expression of the density of states, as we will see in Sect. II.3. Finally, our estimate of the eigenfunction amplitudes based on $W$ can also be used, for instance, to assess the coupling between distant localization subregions.Bondar and Liu (2011)

II.2 Eigenvalue and eigenfunction estimates

In each of the subregions bounded by the valley lines of the LL $u(\vec{r})$, the local fundamental eigenfunction and its corresponding energy can also be accurately determined from $u$ itself inside the localization region (the decay far from the localization region being assessed as exposed in Section II.1). To this end, the landscape $u$, satisfying $\hat{H}u~{}=~{}1$, has to be decomposed on the basis formed by the eigenfunctions $\psi_{i}$ of the Hamiltonian:

[TABLE]

The decomposition coefficients $\alpha_{i}$ can be computed using the self-adjointness of the Hamiltonian:

[TABLE]

From Eqs. (8) and (II.2), one can draw three main remarks. First, the lower energy quantum states contribute more to the landscape $u$ than the high-energy ones (because $E_{i}$ grows in the denominator of Eq. (II.2)). Secondly, in a given localization subregion, the low-energy states ($\psi_{i}$) entering the decomposition of Eq. (8) are essentially the local quantum states of this subregion. Thirdly, in each subregion, the fundamental state has a bump-like shape, while the higher-energy ones, by orthogonality, take positive and negative values which cancel out so that the scalar products $\braket{1}{\psi_{i}}$ of Eq. (II.2) almost vanish. Note that this cancellation also occurs for the high-energy delocalized states of the system. As a consequence, in each localization subregion $\Omega_{m}$, the following relation is deduced:

[TABLE]

$\psi_{0}^{(m)}$ being the local fundamental state of subregion $\Omega_{m}$. This shows that the local fundamental state $\psi_{0}^{(m)}$ is almost proportional to $u$ in $\Omega_{m}$:

[TABLE]

Inserting Eq. (12) into identity (5) and using the fact that $W\equiv u^{-1}$ allows us to evaluate the fundamental energy $E_{0}^{(m)}$ from the landscape only:

[TABLE]

The LL $u(\vec{r})$ therefore provides a direct estimate of the fundamental energy in each of the localization subregions.

II.3 Density of states

Finally, the prediction of the localized energies extends to the prediction of the integrated density of states (IDOS), hence to its derivative, the density of states (DOS). Thanks to the LL theory, these quantities can be computed not only globally for the whole system under consideration, but also locally. We detail here the general case of a 3D system as well as the specific case of a 1+2D-system exhibiting two-dimensional translation invariance.

II.3.1 3D DOS

Due to the uncertainty principle, $\Delta x\Delta k\approx 2\pi$, each three-dimensional one-particle quantum state spreads in phase space $(\vec{r},\vec{k})$ on a volume of order $(2\pi)^{3}$. As a consequence, the number of energy states below a given energy $E$ (i.e., the counting function, also called integrated density of states IDOS) is asymptotically equivalent to $\mathcal{V}(E)/(2\pi)^{3}$ when $E\rightarrow+\infty$, where $\mathcal{V}(E)$ is the volume in phase space determined by $H(\vec{r},\vec{k})\leq E$. This asymptotic behavior is the so-called Weyl’s law.Strauss (2008) In a three-dimensional semiconductor, the Hamiltonian of an electron in the conduction band reads as

[TABLE]

where $E_{c}$ is the conduction band energy and $m_{e}^{*}$ is the effective electron mass. The IDOS deduced from Weyl’s formula is therefore:

[TABLE]

The factor 2 appearing in the numerator accounts for the spin degeneracy, and the integral within parentheses on the second line is simply the volume of a sphere in $k$-space. The IDOS can therefore be written as a space integral of a local quantity that is classically assimilated to a local IDOS, noted LIDOS$(E,\vec{r})$, in the local band structure approximation:

[TABLE]

with

[TABLE]

Differentiating this LIDOS with respect to $E$ gives the local density of states ${\rm LDOS}(E,\vec{r})$:

[TABLE]

One recovers here the classical expression of the local density of states for conduction electrons in a semiconductor.

In a disordered system, Eq. (5) shows that $W(\vec{r})$ acts as an effective potential on the particle. Subsequently, $W$ controls the distribution of energies for localized states in each localization subregion. One of the main consequences of this control, as shown on several examples in Ref. 18, is that $W$ can be used to accurately estimate the integrated density of states, hence of the density of states also, although the approximate DOS obtained after differentiating the IDOS is in essence less accurate. Practically, this is achieved by replacing the original potential (here $E_{c}$) by $W$ in the LDOS:

[TABLE]

This expression physically means that at each point $\vec{r}$, the local density of states is equal to the one of an infinite medium with identical material composition and a parabolic band whose minimal energy would be $W(\vec{r})$ (see Fig. 2). In the classical view of the local band approximation, $\vec{r}$ represents in fact a small volume compared to the typical size of the system, but large enough so that it can contain the local electronic states. This approximation which is routinely used when the wave functions are delocalized Bloch waves, is even more justified when dealing with localized eigenfunctions.

The reason for the quality of the approximate IDOS computed using $W_{c}$ (resp. $W_{v}$) instead for $E_{c}$ (resp. $E_{v}$) has been detailed in Ref. 18, but one can give here a short explanation. Weyl’s law is fundamentally based on an analogy between quantum and classical filling of phase space, the classical being the volume bounded by $H(\vec{r},\vec{k})=E$, and the quantum originating from the uncertainty principle which states that the volume of a quantum state in phase space is approximately constant (due to $\Delta x.\Delta k\approx 1)$. Thus, counting the number of quantum states of energy smaller than $E$ (the IDOS) comes back to assessing the size of the corresponding volume in phase space. It is known that in a disordered medium, localized states contribute to perturb the bulk IDOS and make it depart from the Weyl’s equivalent. Now, the effective potential, which can be seen a smoothed version of the original potential ($E_{c}$ in the case of electrons), exhibits much better defined wells and barriers. The spatial shapes of the quantum states are closer to the classical trajectories in this effective potential than in the original random or disordered potential. In others words, the effective potential is closer to what is “experienced” by the quantum wave if it were a classical particle. As a consequence, the validity of the classical-quantum analogy is strengthened, and the local density deduced from Weyl’s law using $W$ much more accurate than any other available estimate.

One needs here to express one word of caution. While IDOS$(E)$ can be understood as the actual number of states below energy $E$ over the spatial region of integration, the “local” density of states ${\rm LDOS}(E,\vec{r})$ cannot be understood as the actual number of states at point $\vec{r}$ as soon as the spatial fluctuations of this LDOS occur on a scale smaller than the typical spatial extension of the electronic states. It should rather be considered as an “effective” value, a useful tool for assessing physical quantities.

This remark is of particular importance when this LDOS is integrated over all possible energies to compute the local carrier density $n(\vec{r})$:

[TABLE]

This carrier density is apparently a local quantity at point $\vec{r}$, as in the classical DD model of transport,van Roosbroeck (1950) but thanks to the effective potential $W$ appearing in the LDOS, it now also encompasses the quantum confinement induced by the material disorder. As a consequence, this expression is much more accurate than what could be estimated from the classical similar expression with $E_{c}(\vec{r})$ instead of $W(\vec{r})$. At the same time, it can be easily implemented into a PDE model such as a drift-diffusion solver with Poisson-DD-continuity equations, instead of solving the eigenfunctions of the Hamiltonian.

II.3.2 1D DOS

The procedure described above for 3D states can be applied similarly to 1D systems. In this case, the IDOS reads:

[TABLE]

and the local density of states is:

[TABLE]

If one considers a 3D system with translational invariance in the two other directions $x$ and $y$, such as a quantum well, then the quantum states are products of 1D and 2D states in the $z$ direction and the $(x,y)$ plane, respectively. The 3D density of states can thus be deduced from the above 1D density of states along $z$ by convoluting it with the 2D LDOS:

[TABLE]

LDOS${}_{\mbox{2D}}$ being the 2D density of states for free particles, which is constant and equal to $m_{e}^{*}/(2\pi\hbar^{2})$ (we do not count the spin degeneracy here as it is already included in the 1D LDOS). This finally gives the 3D density of states:

[TABLE]

One can notice that the above expression is exactly identical to the 3D local density of states obtained in Eq. (19), except that here $W$ depends only on $z$.

II.4 Setting the potential reference

One needs here to underline a peculiarity of the LL. When solving the Schrödinger equation, the potential $V(\vec{r})$ experienced by the quantum particle can be defined up to a constant value $K$. If one shifts the potential by $K$, then the resulting energies are also shifted by the same constant $K$. However, this invariance does not hold for the landscape $u$. If one considers $u$ being the solution to Eq. (2), then the solution $u_{K}$ corresponding to the same potential shifted by a constant $K$ satisfies:

[TABLE]

If the constant $K$ is much larger than the typical energies of the quantum states, then $Ku_{K}\approx 1$. Therefore, the corresponding effective potential $W_{K}=1/u_{K}$ is very close to $K$. Inserting the effect of the potential shift on the energies into Eq. (3), the amplitudes of the quantum states are controlled by the landscape through the following inequality:

[TABLE]

In the situation where $Ku_{K}$ becomes very close to 1, this inequality is almost trivially satisfied. In other words, the constraint on $\psi$ exerted by the LL $u_{K}$ becomes weaker. This means that the constant $K$ has to be chosen in order to be as small as possible, in such a way that the Hamiltonian remains a positive operator (a condition of applicability of the LL theory). As a consequence, in semiconductor structures where one encounters large, smooth variations of potential superimposed on small scale random potentials, one should resort to solving the LL piecewise in regions over which the variation of the potential $V$ is of the same order of magnitude than the energies of the quantum states.

III The transport model

III.1 The new self-consistent approach

In this section, we present the first implementation of the LL theory in the physical description of processes in semiconductor heterostructure devices. The exposed model uses an hybrid approach where energy levels, wave functions and DOS are computed using the landscape theory, while a standard description is retained for carrier transport and statistics, as described by the DD equations of semiconductor textbooks.Jüngel (2009) Classically, the DD model is described through the following set of coupled equations whose unknowns are the electrostatic potential $\varphi$ and the quasi-Fermi levels $E_{Fn}$ and $E_{Fp}$:

[TABLE]

where $\varepsilon_{r}$ is the medium relative permittivity, $n$ and $p$ are the electron and hole densities, $N_{A}^{-}$ and $N_{D}^{+}$ are the activated doping densities of acceptors and donors; $\rho_{pol}$ is the polarization charge which appears in some semiconductors such as nitrides where electric polarization effects are important; $\vec{J}_{n}$ and $\vec{J}_{p}$ are the electron and hole currents, respectively, and $\mu_{n}$ and $\mu_{p}$ their mobility; $R$ is the recombination rate which includes all types of recombination processes (SRH, radiative, Auger), and $G_{n,p}$ are the carrier generation rates. Finally, LDOSn,p are the bulk local densities of states for electrons and holes, respectively.

The first equation of (27) is the Poisson equation which determines the electrostatic potential from the charge carrier densities. Second and third equations are the continuity equations for both carrier transports. The fourth and fifth equations of (27) are the semi-classical expressions of the current densities, derived from the Boltzmann equation assuming a linear collision kernel. These expressions correspond to the linear response theory in statistical physics. The mobilities $\mu_{n}$ and $\mu_{p}$ are effective parameters that summarize the scattering events and the quantum transport in the bulk materials, either pure compounds or alloys. Finally, the sixth and seventh equations of (27) compute the carrier densities from the quasi-Fermi levels $E_{Fn}$ and $E_{Fp}$ through the densities of states LDOSn,p.

In the usual DD model, Poisson and transport equations are solved self-consistently to obtain the converged electrostatic potential and the quasi-Fermi levels $E_{Fn}$ and $E_{Fp}$. More advanced models, called hydrodynamic, include also energy transport and hot carriers,Jüngel (2009) but at the expense of much larger computational time. In order to account for quantum confinement effects, one has to solve the Schrödinger equation using $E_{c}(\vec{r})$ and $E_{v}(\vec{r})$ as potentials for electrons and holes, respectively. The carrier densities are then deduced from the electron and hole wave functions, $\{\psi_{e,i}\}$ and $\{\psi_{h,j}\}$, and their corresponding energies, $\{E_{e,i}\}$ and $\{E_{h,j}\}$ through summation with Fermi-Dirac distribution:

[TABLE]

In the absence of currents, these carrier densities enter Poisson equation which in turn modifies the electrostatic potential, and then $E_{c}$ and $E_{v}$ (the band offsets $\delta E_{c,v}$ mentioned in Fig. 3 being local properties of the material). This is the usual (without disorder) self-consistent Poisson-Schrödinger scheme, depicted by the blue cycle in Fig. 3. Accounting for quantum transport adds another level of complexity, which consists in solving the Schrödinger equation to determine the equilibrium quantum states, and then compute transitions between these states.Watson-Parris et al. (2011) The dynamic carrier densities now enter Poisson and transport equations, from which one deduces the electrostatic potential and the quasi-Fermi levels. This loop results in lengthy, time-consuming and possibly unstable simulations that can take days of computation for a complex 3D structure, even damagingly longer with compositional disorder. At an even more fundamental level, non-equilibrium Green’s function techniques or atomistic tight-binding models are used to assess fast processes.der Maur (2015); Schulz et al. (2015)

Overall, DD-based models still remain to this date the only models able to simulate large structures and compute carrier transport in realistic semiconductor devices with a computational time compatible with optimization and design (see LL3, Ref. 20).

In the following, the implementation of the LL theory is first presented for the Poisson scheme, then for the full DD model including carrier transport. Introducing the LL $u(\vec{r})$ allows us to entirely bypass solving the Schrödinger equation, a computationally highly demanding step. From the electrostatic potential, we deduce two LLs $u_{e}$ and $u_{h}$, hence, two effective potentials $W_{c}$ and $W_{v}$. According to the theory, these landscapes provide direct estimates of the density of states LDOS${}_{n}(\vec{r})$ and LDOS${}_{p}(\vec{r})$ [see Eq. (22)], and of the carrier densities. This bypass is depicted by the red arrow in Fig. 3. The recombination rates are also assessed by computing the overlap between electron and hole states in the localization subregions (see next subsection). These quantities are then used to solve Poisson and DD equations. This makes the LL theory easily compatible with a classical DD approach. The schematic flow chart of the entire self-consistent simulation process is displayed in Fig. 4.

One has to underline that in this model, the current densities still take their classical DD form expressed in Fig. 4. In this form, the tunneling current caused by the spatial tails of the eigenstates that cross potential barriers is inherently accounted for in the model through the lowering of the effective potential $W_{c}$ (resp. $W_{v}$) as compared to $E_{c}$ (resp. $E_{v}$). However, the quantum current originating from phonon-assisted hopping between different eigenstates or by scattering events, is still accounted for through effective transport parameters, i.e., the mobilities $\mu_{n}$ and $\mu_{p}$. Extension of the theory to account for a full quantum model of transport is currently work in progress.

III.2 Computing generation/recombination processes

The final subsection of the modeling section is dedicated to the computation of the generation and recombination processes in direct gap semiconductors using the LL, the focus being in particular on optical interband transitions, namely absorption and emission. Typically, the computation of these processes in quantum-confined semiconductor structures requires the knowledge of the eigen-energies and eigenfunctions of the quantum states, which determine respectively the energy of the transitions and the overlap between the electrons and holes. As derived in Sec. II.2, these quantities can be in fact directly assessed from the LL. It should be noted that due to the integrals appearing in Eq. (12) and (II.2), these relations cannot be implemented into a local formulation, such as a self-consistent Poisson-DD-landscape loop, but can be evaluated after convergence is reached for the simulated system.

We consider in the following the case of a homogeneous quantum well (QW), for which the classical absorption/emission expressions are rewritten using the LL theory. The extension to the case of QWs with in-plane disorder is briefly outlined, while for more details the reader is referred to the companion papers LL2 and LL3.

III.2.1 Absorption

The quantum mechanical expression of the absorption coefficient $\alpha$ of a homogeneous QW of thickness $L$ as a function of photon energy $h\nu$ is Singh (2007)

[TABLE]

where the prefactor $C=\pi e^{2}\hbar|p_{cv}|^{2}/m_{0}^{2}cn_{r}\varepsilon_{0}h\nu$ depends on the real part of surrounding refractive index $n_{r}$ and the interband momentum matrix element $p_{cv}$. Considering a narrow photon energy range (close to the band gap), we neglect the $k$-dependence of $p_{cv}$ and the $h\nu$ dependence of $C$, and treat this prefactor as a constant; $g_{2D}=m_{r}/\pi\hbar^{2}$ is the joint density of states (JDOS) of a 2-D system including the spin degeneracy with $m_{r}=(1/m_{e}^{*}+1/m_{h}^{*})^{-1}$ being the reduced effective mass. The summation is performed over all electron and hole states (labeled $i$ and $j$, with $i=j=0$ being the fundamental states) with eigen-energies $E_{e,i}$ and $E_{h,j}$. The overlap factor $I_{i,j}$ is determined by the electron and hole eigenstates $\psi_{e,i}$ and $\psi_{h,j}$ as

[TABLE]

Our approach here is to use the LL theory to assess the absorption coefficient directly from the maps of $u_{e}$ and $u_{h}$ computed for electrons and holes using Eq. (2), where $V$ corresponds to the conduction and valence band potential, respectively. The valleys of the landscapes (which are also the crest lines of the $W_{c,v}$ potentials) partition the domain into localization subregions. Note that in the case of a homogeneous QW, only one localization subregion exists for each type of carrier, determined by the confinement along the growth direction, while in the case of a disordered QW several subregions can be found.Piccardo et al. (2017); Li et al. (2017) The determination of the valley lines in a homogeneous QW is discussed in Sec. IV.1.

In each subregion the eigenfunctions of the fundamental states, $\psi_{e,0}$ and $\psi_{h,0}$, and the corresponding eigen-energies, $E_{e,0}$ and $E_{h,0}$, can be calculated using Eq. (12) and (II.2). From the estimation of the eigenstates, the overlap factor $I_{0,0}$ of Eq. (30) can then be directly computed. To sum over all interband transitions, Weyl’s law can be used to estimate the system DOS instead of summing over all possible $i,j$ as in Eq. (29). In Sec. II.3 it was shown that in a 3D system exhibiting confinement along one direction and translational invariance along the others, Weyl’s law predicts an LDOS exactly identical to that of the bulk case. Although the values of $E_{e,0}$ and $E_{h,0}$ account for the quantum confinement of the carriers in the disordered potential, bulk asymptotic Weyl’s law is only a continuous yet good approximation of the discrete energy spectrum. Future works should allow us to provide a discrete and even better estimate of the spectrum in each localization subregion, based on $W$. To compute the JDOS of the homogeneous QW, the following well-known expression for bulk is thus used considering an effective band gap $E_{g}+E_{e,0}+E_{h,0}$:

[TABLE]

Finally the absorption coefficient of the homogeneous QW can be rewritten as

[TABLE]

where all quantities can be derived from the landscapes without solving the Schrödinger equation. (The factor 2/3 in Eq. (32) and the absence of the $1/L$ prefactor appearing in Eq. (29) are due to the use of the bulk expression of the JDOS in the absorption coefficient.Singh (2007))

The procedure for determining the absorption coefficient in the case of a disordered QW remains essentially the same. When superimposing the maps of the two landscapes, one can define the subregions which are the intersections between the various localization subregions for electrons and holes (Fig. 5). For each of these electron-hole “overlapping” subregions a local value of $\alpha$ can be computed using Eq. (32). The overall absorption coefficient of the disordered QW is obtained by summing over all electron and hole subregions, $\Omega_{m}$ and $\Omega_{n}$, as

[TABLE]

where the JDOS and overlap factors depend on the fundamental state and corresponding energy of the considered subregion, namely $\psi_{e,0}^{(m)}$ and $E_{e,0}^{(m)}$ for electrons, and $\psi_{h,0}^{(n)}$ and $E_{h,0}^{(n)}$ for holes, which again can be calculated from the landscapes using Eq. (12) and (II.2). Such an accurate accounting of carrier localization in a disordered material is crucial for a precise assessment of the below-gap absorption processes. A detailed study on this topic is presented in Ref. 19.

III.2.2 Emission

Let us first discuss the homogeneous case. While in the previous derivation of optical absorption [Eq. (32)] the conduction and valence bands were considered as completely empty and filled, respectively, in the case of optical emission the occupation of the states depends on the specific injection conditions of the QW. Assuming that the quasi-Fermi levels for electrons and holes are known, the electron and hole distributions $n(z)$ and $p(z)$ in the QW can be determined using Eq. (20). The radiative recombination rate then reads as

[TABLE]

corresponding to the transition energy $E_{g}+E_{e,0}+E_{h,0}$ calculated from the landscapes as in Sec. III.2.1. $B_{0}(z)$ is the local radiative recombination coefficient. If the transition is homogeneously broadened then a Lorentzian function must be used to determine the emission spectrum of the non disordered QW as

[TABLE]

In the case of a disordered QW, the emission spectra have to be summed over all possible transitions between localized states to produce the inhomogeneously broadened luminescence spectrum of the QW (see companion paper LL3, Ref. 20).

The localization landscapes can also be used to compute non-radiative recombination processes which may be strongly affected by disorder. In Auger recombination the confinement increases the carrier momentum compared to free carriers, leading to a better overlap with the wave function of the final high-energy carrier, and thus an enhanced transition probability. Let us consider for instance the “hhe” Auger process: an electron and a hole recombine and transfer their energy through Coulomb interaction to a second hole which becomes highly energetic (hot). Using the Fermi golden rule, one expresses the Auger recombination rate in one overlapping region:

[TABLE]

where $\rho(E_{f})$ is the density of final states and $M_{if}$ is the matrix element defined as:

[TABLE]

$r_{i}$ being the initial and final positions, and $V$ the Coulomb interaction potential. $\psi_{h,0}(r_{1})$ and $\psi_{h,0}(r_{2})$ are the initial states of the two holes, while $\psi_{e,0}(r_{1})$ is the initial electron state. Finally, $\varphi_{h,f}(r_{2})$ is the final state of the second hole. After dividing the entire system into the aforementioned overlapping regions, the above integral is computed over each region separately. For instance, in the above integral, the final hole wave function $\varphi_{h,f}$ is assumed to be a simple plane wave of wave vector $k_{f}$, while the initial states are approximated by the local landscapes in the considered overlapping region (see Section II.2). All these local integrals can then be assembled to build maps of the Auger recombination times in the whole system. The influence of compositional disorder on Auger recombination tackled by the localization landscape theory is currently under study.

IV Applications of the landscape theory to simple 1D heterostructures

In the following sections, we study several structures to show in various configurations the ability of our model to capture the main features of quantum devices. We successively study the eigen-energies and the overlap between quantum states in single QW (SQW) and 3-QW structures (Section IV.1). We then apply our model to the computation of density of electronic states in 1-D periodic and disordered superlattices (Section IV.2). Finally in Section IV.3, a 1-D SQW structure is simulated to perform a detailed comparison of the carrier distributions computed by the classical Poisson equation, the Poisson-landscape model (LL theory), and the Poisson-Schrödinger equation, respectively. All modeled structures presented in the following are 1-D with homogeneous layers. A full 3-D modeling of carrier transport in a semiconductor device including the effect of compositional disorder is presented in the companion paper LL3 (Ref. 20).

IV.1 Eigen-energies and overlap

In GaN-based materials, the polarization field at hetero-interfaces can strongly depend on the growth direction due to the wurtzite crystal structure. Conventionally, GaN grown along the (0001) direction ($c$-plane) cannot avoid a lattice mismatch which in turn induces a piezoelectric polarization. In addition, the relative displacement of the cation and anion sub-lattices from the ideal wurtzite position generates a net spontaneous polarization. In contrast, GaN grown along the (1$\bar{1}$00) direction ($m$-plane) is non-polar and does not exhibit any strain-induced piezoelectric polarization or spontaneous polarization field. In the simulation framework developed here, the density of polarization charges $\rho_{pol}$ is computed by taking the divergence of the total polarization ($\nabla\cdot$Ptotal), including spontaneous and piezoelectric polarization fields to account for the internal electric field at the hetero interfaces. Spontaneous and piezoelectric polarizations are computed using Eq. (39) and (40), while parameters of polarization values and piezoelectric coefficients can be found in Table 2 and 3, respectively (see Appendix).

First, we test the ability of the Poisson-LL model to accurately predict the energies and the spatial extension of quantum states in QW structures. Two different types of structures are simulated (Fig. 6). The first one is a 3 nm SQW of $m$-plane GaN enclosed between 50 nm barriers of AlxGa1-xN, several values of Al fraction $x$ being computed in a range from 0.1 to 0.5. The second structure is formed of 3 QWs of $c$-plane GaN of 3 nm width separated by 7 nm barriers of AlxGa1-xN, also enclosed between 50 nm barriers of AlxGa1-xN. This allows us to test the confinement effects in the wells in our model.

The LLs and the corresponding localization potentials are computed from the conduction and valence bands for both structures without external applied bias. Figure 7 displays the corresponding band structures, together with the effective potentials $W_{c}$ and $W_{v}$. Despite the sharp boundaries of the conduction and valence bands at the well-barrier interface, one can see that the variations of the potentials $W_{c}$ and $W_{v}$ extend much further than the width of the well. The fundamental energies for electrons and holes in the well are therefore computed from Eq. (II.2) on a larger integration region than the well itself. The integration boundaries (different for electron and hole) are represented by the dashed lines in Fig. 7. The boundaries of the integration region outside from the QW are first set to the inflection points of the local effective potential $W$. Then the integration boundaries are extended from these points by the characteristic decay length of a wave function of energy $E$ in a barrier of height $W$, i.e., $\hbar/\sqrt{2m(W-E)}$. The lighter effective mass of electrons therefore leads to a larger integration domain than the one of the holes which have a heavier effective mass. Eigen-energies as well as the overlap between electron and hole fundamental states $\psi_{e}$ and $\psi_{h}$ are computed on these integration regions, $\psi_{e}$ and $\psi_{h}$ being approximated by $u_{e}$ and $u_{h}$ normalized on the same region (cf. Eq. (12)).

Figures 8 and 9 display the comparisons between the fundamental energies computed by solving the Schrödinger equation directly, and the energies computed using Eq. (II.2), for $m$-plane and $c$-plane cases with various values of the aluminum content. In both figures, the top frames exhibit a remarkable agreement between the two different computations. The bottom left frames show the differences between the electron and hole energy, in other terms the energy of the smallest radiative transition. The bottom right frames display the value of the overlap integral between electron and hole fundamental states (in the case of the 3-QW structure, the states are taken in the central quantum well) (dashed line), compared with the estimates obtained using Eq. (11) for the wave functions. The overlap region used for the calculation are the intersections of the regions $\Omega_{e}$ and $\Omega_{h}$ displayed in Fig. 7. Although one observes here a slight deviation by a few percents of the approximated value using the landscapes, the agreement, both in absolute value and in trend, remains very good in the two structures, for all values of the Al content.

IV.2 Density of states, overlap and absorption computation

IV.2.1 Periodic superlattice

AlGaN/GaN superlattice (SL) structures are widely used in commercial LEDs to prevent the electron current leakage and improve lateral current spreading.Han et al. (2009); Song et al. (2013) However, modeling SL structures still remains a challenge for classical Poisson-DD solvers. In the classical picture, the resistance experienced by carriers is determined by the barrier height. In SLs, this leads to an overestimation of its value. Schrödinger-based solvers can model SLs, accounting for wave function coupling and tunneling effects. However, this approach cannot be applied in the case of multidimensional devices to study current crowding effects of disordered systems, due to the high demand of computation time. We show here that the LL theory allows us to overcome this constraint.

We model a 20-pairs $n$-type $m$-plane Al0.4Ga0.6N/GaN SL structure with flat band conditions at both ends, as shown in Fig. 10(a), to analyze the electron transport behavior. Two structures with different periodicities, 5 nm/5 nm and 1 nm/1 nm QW/barrier thickness, are simulated. A nonpolar $m$-plane orientation is considered at first to study the intrinsic transport properties of SLs as described by the landscape model. In the thicker SL (5 nm/5 nm), quantum effects are weaker and $W_{c}$ deviates only slightly from $E_{c}$, as shown in Fig. 10(b). In contrast, in the thinner SL (1 nm/1 nm) $W_{c}$ is significantly different from $E_{c}$, as shown in Fig. 10(c). This difference is a manifestation of the coupling between wells when the barriers are thin, and expresses how the effective potential is able to translates quantum tunneling into a shift of the effective conduction band edges.

To evaluate the effective barrier seen by carriers in a complicated structure such as SLs, Weyl’s asymptotic formula is used Arnold et al. (2016) to obtain the 1D LIDOS [see Eq. (II.3.2)], which is then integrated over the entire system:

[TABLE]

where IDOS$(E)$ is the integrated density of available states. As seen in Section II.1, $W_{c}(z)=1/u_{e}(z)$ can be understood as the effective conduction band edge.

For large SL barrier and QW thicknesses, quantum effects are weak and classical. Therefore original potential- and landscape-based models give very similar estimates of the real IDOS (blue and black lines), as displayed in Fig. 11(a). One can observe that the take-off energy is better approximated by the landscape model (continuous black line) when compared to the computation from Schrödinger equation (red dashed line). In this situation of very weak coupling between wells, the similarity of all wells (no disorder here) induces a strong energy degeneracy which appears in the real IDOS. This IDOS increases through large steps of height 20 (the number of wells) which underlines the discrete nature of the spectrum, while the landscape-based (or $W_{c}$-based) approximation shows a continuous line following these steps.

As the SL thicknesses of QWs and barriers becomes smaller, wells become coupled, the degeneracy is lifted, creating subbands in the entire structure, and the corresponding IDOS (red dashed line) proceeds by smaller steps of height 1 which are almost invisible at the scale of the figure. The landscape-based IDOS ($W_{c}$-based Weyl’s law) shows here an excellent agreement with the exact calculation from the Schrödinger equation [Fig. 11(b)], showing that it takes well into account the coupling between wells, in other words, the tunneling phenomenon. The approximation based on the original potential $E_{c}$, however, falls very far off the true IDOS, especially for lower energy states. Going into finer detail, one can also observe that in this last case, the $W$ approximation of the IDOS becomes positive about 50 meV earlier than the actual IDOS of the Schrödinger equation. This slight discrepancy comes from the fact that, by definition, the $W$ approximation of the IDOS is continuous while the actual IDOS is step-wise. Therefore, to reach a value close to 1 at an energy corresponding to the fundamental energy of the system, the approximate IDOS has to take off and be positive at a smaller energy.

To test further the LL model, we analyze more complicated structures such as 20-pairs $n$-type SL with double periods composed of 10-pairs 1 nm/1 nm and 10-pairs 2 nm/2 nm SLs, respectively, including $m$-plane and $c$-plane cases. Figures 12(a) and 12(c) show the conduction band edge $E_{c}$ and $W_{c}$ for $m$-plane and $c$-plane conditions. Here again, the difference between $E_{c}$ and $W_{c}$ increases as the period become smaller, exhibiting stronger tunneling effects. The IDOS displayed in Figs 12(b) and 12(d) shows the quality of the approximation provided by $W_{c}$, regardless of the structure complexity. We found the $c$-plane case exhibits a better agreement with the Schrödinger model due to distinct potential energy distributions. As randomness and inter-coupling increase in the system, the prediction of $W_{c}$ becomes closer to the solution of the Schrödinger equation. This last result was in fact already partly demonstrated in Ref. 18 for very random systems.

IV.2.2 Disordered superlattice

To illustrate the efficiency of the landscape model in the case of disordered systems, we compute the band structure and the density of states in nitride $m$-plane and $c$-plane disordered SL. Two types of SL are investigated: the first is composed of 20 pairs of well/barrier layers, while the second is composed of 200 pairs. In both cases, wells and barriers have 1-nm thickness. The barrier material is Al0.4Ga0.6N and the well material is AlxGa1-xN, where $x$ is randomly and independently determined in each well, using a uniform law between 0 and 0.4.

For each type of SL, the band structure and the IDOS are computed (Figs 13 and 14). In both cases, although the edge of the conduction band $E_{c}$ now exhibits large fluctuations across the structure due to the compositional disorder, the quantum coupling between wells translates into a much smoother effective potential $W_{c}$. The value of $W_{c}$ is significantly larger than $E_{c}$ and can be interpreted as a local fundamental energy from the expression of the local density of states in Eq. (22). Looking at the IDOS (right column of Figs 13 and 14), we observe a very good agreement between the actual counting functions and its approximation obtained using $W_{c}$. One has to note the same slight discrepancy of the take-off energy between both functions that was observed in the periodic case. This discrepancy, which comes from the different natures of the two functions (continuous versus stepwise), appears to be of the same order of magnitude independently on the number of wells in the structure.

We also use the $c$-plane 200-pairs disordered SL structure (bottom of Fig. 14) to test the quality of the local integrated density of states approximation. To this end, the SL structure is divided into four disjoint regions of same length. The local integrated density of states (LIDOS) of each region is computed in two different ways: first, by directly solving the Schrödinger equation in the entire system, each quantum state is then assigned a unique region among the four already defined, which is the region where the state reaches its maximum amplitude. Secondly, by integrating the LDOS of Eq. (22) on each regio, comparisons of the two methods for all four regions are presented in the four bottom frames of Fig. 15. Here also one can observe the accuracy of the LIDOS computed using the landscape approach.

We finally test the accuracy of the overlap and joint density of states estimates. To this end, the absorption spectrum between electrons and heavy holes (of respective masses 0.20 and 1.87) of both disordered SL (20 and 200 wells) are computing summing over all localization subregions the formula of Eq. (33), and compared to exact computations using the quantum states computed by Schrödinger equation. Comparisons between the two types of calculations are displayed in Fig. 16, both in linear and log scale. One can see that the absorption curves obtained using the LL (plain black lines) are very close to the ones computed using Schrödinger equation (red dashed lines), over almost three decades of absorption rate, confirming the quality of the landscape approach. Our results show that $W$ appears as a very likely candidate for inserting into a model of carrier transport able to account for quantum effects in complicated systems.

IV.3 Carrier distribution

To study and compare in detail the carrier distributions predicted by the different models, a 1D SQW is simulated. The structure is composed of a 3-nm active GaN layer, enclosed between two barriers of AlxGa1-xN, where $x=0.2$ will be the reference case. The band offsets between GaN/AlGaN conduction bands are assumed to be 63$\%$ of the band-gap discontinuity. The detailed band structure parameters for GaN and AlN are provided in Table 1 (see Appendix). All parameters of AlxGa1-xN alloys are obtained by interpolation, where the bandgap alloy bowing parameter is assumed to be 0.8 eV.

The $m$-plane and $c$-plane cases (without and with the polarization charge induced at the interface) are both discussed. The detailed dimension and material doping level are shown in Fig. 17(a), where the electron effective mass of Al0.2Ga0.8N and GaN is 0.214m0 and 0.20m0, respectively. The doping is assumed to be fully activated to ensure the same activation condition for comparison.

The carrier density distribution is solved through three different methods: (1) the classical Poisson model solving Poisson equation for the charge distribution of ionized donors across the QW structure; (2) self-consistent Poisson-LL approach; and (3) self-consistent Poisson-Schrödinger approach. Figure 17(b) displays the computed carrier distribution and the potential energy for the $m$-plane case with a symmetric potential energy by the three methods. We can see that the carrier distribution obtained from the classical Poisson equation (blue line) is almost constant across the well and dropping sharply outside, not accounting for the quantum nature of the electron. The Poisson-Schrödinger solver provides a very different outcome: the carrier distribution is more confined in the center of the well [red line of Fig. 17(b)] and extends smoothly outside of it, revealing the wave-function shape of the carrier. Turning to the landscape model, the effective quantum potential $W_{c}$ (black line) exhibits a smoother behavior than the original conduction band (blue line). The carrier density computed after replacing the original conduction band edge $E_{c}$ with $W_{c}$ appears very similar to that computed from the Schrödinger equation.

In Fig. 17(c), we evaluate the $c$-plane case with an asymmetric potential profile induced by the polarization charge. The carrier distribution computed by the classical Poisson equation is sharp and mostly located at the minimum of the conduction band. This does not match the result obtained from the Poisson-Schrodinger solver. This excessively large carrier density located at the interface might be the reason why some numerical studiesKuo et al. (2010) using Poisson and DD solvers to study GaN-based polar QWs adopt only $\sim$50% of the theoretical polarization charge, in order to reproduce a carrier distribution consistent with the Schrödinger equation. Here again, the Poisson-LL model results in a smoother carrier distribution, as shown in Fig. 17(c), much closer to the Schrödinger solution.

In the case of a deeper $m$-plane QW with $x=0.4$, carrier distributions computed with the Poisson-LL and Poisson-Schrödinger models are also in good agreement, as displayed in Fig. 17(d). To quantify this agreement, we simulate a series of Al composition from 10% to 50% ($0.1\leq x\leq 0.5$). For each simulation, we compute the dimensionless quantity $\Delta$ defined as the relative difference of peak carrier densities between the Poisson-landscape and Poisson-Schrödinger models. As displayed in Fig. 18, $\Delta$ becomes smaller when the carriers are well localized within a deeper potential (larger Al composition) and almost vanishes above 40% Al composition.

We can conclude from this that the LL model matches very well the solution of the Schrödinger equation when the system is strongly localized. Even in systems with a lower degree of localization (here, small $x$) the prediction of the Poisson-LL model gives an acceptable agreement with the exact Poisson-Schrödinger solution and provides an overall description of the carrier distribution much more accurate than the result of a classical Poisson solver. The discrepancy observed between the carrier distributions computed in Poisson-LL and Poisson-Schrödinger models can be attributed to the fact that the local DOS used in Poisson-LL does not fully describe the exact shape of the wave function. If required, an even more accurate approximation might be achieved by using the effective confining potential to compute a correction to this carrier distribution without any adjustable parameter.Ancona and Tiersten (1987); Ancona and Iafrate (1989)

V Conclusion

In this work, we applied the LL theory,Filoche and Mayboroda (2012); Arnold et al. (2016) until now a purely mathematical framework, to build a new model of quantum and disordered semiconductors devices. In the LL theory, the Schrödinger equation is replaced by an associated Dirichlet equation whose solution is called the localization landscape, $u$. The reciprocal of this landscape, $W=1/u$, acts as an effective quantum classical confining potential which governs localization of quantum states. Eigenfunction profiles and geometry of localization subregions can be retrieved from a direct analysis of this effective potential. Inserting $W$ into Weyl’s law also provides a very good approximation of the density of states, the carrier concentrations, and the spatial distributions of charges, especially in the case of strong localization. Besides, as seen here in 1D cases, and shown to play a decisive role in the 3D modeling of LEDs (see LL3, Ref. 20), the LL theory simulates to an excellent approximation two major effects of quantum mechanics, namely, the reduction of barrier heights (tunneling effect) and the raising of energy ground states (quantum confinement effect).

We have presented here how to compute localized states, energies, density of states using the landscapes, and how to couple them to Poisson and DD equations to model carrier transport in SL structures and carrier localization in nitride-based systems. In principle, this method is not restricted to modeling nitride-based devices (as the examples presented here or in companion papers dealing with the simulation of the absorption edge (LL2, Ref. 19) and LEDs (LL3, Ref. 20) based on InGaN alloy materials), but can expand to other semiconductor materials and to any electronic or optoelectronic properties requiring the knowledge of electron and hole quantum states. Moreover, except for the transport and carrier distribution issues discussed in this paper, the emission and absorption in the disordered system can be properly modeled in terms of eigen-energy calculation (LL2, Ref. 19). Finally, as observed in real modeling exercises (LL3, Ref. 20), the computation time using the landscape model is considerably reduced compared to a conventional Schrödinger solver, which makes this model ideal for simulating and designing quantum 3D real-world devices.

Appendix A Band structures and polarization in AlGaN alloys

Table 1 displays the detailed band structure parameters for GaN and AlN, respectively. The parameters of all AlxGa1-xN alloys are obtained by interpolation, the bandgap alloy bowing parameter being taken equal to 0.8 eV. Spontaneous and piezoelectric polarizations are computed using Eqs. (39) and (40), while parameters of polarization values and piezoelectric coefficients are provided in Tables 2 and 3, respectively.

[TABLE]

[TABLE]

Acknowledgements.

This work was supported by the project CRIPRONI (ANR-14-CE05-0048-01/MOST-104-2923-E-002-004-MY3) of the French National Research Agency (ANR) and Taiwanese Ministry of Science and Technology (MOST). 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 and Claude Weisbuch was provided by the DOE Solid State Lighting Program under Award #DE-EE0007096.