to 1-D simulations, since the extension to a 3-D system is computationally basically unfeasible. Instead of solving the large eigenvalue problem connected to evaluating the Schrödinger equation, we have implemented quantum corrections via localization landscape theory (LLT) [132]; this approach is numerically much more efficient and gives results similar to a full self-consistent Schrödinger-Poisson solver, as we will also discuss below. From LLT we extract an effective confining potential for the conduction and valence band edge starting from the TB energy landscape.

An example of the resulting quantum corrected energy landscape is given in Fig. B.3 (a) and (b). Here, in-plane band edge profiles for a single atomic plane through an In<sub>0.1</sub>Ga<sub>0.9</sub>N QW, after LLT has been applied, are shown. As Fig. B.3 (a) reveals, the fluctuations in the valence band edge energy due to alloy fluctuations are of the order of 100 meV. In combination with the high effective hole mass, these fluctuations are large enough to give rise to strong carrier localization effects as seen in other studies already [260, 242, 102]. We therefore expect that, especially for holes, the inclusion of random alloy fluctuations in the simulation will impact the carrier distribution. Consequently recombination rates are also expected to be noticeably influenced.

The variation in the conduction band edge energy is significantly smaller (order of 30 meV), as can be seen in Figure B.3 (b). Since the effective electron mass is much lower in comparison with the holes, electron wave functions are less strongly perturbed by alloy fluctuations. The impact that these fluctuations in the band edge energies have on the radiative recombination is also seen in Fig. B.3 (c); the radiative recombination is calculated with *ddfermi* as will be described in section B.3.2. The correlation between the valence band edge maxima and regions of high radiative recombination can be clearly identified; similar spatial profiles can be seen for non-radiative (Auger) recombination (not shown).

In order to highlight the impact of random alloy fluctuations on carrier transport and the distribution of carriers across a MQW system, we compare our atomistic calculations with the outcome of a VCA. In the latter case a homogeneous effective crystal is constructed where material properties are chosen to be interpolated properties of the binaries InN and GaN within the QW region. Here a linear, composition weighted interpolation scheme is employed. A bandgap bowing of -2.0 eV is used, consistent with the underlying atomistic TB model [60]. The VCA description, without any Gaussian broadening, is similar to commercially available packages. However, and in contrast to commercial software packages, quantum corrections via LLT can also be taken into account in our VCA simulations, following the approach used for the random alloy case.

## B.3.2 Device simulation

Having outlined above the generation of the energy landscape of the active region, e.g. the (In,Ga)N/GaN MQW system, a full device mesh, including the *n*- and *p*-doped regions, needs to be constructed on which the DD equations are solved. To achieve this we proceed as follows and divide the device mesh into two regions: an atomistic and a macroscopic one. The atomistic region is used to describe the MQW region and has as many grid points as atoms in the system. These points contain information about the conduction and valence band edge energies calculated from TB, as discussed above.

In order to capture the effects of carrier localization in the calculations, the in-plane dimensions of our 3-D simulation cell should be larger than the localization length of the holes, given that electrons are less strongly affected by alloy fluctuations [252]. In our atomistic calculations we use a system with in-plane dimensions of 12.8 × 11.0 nm². This is large enough to see the effects of hole localization as the in-plane hole localization length for In<sub>0.1</sub>Ga<sub>0.9</sub>N QWs is of
the order of 1 nm [253]. The in-plane dimensions can be seen in Fig. B.3 (a) and (b) where the in-plane valence and conduction band edges of an In$_{0.1}$Ga$_{0.9}$N QW are shown. In case of the VCA, given that there are no variations in material properties (band edge energies) within the growth plane (c-plane), a much smaller in-plane area is sufficient (1.3 × 1.1 nm$^2$), which reduces the numerical effort.

LLT is solved on this (finite element) mesh using a finite element method (FEM). The DD calculations are carried out employing a Voronoi finite volume method (FVM) [128]. Therefore, the generated FEM mesh must be transferred to an appropriate FVM mesh. Every point from the FEM mesh is included on the FVM mesh, as well as extra points required to produce a boundary-conforming Delaunay tetrahedral mesh; conduction and valence band data are then interpolated onto these additional nodes. This mesh is then embedded within a macroscopic device mesh which contains information about the p- and n-doped regions. In our atomistic transport studies here, we focus on systems without an (Al,Ga)N electron blocking layer (EBL). In principle, an atomistic mesh resolution would be required for the EBL, given the alloy fluctuations in (Al,Ga)N. However, and as we will discuss below, an (Al,Ga)N EBL is of secondary importance for the questions targeted in the present study. Therefore, outside active MQW region of the system, pure GaN is assumed. Thus, the conduction and valence band edge values are position independent (except for changes due an applied bias). The absence of strongly fluctuating band edges in the macroscopic mesh region allows us to use a sparse mesh and scale the simulation to a full device. The mesh is created using TetGen [247] and the interpolation is handled via WIAS-pdelib [141]. More details on the mesh generation can be found in Ref. [221].

Equipped with knowledge about the mesh generation, we turn now to the DD simulations. To do so we build on the van Roosbroeck system of equations [258]:

$$-\nabla \cdot (\varepsilon_s(\mathbf{r})\nabla\psi(\mathbf{r})) = q(p(\mathbf{r}) - n(\mathbf{r}) + C(\mathbf{r})), \qquad (B.3a)$$

$$\nabla \cdot \mathbf{j}_n(\mathbf{r}) = qR(\mathbf{r}), \qquad (B.3b)$$

$$\nabla \cdot \mathbf{j}_p(\mathbf{r}) = -qR(\mathbf{r}), \qquad (B.3c)$$

$$\mathbf{j}_n(\mathbf{r}) = -q\mu_n n(\mathbf{r})\nabla\varphi_n(\mathbf{r}), \qquad (B.3d)$$

$$\mathbf{j}_p(\mathbf{r}) = -q\mu_p p(\mathbf{r})\nabla\varphi_p(\mathbf{r}). \qquad (B.3e)$$

In the above equations, $q$ is the elementary charge, $\varepsilon_s(\mathbf{r}) = \varepsilon_0\varepsilon_r(\mathbf{r})$ is the dielectric permittivity, $\psi(\mathbf{r})$ is the device electrostatic potential, $p(\mathbf{r})$ and $n(\mathbf{r})$ are the hole and electron densities, $C(\mathbf{r}) = N_D^+(\mathbf{r}) - N_A^+(\mathbf{r})$ is the net activated dopant density, $\mathbf{j}_n(\mathbf{r})$, $\mathbf{j}_p(\mathbf{r})$, $\varphi_n(\mathbf{r})$ and $\varphi_p(\mathbf{r})$ are the electron and hole current densities and the respective quasi-Fermi potentials. The total recombination rate is denoted by $R$.

The carrier densities are related to the band edge energies, $E_c^{dd}(\mathbf{r})$ and $E_v^{dd}(\mathbf{r})$, and quasi-Fermi potentials, $\varphi_n(\mathbf{r})$ and $\varphi_p(\mathbf{r})$, via the state equations [128]

$$n(\mathbf{r}) = N_c \mathcal{F} \left( \frac{q(\psi(\mathbf{r}) - \varphi_n(\mathbf{r})) - E_c^{dd}(\mathbf{r})}{k_B T} \right), \qquad (B.4a)$$

$$p(\mathbf{r}) = N_v \mathcal{F} \left( \frac{E_v^{dd}(\mathbf{r}) - q(\psi(\mathbf{r}) - \varphi_p(\mathbf{r}))}{k_B T} \right). \qquad (B.4b)$$

Here, $N_c$ and $N_v$ are the effective density of states for the conduction and valence band, respectively. For the distribution function, $\mathcal{F}$, we use Fermi-Dirac statistics, and $k_B$ is the Boltzmann constant. A temperature of $T = 300$ K has been used in all calculations. It is to
note that the valence band, $E_{n}^{dd}(\mathbf{r})$, and conduction band edge energy, $E_{c}^{dd}(\mathbf{r})$, can be described either by a VCA, a VCA including quantum corrections via LLT, or an atomistic random alloy calculation including LLT-based quantum corrections.

The total recombination rate in Eqs. (B.3b) and (B.3c), $R$, is calculated using the ABC model [178, 226]. Here, $R$ is the sum of (defect related) Shockley-Read-Hall, $R^{SRH}$, radiative, $R^{RAD}$ and (non-radiative) Auger recombination rate, $R^{AUG}$. The SRH rate is obtained from:

$$R^{SRH}(\mathbf{r}) = \frac{r(n,p)}{\tau_p(n(\mathbf{r}) + n_i(\mathbf{r})) + \tau_n(p(\mathbf{r}) + n_i(\mathbf{r}))}, \qquad (B.5)$$

the radiative part via

$$R^{RAD}(\mathbf{r}) = B_0 r(n,p), \qquad (B.6)$$

and the Auger rate is calculated as

$$R^{AUG}(\mathbf{r}) = (C_n n(\mathbf{r}) + C_p p(\mathbf{r}))r(n,p). \qquad (B.7)$$

In Eqs. (B.5) to (B.7)

$$r(n,p) = n(\mathbf{r})p(\mathbf{r}) - n_i^2(\mathbf{r})$$

and

$$n_i^2(\mathbf{r}) = n(\mathbf{r})p(\mathbf{r}) \exp \left( \frac{q\varphi_n - q\varphi_p}{k_B T} \right).$$

The above equations require further input, namely the radiative recombination coefficient $B_0$, the Auger recombination coefficients $C_p$ and $C_n$ as well as the the SRH lifetimes $\tau_p$ and $\tau_n$. All these parameters will in principle carry a composition dependence [210, 211, 181]. Furthermore, $B_0$, $C_p$ and $C_n$ will also be carrier density dependent [93, 94, 169]. We follow here the widely made assumption that these coefficients are constant across the InGaN MQW region [194, 203]. In the following we take a weighted average of parameters calculated in Ref. [210] for an electron and hole density of $3.8 \times 10^{18} \text{ cm}^{-3}$, which is a good approximation for the average carrier densities in the QWs at a current density of $50 \text{ A/cm}^2$. As our active region consists of two In$_{0.1}$Ga$_{0.9}$N QWs and one In$_{0.125}$Ga$_{0.875}$N QW we evaluate the different recombination coefficients as follows:

$$R_i^{\text{eff}} = \frac{2 \times (R_i^{10\%}) + 0.5 \times (R_i^{15\%} + R_i^{10\%})}{3}. \qquad (B.8)$$

Here, $R_i \in \{B_0, C_n, C_p\}$ are the radiative recombination, electron-electron-hole and hole-hole-electron Auger recombination coefficients, respectively. As there are no values for an In$_{0.125}$Ga$_{0.875}$N QW in Ref. [210], a linear average of the coefficients in In$_{0.1}$Ga$_{0.9}$N and In$_{0.15}$Ga$_{0.85}$N wells has been used. A summary of the material parameters employed in all simulations is given in Table B.1.

The numerical approximation of the van Roosbroeck system is implemented (in 3-D) in ddfermi [104]. As already mentioned above, we employ the FVM and the current is discretized using the SEDAN (excess chemical potential) approach [265, 50, 2], which yields a thermodynamically consistent flux approximation in the sense of Ref. [128].

To simulate the devices under study, we also used the commercial software nextnano [30], which relies on the simulation of a self-consistent Schrödinger-Poisson-DD system. In this work we use nextnano to simulate the carrier transport in the above discussed MQW systems within a 1-D approximation. In nextnano we utilize the same parameter set as in the ddfermi simulations. Therefore, the obtained results can be directly compared to our 3-D VCA model.
<table><caption>Table B.1 – Material parameters used in the different regions of the simulation supercell. Parameters denoted with † are taken from [194]; parameters denoted with ‡ are derived from [210] as described in the main text.</caption><thead><tr><th colspan="2">Parameter</th><th colspan="3">Value in each region</th></tr><tr><th>Name</th><th>Units</th><th>p-GaN</th><th>i-InGaN</th><th>n-GaN</th></tr></thead><tbody><tr><td>Doping</td><td>cm<sup>-3</sup></td><td>5 &times; 10<sup>18</sup></td><td>1 &times; 10<sup>16</sup></td><td>5 &times; 10<sup>18</sup></td></tr><tr><td>&mu;<sub>h</sub></td><td>cm<sup>2</sup>/Vs</td><td>5</td><td>10</td><td>23</td></tr><tr><td>&mu;<sub>e</sub></td><td>cm<sup>2</sup>/Vs</td><td>32</td><td>300</td><td>200</td></tr><tr><td>&tau;<sub>p</sub></td><td>s</td><td>10</td><td>1 &times; 10<sup>-7</sup></td><td>7 &times; 10<sup>-10</sup></td></tr><tr><td>&tau;<sub>n</sub></td><td>s</td><td>6 &times; 10<sup>-10</sup></td><td>1 &times; 10<sup>-7</sup></td><td>10</td></tr><tr><td>B<sub>0</sub></td><td>cm<sup>3</sup>/s</td><td>2.8 &times; 10<sup>-11</sup></td><td>2.8 &times; 10<sup>-11</sup></td><td>2.8 &times; 10<sup>-11</sup></td></tr><tr><td>C<sub>p</sub></td><td>cm<sup>6</sup>/s</td><td>5.7 &times; 10<sup>-30</sup></td><td>5.7 &times; 10<sup>-30</sup></td><td>5.7 &times; 10<sup>-30</sup></td></tr><tr><td>C<sub>n</sub></td><td>cm<sup>6</sup>/s</td><td>1 &times; 10<sup>-31</sup></td><td>1 &times; 10<sup>-31</sup></td><td>1 &times; 10<sup>-31</sup></td></tr></tbody></table>

When including quantum corrections in ddfermi, LLT is used. In nextnano a self-consistent Schrödinger-Poisson-DD calculation is performed where a k·p Hamiltonian is used to calculate eigenstates across the full simulation domain. Following the ddfermi set up, in nextnano we employ also a 1-band model for the calculation of the electron and hole densities.

## B.4 Results

In this section we present the results of our study on the carrier distribution in the above described (In,Ga)N/GaN MQW systems. To understand the impact of the alloy microstructure on the carrier distribution, in Section B.4.1 we start with standard 1-D calculations building on the commercial software package nextnano [30]. We use this entirely continuum-based description of the QWs also to determine the impact (i) of an EBL and (ii) a self-consistent Schrödinger-Poisson-DD treatment on the transport properties. Moreover, and as already mentioned above, (ii) can also be used as a benchmark for our 3-D ddfermi solver. In Section B.4.2 we then proceed to investigate the influence of random alloy fluctuations on the carrier distribution in the (In,Ga)N/GaN MQW stack.

### B.4.1 Continuum-based simulations of the carrier transport in (In,Ga)N-based LEDs

To examine the impact of random alloy fluctuations on the carrier distribution in an (In,Ga)N/GaN MQW stack, we start with a 'standard' 1-D simulation approach that is widely applied in the literature. In a first step we begin with nextnano calculations and as outlined above, compare the results to our ddfermi data.

#### nextnano simulations

To study how the presence of an EBL affects the ratio of radiative recombination $\rho$, Eq. (B.1), the systems outlined in Section B.2 are simulated with and without a 20 nm Al$_{0.15}$Ga$_{0.85}$N EBL using nextnano. The EBL is separated from the p-side QW (position 3 (p-side) in Fig. B.2) by a 10 nm GaN barrier. Similar settings for an (Al,Ga)N EBL have been used in previous studies [194]. The nextnano calculated ratio of radiative recombination $\rho$, when varying the position of the deep QW in the MQW stack, are depicted in Fig. B.4 (a). Turning first to the data without
<figure><img src="image_52.png" /></figure>

Figure B.4 – Ratio of radiative recombination $\rho$, Eq. (B.1), from the shallow wells ($\text{In}_{0.1}\text{Ga}_{0.9}\text{N}$) to recombination from the deep well ($\text{In}_{0.1}\text{Ga}_{0.9}\text{N}$) calculated as a function of the position of the deep well in the multi-quantum well stack. Here $\rho$ is evaluated using (a) nextnano excluding (purple) and including (green) quantum corrections via a self-consistent Schrödinger-Poisson-drift diffusion solver; results are shown when excluding (solid, filled circles) and including (dotted, open circles) an $\text{Al}_{0.15}\text{Ga}_{0.85}\text{N}$ blocking layer, and (b) ddfermi excluding (purple), including (green) quantum corrections via localization landscape theory (LLT) using a virtual crystal approximation (VCA) and a random alloy calculation including LLT-based quantum corrections (blue); these calculations neglect the AlGaN blocking layer.

quantum corrections, we find that in the case of the employed 1-D VCA-like continuum-based description, $\rho$ is small when the deep QW is at the n-side (position 1 (n-side) Fig. B.2) and larger when the deep well is at the p-side (position 3 (p-side) Fig. B.2). Thus, the 1-D model predicts the opposite trend when compared to experiment [145]. This trend is only slightly changed when including quantum corrections via a self-consistent Schrödinger-Poisson-DD model. In this case a much weaker dependence of the results on the position of the deep QW in the MQW stack is observed. However, even when including quantum corrections, the nextnano results for $\rho$ are not reflecting the experimentally observed behavior (see discussion above). Figure B.4 (a) reveals also that qualitatively the results do not depend on the presence of the EBL, indicating that for the structures considered, this feature of an LED is of secondary importance for the aims of this work.

### ddfermi simulations

Since we are also able to use the atomistic framework in a VCA setting, we compare our ddfermi results, cf. Fig. B.4 (b) (purple), with those from nextnano, cf. Fig. B.4 (a) (purple, solid). We focus on structures which neglect the EBL as we have found above that it does not impact results in a VCA. In both nextnano and ddfermi a similar trend is found: the deep QW dominates recombination only when it is located at the n-side. This is illustrated further in Fig. B.5 (a), which displays the contribution (in percent) to the radiative recombination rate from each QW (colors) in the MQW stack. The data are shown as a function of position of the deep QW in the MQW system. This confirms that it is always the QW which is closest to the n-doped side (position 1) that dominates the recombination process; the n-side QW contributes $\approx 95\%$ when the deep QW is at position 1, $\approx 70\%$ when the deep QW is at position 2 or 3. Again, we stress that this is the opposite trend to the experimental findings in Ref. [145].
<figure><img src="image_53.png" /></figure>

Figure B.5 – Contribution of each quantum well (n-side; centre; p-side) in the (In,Ga)N multi-quantum well system to the total radiative recombination $R_{\Omega_i}^{RAD}$ for $i \in \{n-side, center, p-side\}$ as a percentage of the total radiative recombination from all 3 quantum wells for (a) virtual crystal approximation (VCA), (b) virtual crystal approximation with quantum corrections included via localization landscape theory (VCA + LLT) and (c) a random alloy calculation including localization landscape theory based quantum corrections (Random alloy + LLT). That data are shown as a function of the position of the deep quantum well (x-axis). Each bar contains the percentage recombination from the n-side quantum well (purple), the center quantum well (green) and the p-side quantum well (blue). Labelling is consistent with that introduced in Fig. B.2.

To shed more light on this result, the upper row in Figure B.6 depicts the average hole (black, solid), electron (black, dashed) and radiative recombination (red) rate along the c-axis when the deep QW (In$_{0.125}$Ga$_{0.875}$N well) is (a) closest to the n-side (position 1), (b) in the centre of the MQW stack (position 2) and (c) closest to the p-side (position 3). Focusing on the VCA data, Figs. B.6 (i) (a-c), we see the cause of the dominant recombination from the n-side QW: the hole density is always high in this region, independent of which well is closest to the n-side. In particular, the p-side QW fails to capture holes effectively and consistently has the lowest hole density. We note that a similar behavior is also found in the nextnano calculations discussed in Sec. B.4.1.

Given that our VCA ddfermi approach and nextnano treat (In,Ga)N as a homogeneous alloy that can be described by averaged material parameters which do not vary throughout the wells (no alloy fluctuations included), it allows us also to compare the implemented methods for quantum corrections in DD simulations. Here, as discussed above, nextnano builds on the widely used Schrödinger-Poisson-DD model while ddfermi utilizes the recently developed LLT method. It has been discussed and shown in the literature that the LLT method can produce results in good agreement with the solution of the Schrödinger equation in the case of a 1-D effective mass approximation [132, 74]. Looking at Fig. B.4 (a) (green, solid) and Fig. B.4 (b) (green) we see that the results from our in-house developed ddfermi-based 3-D model, which employs LLT (3-D, ddfermi), are very similar to the standard self-consistent 1-D Schrödinger-Poisson-DD calculation underlying nextnano. This gives confidence that our LLT treatment is providing a comparable description of the quantum corrections in the system.

Overall, Fig. B.4 (a) reveals that when including quantum corrections in the VCA calculations, the position of the deep QW has little impact on the ratio of the relative radiative recombination, $\rho$. From Fig. B.5 (b) one can also gain more insight into this behavior and how quantum
corrections impact the carrier distribution in the MQW stack. In the absence of quantum corrections but utilizing a VCA, Fig. B.5 (a), the well closest to the *n*-side dominates the relative radiative recombination ratio *ρ* independent of the position of the deep well in MQW systems. When including quantum corrections this situation is now changed: the deep QW is now the dominant emitter independent of its position in the MQW stack.

This behavior becomes clear when analyzing the electron and hole densities as a function of the position of the deep well in the (In,Ga)N/GaN MQWs, as shown in Fig. B.6 (ii). Looking at the electron densities first, we find that electrons preferentially occupy the well closest to the *p*-side. This effect is enhanced when the deep QW is closest to the *p*-side (cf. Fig. B.6 (ii) (c)). In our previous study on uni-polar electron transport [221], we have already seen that including quantum corrections leads to a softening of the potential barrier at the QW barrier interfaces. This in turn can lead to an increased electron current at a fixed bias point, when compared to a VCA system without LLT treatment, and thus the electrons can more easily ‘overshoot’ the wells in the MQW system. As a consequence, a lower electron density in the well closest to the *n*-side is observed. Turning to the hole density, the situation is different. Here, we find that holes preferentially populate the well closest to the *n*-side. Only when the deep QW is closest to the *p*-side, the hole density in this well is noticeably increased. However, when comparing the distribution of holes in the MQW as a function of the position of deep well in absence (Fig. B.6 (i)) and presence (Fig. B.6 (ii)) of quantum corrections, the results are not very different. This indicates that quantum corrections, at least when employing a VCA, are of secondary importance for the hole distribution. This finding is consistent with our previous results on uni-polar hole transport [222], where we have discussed that due to the high effective hole mass and the small valence band offset, quantum corrections have a smaller impact on the hole transport when compared to electrons. As a consequence, the distribution of holes follows a similar pattern to that of the VCA where quantum corrections are neglected. Finally, when looking at the ratio of radiative recombination *ρ*, it is important to note that this quantity is not only determined by having both large electron and hole densities in the same well but also by their spatial overlap. As one can infer from Fig. B.6 (ii) (a-c), the largest radiative recombination rate is always observed in the deepest well. This indicates also that the spatial overlap of electron and hole densities largest in the deep QW regardless of its position across the MQW system. We stress again that even when including quantum corrections in the VCA calculations, the resulting trend in *ρ* is not reflecting the trend observed in experimental studies [145].

## B.4.2 Impact of a random alloy fluctuations on the carrier transport in (In,Ga)N/GaN MQWs

In the last step, we move away from the VCA description of the system and include, in addition to quantum corrections, also random alloy fluctuations in the calculations. Figure B.4 (b) (blue) shows that, and this time in line with the experimental results by Galler et. al [145], the deep QW only contributes significantly to the radiative recombination when it is closest to the *p*-side (position 3). In fact, when including random alloy fluctuations in the calculations, the well closest to the *p*-side always has the largest contribution to total radiative recombination, as can be seen in Fig. B.5 (c).

To understand this behavior, Fig. B.6 (iii) depicts the electron and hole densities in the different wells as a function of the position of the deep well in the MQW systems. Looking at the electron density first, in comparison to the VCA calculations both including and excluding quantum corrections, random alloy fluctuations lead a reduction in electron density at the *n*-side QW. As discussed above and previously, quantum corrections can lead to increased electron transport, and including alloy fluctuations adds further to this effect due to the softening of the
barrier at the well interfaces [221]. As a consequence, the electrons can more easily ‘overshoot’
the wells in the MQWs, which can also be seen in the increased electron density beyond the *p*-side QW when alloy fluctuations and quantum corrections are included. However, in comparison
to the VCA result including quantum corrections, the electron density in the *p*-side well is only
slightly affected by alloy fluctuations.

In contrast, hole densities in the *p*-side QW are more dramatically changed by alloy fluctua-
tions. As Figs. B.6 (ii) and B.6 (iii) show, in comparison to the VCA description, alloy fluctuations
lead to an increase in the hole density in the *p*-side QW (position 3) even when the deep QW
is closest to the *n*-side (position 1) or in the centre (position 2) of the MQW system. While the
smoothing of the well barrier interface can increase hole transport, as in the case of electrons,
there are now also alloy disorder induced localization effects to contend with. As discussed in
our previous work, these localization effects are *detrimental* to hole transport [222] and result
in an increased hole density in the *p*-side QW. As a consequence, the well closest to the *p*-side
dominates radiative recombination

We note that there is still a reasonable hole density present in the *n*-side QW (Fig. B.6 (iii) (a-c)). In general, the distribution of carriers will also depend on the GaN barrier width and a 5 nm barrier is narrow enough to allow for some hole transport across the MQW [199]; a similar dependence of hole transmission on the barrier width has been seen in previous non-equilibrium Green's function studies [220]. Thus we expect that increasing the barrier width will mainly lead to a reduction of the hole density in the well furthest away from the *p*-side, but should to a lesser extent affect the hole density in the well closest to the *p*-side. Therefore, even for a larger barrier width than the here considered 5 nm, we expect that the recombination will still be dominated by the *p*-side QW.

We note that based on the VCA results we did not consider the EBL in the atomistic calcu-
lations. In general the EBL needs to be treated with an atomistic resolution. Previous studies
of (Al,Ga)N barriers in uni-polar device settings have found that the impact of these barriers is
lower than what is expected from a 1-D simulation for both electrons [45] and holes [231]. Thus
given that our VCA calculations show that the presence of the EBL is of secondary importance
for our study, we expect a similar conclusion when treating the EBL with alloy fluctuations.
Therefore, it is unlikely that the EBL impacts the here presented result, however this question
may be targeted in future studies.

## B.5 Conclusions

In this work we apply a 3-D quantum corrected multiscale simulation framework to gain
insight into the impact of random alloy fluctuations on the electron and hole distribution across
the active region of an (In,Ga)N/GaN LED. To study the spatial distribution of carriers we have
followed literature experimental studies [145] and analyzed the radiative recombination ratio in
a multi-quantum well system, where one of the wells in the system has a higher indium content
(deeper well) and its position is varied within the stack.

The here considered MQW systems are not only of interest for a comparison with experiment,
they provide also the ideal opportunity to benchmark and validate results from our in-house de-
veloped 3-D multiscale simulation framework against commercially available software packages.
To do so we treat the QWs in a virtual crystal approximation (VCA), to mimic the 1-D simulation
widely used in the literature for (In,Ga)N QWs and implemented in the commercial software
package nextnano. In addition, this study allows us also to compare the different schemes to
account for quantum corrections (localization landscape theory vs. Schrödinger-Poisson-DD
simulations) in the simulations. Overall, this analysis showed very good agreement between
results obtained from our in-house software and nextnano (without random alloy fluctuation).

Equipped with this benchmarked model, our analysis reveals that including (random) alloy
fluctuations in the calculations is vital for reproducing trends seen in experiment. More specifi-
cally, when using the widely employed virtual crystal approximation (VCA), the hole density
in the well closest to the *p*-doped region of the device is significantly reduced compared to our
atomistic random alloy calculation. As a consequence, and in contrast to the experiment, in VCA
the well closest to the *p*-side contributes very little to the radiative recombination process, an
effect that can be reduced by accounting for quantum corrections. While this leads to enhanced
radiative recombination from the well closest to the *p*-side, at least when this well is the deep
well, it still does not reflect the trends observed in the experimental studies. However, when
including random alloy fluctuations and quantum corrections in our 3-D simulation framework,
these effects lead to an increase in the hole density in the well closest to the *p*-side. Consequently,
this well dominates the radiative recombination process in line with the experimental data. We
note that in addition to quantum corrections and alloy fluctuations no further ingredients are
required (e.g. multi-population model) to explain the experimentally observed trends. Therefore,
our calculations highlight that alloy fluctuations are a key ingredient in simulations guiding the
design of III-N based devices. Thus, the here developed model presents an ideal starting point
for future calculations.
<figure><img src="image_54.png" /><figcaption>Figure B.6 – Hole density (black, solid), electron density (black, dashed), and radiative recombination rate (red, solid) averaged over each atomic plane along the transport direction. Results from calculations building on (i) a virtual crystal approximation (top), (ii) a virtual crystal including quantum corrections via localization landscape theory (LLT) (center) and a (iii) random alloy description including LLT-based quantum corrections (bottom); the deep well is located at (a) the *n*-side (left), (b) the center (middle) and (c) the *p*-side (right). The data are shown on a log scale.</figcaption></figure>
# Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities

## Outline of the current chapter

<table><tr><td>C.1 Introduction</td><td>190</td></tr><tr><td>C.2 Drift-diffusion equations and diffusion enhancement</td><td>191</td></tr><tr><td>C.3 Finite volume space discretization</td><td>192</td></tr><tr><td>C.3.1 Discrete thermodynamic consistency</td><td>193</td></tr><tr><td>C.3.2 The Scharfetter-Gummel scheme</td><td>193</td></tr><tr><td>C.3.3 SEDAN scheme</td><td>193</td></tr><tr><td>C.4 Simulations</td><td>194</td></tr><tr><td>C.5 Conclusion</td><td>197</td></tr></table>

This Appendix is a submitted [126] work, in collaboration with Patricio Farrell, Michael O'Donovan, Stefan Schulz and Thomas Koprucki.

We show the importance of using a thermodynamically consistent flux discretization when describing drift-diffusion processes within light emitting diode simulations. Using the classical Scharfetter-Gummel scheme with Fermi-Dirac statistics is an example of such an inconsistent scheme. In this case, for an (In,Ga)N multi quantum well device, the Fermi levels show an unphysical hump within the quantum well regions. This result originates from neglecting diffusion enhancement associated with Fermi-Dirac statistics in the numerical flux approximation. For a thermodynamically consistent scheme, such as the SEDAN scheme, the humps in the Fermi levels disappear. We show that thermodynamic inconsistency has
far reaching implications on the current-voltage curves and recombination rates.

## C.1 Introduction

In recent years drift-diffusion simulations have offered a numerically attainable method for studying carrier dynamics in a wide variety of devices including light-emitting diodes (LEDs) [237, 225, 219], transistors [251, 91] and solar cells [232, 255]. The physical interpretation of the model is quite straightforward: carriers in a device tend to diffuse from regions of high carrier density to low, and have drift motion due to an applied force such as an electric field in the device. On the other hand, the numerical implementation of the model can carry pitfalls which lead to an incorrect description of the device behaviour as we will highlight below in detail.

The purpose of this publication is therefore to show the practical unphysical implications of disobeying thermodynamic consistency. In the stationary case, thermodynamic consistency for discretized drift-diffusion equations can be defined by the demand that the zero bias solution coincides with the thermodynamic equilibrium. In the transient case, it is closely related to the fact that, for boundary conditions compatible with the thermodynamic equilibrium, the solution converges to this equilibrium when time tends to infinity. It is already known that disobeying this property causes non-physical dissipation in the steady state, see [26]. However, the inconsistent discrete approximation of the numerical fluxes has also more direct consequences for the quasi Fermi potentials and is thus important for accurately describing the physics of a device in the frame of drift-diffusion simulations. If one inconsistently approximates the fluxes the quasi Fermi level will show a completely wrong behavior in for instance quantum well regions, which are at the heart of modern LEDs. This has a knock-on effect for the description of the carrier densities and thus also recombination and current-voltage (IV) curves, which we will illustrate with an (In,Ga)N quantum well structure, a material system of strong interest for energy efficient solid state lighting [158] and for which considerable effort has been undertaken to develop advanced carrier transport models [221].

For the Boltzmann distribution, the classical Scharfetter-Gummel scheme [239] presents such a thermodynamically consistent scheme. Strictly monotonically increasing non-Boltzmann distribution functions lead to *diffusion enhancement*. Various extensions of the Scharfetter-Gummel scheme have been suggested to account for this effect, see [228, 171, 249]. Unfortunately, these schemes are not thermodynamically consistent. In [184, 185] a thermodynamically consistent generalization for Blakemore statistics (which is itself a special case of [120]) is presented in the spirit of [239] by solving local Dirichlet problems. But this generalization requires solving local nonlinear equations during assembly and the iterative solution of the coupled system. It is therefore computationally prohibitively expensive. A computationally more affordable approach is presented in [67]. On the other hand, in [26] the author presents another extension of the Scharfetter-Gummel scheme using a proper average of the nonlinear diffusion guaranteeing thermodynamic consistency for a specific choice of the distribution function. An alternative interpretation of this approach based on averaging the diffusion enhancement for a very general class of statistical distribution functions was given in [186]. Finally, the SEDAN scheme [265, 50, 2] includes the nonlinearity in the drift instead of the diffusion part of the flux and thus also yields a thermodynamically consistent scheme.

The remainder of this paper is organized as follows: In Section C.2, we describe the bipolar drift-diffusion model for charge transport in semiconductors. Its finite volume discretization including the flux discretizations is described in Section C.3. The formal definition of discrete
thermodynamic consistency is also presented in this section. We compare thermodynamically consistent and inconsistent schemes in Section C.4 by studying the distribution of densities and quasi Fermi levels within an (In,Ga)N quantum well (QW) system, which is embedded in a p-i-n junction. Finally, we conclude in Section C.5.

## C.2 Drift-diffusion equations and diffusion enhancement

We briefly introduce a model based on nonlinear partial differential equations which describes bipolar charge transport in a semiconductor. More details can be found in [128]. The dependence of the carrier densities $n$ and $p$ on the chemical potentials for electrons and holes $\eta_n$ and $\eta_p$ are described by a statistical distribution function $\mathcal{F}$ as well as conduction and valence band densities of states $N_c$ and $N_v$ via the state equations $n = N_c\mathcal{F}(\eta_n)$ and $p = N_v\mathcal{F}(\eta_p)$. Typical choices for the distribution function are $\mathcal{F}(\eta) = \exp(\eta)$, the so-called Boltzmann approximation, or $\mathcal{F}(\eta) = F_{1/2}(\eta) = \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{E^{1/2}}{e^{E-\eta+1}} dE$, namely the Fermi-Dirac integral of order 1/2 describing degenerate semiconductors.

The chemical potentials are related to the quasi-Fermi potentials of electrons and holes $\varphi_n$ and $\varphi_p$ via

$$ \eta_n = (q(\psi - \varphi_n) - E_c)/(k_B T) \text{ and } \eta_p = (q(\varphi_p - \psi) + E_v)/(k_B T). $$

Here $q$ denotes the elementary charge, $\psi$ the electrostatic potential, $k_B$ the Boltzmann constant, $T$ the temperature and $E_c$ and $E_v$ the conduction and valence band-edge energies. This model assumes that charge carriers behave as if they are in a bulk material, described by a 3-D density of states. In a slowly varying potential this is a valid description, however in a quantum well system the abrupt interface requires a more advanced treatment. To account for the quantum mechanical nature of electrons and holes one option is to solve the Schrödinger equation for the confining potential energy formed by $V_{c,v} = E_{c,v} - q\psi$. This is a numerically demanding approach, as the Schrödinger equation is an eigenvalue problem which would need to be solved self-consistently coupled to Poisson and drift-diffusion equations. While such a calculation is feasible in 1D, extending this to 2D or 3D structures is numerically prohibitive. Therefore, in recent years significant efforts have been undertaken to establish methods that account for quantum corrections but are computationally cheaper [129, 225]. One of these approaches is based on the so-called localization landscape theory (LLT), which allows to extract a (non-local) effective potential. It has been shown that LLT provides a good approximation of the single particle ground states, not only in square wells but also triangular wells which are relevant for systems with a polarization field – such as (In,Ga)N QWs [74]. The resulting (effective) confining potentials, $E_{c,v}^{\text{eff}}$, exhibit band edges that are softened and approximate the finite extent of carrier wavefunctions<fnref n="1" />. In the framework of LLT a linear system of equations is solved to determine the effective potential without introducing extra free parameters [132]. Our implementation of this method in conjunction with a drift-diffusion based carrier transport solver is discussed in more detail in [219] and [221]. In [219] we have compared results of our quantum corrected drift-diffusion model employing LLT with the results of a commercially available software package utilizing a ‘standard’ Schrödinger-Poisson solver. Our findings indicate that LLT can produce results in good agreement with the fully coupled Schrödinger-Poisson-drift-diffusion solver, highlighting that LLT captures quantum mechanical corrections sufficiently. In recent years this method has been used to study transport behaviour of numerous semi-conductor structures include LEDs [202, 204, 206, 211], blocking layers [231] and superlattices [256]; moreover, these

<footnote n="1">1. From here on within this work $E_{c,v}$ shall refer to the band edge values which are modified to account for the effective potential, $E_{c,v} = E_{c,v}^{\text{eff}}$.</footnote>
methods have been used successfully alongside experimental studies to gain insight into device behaviour [203, 236].

We model a bipolar semiconductor device as a domain $\Omega \subset \mathbb{R}^d$ where the carrier transport in a self-consistent electrical field is described by a system of partial differential equations. In the steady-state case this drift-diffusion system consists of Poisson's equation for $\psi$ and continuity equations for electrons and holes:

$$ -\nabla \cdot (\varepsilon_0 \varepsilon_r \nabla \psi) = q(C + p - n), \quad \vec{x} \in \Omega, \qquad (C.1) $$

$$ -\nabla \cdot \mathbf{j}_n = -qR, \quad \nabla \cdot \mathbf{j}_p = -qR, \quad \vec{x} \in \Omega. \qquad (C.2) $$

Here, $\varepsilon_r$ is the relative permittivity, $C$ is the net doping profile, and $R = R(n,p)$ describes carrier recombination. Electron and hole current densities can be expressed in terms of quasi-Fermi potentials by

$$ \mathbf{j}_n = -q\mu_n n \nabla \varphi_n, \quad \mathbf{j}_p = -q\mu_p p \nabla \varphi_p, \qquad (C.3) $$

or for any strictly monotonic Fermi-like distribution function $\mathcal{F}(\eta)$ in drift-diffusion form

$$ \mathbf{j}_n = q\mu_n \left[ U_T g\left(\frac{n}{N_c}\right) \nabla n - n \nabla \left(\psi - \frac{E_c}{q}\right) \right] \\ \text{and } \mathbf{j}_p = -q\mu_p \left[ U_T g\left(\frac{p}{N_v}\right) \nabla p + p \nabla \left(\psi + \frac{E_v}{q}\right) \right], \qquad (C.4) $$

where $\mu_n$ and $\mu_p$ denote the electron and hole mobilities, respectively, and $U_T = k_B T/q$ is the thermal voltage. The factor $g$ can be defined in terms of densities, $g(x) = x(\mathcal{F}^{-1})'(x)$, for $x \in \mathbb{R}$. This factor is the so-called diffusion enhancement appearing as a density-dependent modification factor in the generalized Einstein relation, see [257], leading in general to a non-linear diffusion coefficient. For the Boltzmann distribution, $\mathcal{F}(\eta) = \exp(\eta)$, we have $g \equiv 1$ and the current expressions (C.4) reduce to the usual ones with linear diffusion.

### C.3 Finite volume space discretization

We discretize the domain $\Omega$ using the Voronoï box based finite volume method introduced in [204], also known as “box method” due to [15]. It uses a simplicial boundary conforming Delaunay grid ([248]) which allows to obtain control volumes surrounding each given collocation point $\mathbf{x}_K$ by joining the circumcenters of the simplices containing it, see [128] for details.

Let $\partial K$ denote the boundary of the control volume $K$, and $|\xi|$ the measure of a geometrical object $\xi$. For each control volume $K$, we integrate the continuity equation (C.2) and apply Gauss's theorem to the integral of the flux divergence. Restricting our considerations to the electron transport equation, we obtain

$$ \begin{aligned} 0 &= \int_{\partial K} \mathbf{j}_n \cdot \mathbf{n} \, ds - \int_K qR \, dx = \sum_{L \text{ neighbor of } K} \int_{\partial K \cap \partial L} \mathbf{j}_n \cdot \mathbf{n}_{KL} \, ds - \int_K qR \, dx \\ &\approx \sum_{L \text{ neighbor of } K} |\partial K \cap \partial L| \mathbf{j}_{n,KL} - q|K|R(n_K, p_K), \end{aligned} \qquad (C.5) $$

where $\mathbf{n}$ is the internal unit normal to $\partial K$ and $\mathbf{n}_{KL}$ is the internal unit normal to the interface
$\partial K \cap \partial L$ for each neighbor $L$ of $K$. The values $n_K, p_K$ are the numerical approximations of the densities $n$ and $p$ at the collocation points $x_K$, and $j_{n,KL}$ are approximations of the normal currents through $\partial K \cap \partial L$. In the same manner the discretization of the Poisson equation can be obtained. A more detailed discussion of this method can be found in [128].

### C.3.1 Discrete thermodynamic consistency

One property which holds on a continuous level to avoid unphysical state dissipation is the preservation of thermodynamic equilibrium [128]. Mathematically, this means that vanishing fluxes shall imply constant quasi Fermi potentials. The classical discrete counterpart of this property is formulated as below (see for example [186, 128]): a numerical flux $j = j_{KL}$ is said to be thermodynamically consistent if it satisfies an analogous discrete relation, i.e.

$$j = 0 \text{ implies } \delta\varphi_{KL} = 0, \qquad (C.6)$$

where $\delta\varphi_{KL} = (\varphi_L - \varphi_K)/U_T$. Similarly, we define $\delta\eta_{KL} = \eta_L - \eta_K$ and $\delta\psi_{KL} = (\psi_L - \psi_K)/U_T$ and $\delta E_{KL} = (E_{c,L} - E_{c,K})/(qU_T)$. We point out that the condition (C.6) holds in equilibrium. Here, we introduce a stronger notion of thermodynamic consistency, which holds outside of equilibrium, namely

$$j \le 0 \text{ implies } \delta\varphi_{KL} \ge 0 \quad \text{and} \quad j \ge 0 \text{ implies } \delta\varphi_{KL} \le 0. \qquad (C.7)$$

An important property of defining thermodynamical consistency like above is that the sign of the numerical current is consistent with that of its continuous counterpart (C.3). Thermodynamic consistency is also important, when coupling the van Roosbroeck system to heat transport models [128]. We discuss now different numerical fluxes that may be used within a Voronoï finite volume framework.

### C.3.2 The Scharfetter-Gummel scheme

First, we introduce the well known, classical Scharfetter-Gummel flux approximation [239] given by

$$j_{sg} = j_0 \{B(\delta\psi_{KL} - \delta E_{KL})\mathcal{F}(\eta_L) - B(-\delta\psi_{KL} + \delta E_{KL})\mathcal{F}(\eta_K)\}, \qquad (C.8)$$

where the constant $j_0$ is given by $j_0 = q\mu_n N_c \frac{U_T}{h_{KL}}$ for $h_{KL} = |\vec{x}_K - \vec{x}_L|$, and $B$ is the Bernoulli function, $B(x) = \frac{x}{\exp(x)-1}$. It is important to point out that Scharfetter and Gummel introduced this numerical flux only in the Boltzmann regime, i.e. $\mathcal{F} = \exp$. In this case, the flux is thermodynamically consistent in the sense of (C.6). However, once we leave the Boltzmann regime, i.e. $\mathcal{F} \neq \exp$, and continue using (C.8) this numerical flux will no longer be thermodynamically consistent.

### C.3.3 SEDAN scheme

Next, we present the SEDAN scheme, which yields a thermodynamically consistent approach also for state equations which do not necessarily rely on the Boltzmann approximation. The earliest reference we could find for such an excess chemical potential scheme is the source code of the SEDAN III simulator [265], which explains the reason we use this name. A numerical analysis focused comparison of this flux approximation is given in [50] and simulation results are presented in [2]. The scheme is motivated by rearranging the drift part to include the excess
<figure><img src="image_55.png" /><figcaption>Figure C.1 – Conduction band edge (black) and quasi Fermi energy (red) at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac (FD) statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.</figcaption></figure>

chemical potential, $v^{ex} = \ln \mathcal{F}(\eta) - \eta$, yielding

$$j_{\text{sedan}} = j_0 \{B(Q_{\text{KL}}) \mathcal{F}(\eta_L) - B(-Q_{\text{KL}}) \mathcal{F}(\eta_K)\} \quad (C.9)$$

with

$$Q_{KL} = \delta\psi_{KL} - \delta E_{KL} + v_L^{ex} - v_K^{ex} = \delta\psi_{KL} - \delta E_{KL} - \delta\eta_{KL} + \ln \frac{\mathcal{F}(\eta_L)}{\mathcal{F}(\eta_K)}. \quad (C.10)$$

Note that, using the definition of $Q_{KL}$ and the fact that $e^x B(x) = B(-x)$, one can reformulate the SEDAN flux as

$$j_{\text{sedan}} = j_0 B(Q_{\text{KL}}) \mathcal{F}(\eta_L) (1 - e^{\delta \phi_{\text{KL}}}).$$

Therefore, it is easy to see that the SEDAN flux satisfies (C.7), since both $B(Q_{KL})$ and $\mathcal{F}(\eta_L)$ are positive: it is a thermodynamically consistent numerical flux.

Note that when applying Boltzmann statistics $v^{ex} = 0$ and the SEDAN flux becomes equivalent to the Scharfetter-Gummel expression. Therefore, in the next section, when displaying results using Boltzmann statistics, we only show results from one numerical scheme.

## C.4 Simulations

To illustrate the importance of using thermodynamically consistent flux approximations, we study a simple (In,Ga)N multi QW (MQW) system. In particular, we consider three QWs and the same set of parameters as in [219]. For large negative values of $\eta$ ($\eta \le -2$, which correspond roughly to densities below 14% of the effective density of states $N_c$, thus a low carrier density regime in the conduction band of the wells) the Boltzmann approximation provides a good estimate of the Fermi-Dirac statistics. Therefore in certain cases Scharfetter-Gummel scheme can offer a good description of the drift-diffusion model (e.g. Ref. [222]).

On the other hand there are situations where Boltzmann statistics will not suffice and the importance of using a thermodynamically consistent scheme becomes apparent. Figure C.1 (a) displays the conduction band edge and quasi Fermi energies of the MQW system when treated using Boltzmann statistics at a bias of 3.3V. Applying the Scharfetter-Gummel scheme (C.8) to Fermi-Dirac statistics for a bias of 3.3V leads to humps in the quasi Fermi energy within each
QW, see Figure C.1 (b).

Recalling the condition (C.7), for thermodynamically consistent schemes, the discrete gradient of the discrete quasi Fermi level should indicate the direction of the discrete electron flow (as it is in the continuous case, according to (C.3)). However in Figure C.1 (b), one can see that, when the Scharfetter-Gummel scheme is applied to Fermi-Dirac statistics, resulting in a thermodynamically *inconsistent* scheme, the direction of electron flow is to the right outside the QW regions (e.g. between -5 nm and 3 nm) but to the left inside the QW regions (e.g. between 10 nm and 15 nm), which is not in accordance with (C.7). Moreover, as there is no generation of carriers in the system this change in direction of the electron flux is highly unphysical.

<figure><img src="image_56.png" /><figcaption>Figure C.2 – Numerical electron flux averaged over each atomic plane at a bias of 3.3V shown for Boltzmann statistics (black), Fermi-Dirac statistics incorrectly using the Scharfetter-Gummel (SG) flux discretization (red) and Fermi-Dirac statistics correctly using the SEDAN flux discretization (blue).</figcaption></figure>

These humps in the quasi Fermi level are not present if one uses the Scharfetter-Gummel scheme with the Boltzmann approximation, which is a thermodynamically *consistent* scheme (Figure C.1 (a)). Similarly, using Fermi-Dirac statistics with the SEDAN scheme does not exhibit the unphysical humps in the Fermi level (Figure C.1 (c)).

Another perspective to interpret thermodynamic inconsistency is to note the incorrect interplay between quasi Fermi energies and local current fluxes. We see in Figure C.2 that the local numerical electron fluxes are positive and decrease monotonically across the QW regions. This is true for all three settings, that is, Boltzmann statistics in combination with the classical Scharfetter-Gummel scheme, as well as Fermi-Dirac statistics using both consistent and inconsistent numerical schemes. By our definition of a thermodynamically consistent scheme (C.7), a positive local numerical flux should imply a negative quasi Fermi potential discrete gradient. In fact, this is true for both consistent settings as one can see in Figure C.1 (a) and (c). However, for the inconsistent case the derivative of the quasi Fermi energies can become positive inside the wells, Fig. C.1 (b), also in disagreement of how the flux should behave at the continuous level (see the relation between the flux and the quasi Fermi potential (C.3)).

Previous studies of the numerical flux approximations have shown that a thermodynamically inconsistent scheme can result in the incorrect sign of the particle flux [186]. This is also reflected by the fact that (C.6) holds only for consistent schemes, which guarantee the physically correct sign of the current also far from equilibrium.
<figure><img src="image_57.png" /><figcaption>Figure C.3 – Electron density at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.</figcaption></figure>

The physical reason of why an inconsistent scheme produces humps within the QWs becomes apparent when looking at the corresponding densities at a bias of 3.3V. The Fermi-Dirac function grows like a polynomial while the Boltzmann approximation grows for large densities exponentially, see for example Figure 50.9 in [128]. This different behaviour leads to nonlinear diffusion, the so-called diffusion enhancement, for non-Boltzmann statistics, see (C.4), or [186]. The Scharfetter-Gummel scheme neglects the diffusion enhancement assuming only linear diffusion, which has a knock-on impact on the carrier density: the densities calculated using Boltzmann statistics (Fig. C.3 (a)) and using Fermi-Dirac statistics with the Scharfetter-Gummel scheme (Fig. C.3 (b)) are visibly indistinguishable. However, in order to produce the same density between a Boltzmann and Fermi-Dirac calculation, the quasi Fermi levels must differ. This results in the unusual behaviour of the quasi Fermi level exhibited in Figure C.1 (b). Comparing these densities with the correctly calculated density using Fermi-Dirac statistics in combination with the thermodynamically consistent SEDAN scheme (Figure C.3 (c)) we see that the choice of statistics function will impact carrier density in the well, and more strongly in the barrier regions at the here chosen example bias of 3.3V.

Because it influences the carrier density, thermodynamic inconsistency has direct implications for the computed recombination rates as well as the current-voltage (IV) curves. Next, we compare the recombination rates calculated using Boltzmann statistics with those calculated while incorrectly using the Scharfetter-Gummel scheme implementing Fermi-Dirac statistics.

This is highlighted in Figure C.4 (a), where the differences between Fermi-Dirac and Boltzmann-like behaviour are shown for the three recombination rates, calculated via

$$ \Delta \log(\text{recomb.}) = \log(R_i^{\text{SEDAN}}) - \log(R_i^{\text{Boltzmann}}), $$

where $R_i^s$ is the recombination rate associated with the process $i$ ($i \in \{\text{Shockley-Read-Hall, Radiative, Auger}\}$) calculated with the scheme $s$ ($s \in \{\text{Boltzmann, SEDAN}\}$). From this figure it becomes clear that the Boltzmann behaviour overall underestimates the recombination across the multi QW region of the device. In particular, the Auger recombination is underestimated by up to two orders of magnitude at a bias of 3.3V. These differences increase as the bias is increased (not shown). This can have consequences for overall device behaviour such as the internal quantum efficiency and the IV curves. The latter are shown in Figure C.4 (b), where the decreased recombination current displayed in the Boltzmann and thermodynamically inconsistent Fermi-Dirac scheme leads to an underestimate of the current density by close to an
<figure><img src="image_58.png" /></figure>

Figure C.4 – (a) Difference in magnitude of the Shockley-Read-Hall (SRH, black), radiative (red) and Auger (blue) recombination rates between Fermi-Dirac and Boltzmann statistics at a bias of 3.3V, calculated as described in the main text. (b) Current density-voltage (IV) curves using Boltzmann statistics (black, dashed), Fermi-Dirac statistics using the Scharfetter-Gummel scheme (SG, red) and Fermi-Dirac statistics using the SEDAN scheme (black, solid), shown on a log scale.

order of magnitude at 3.6V.

The results highlighted above indicate that Fermi statistics implemented using a thermodynamically inconsistent scheme will result in Boltzmann-like behaviour in LED simulations – at least in terms of carrier and current densities. If this is extended to laser simulations the consequences can be even more dramatic, as the gain calculation depends on the difference between the electron and hole quasi Fermi energies [14], expressed by the so-called Fermi voltage. In this case the unphysical humps seen in Figure C.1 (b) will lead to an incorrect prediction of the transparency density.

## C.5 Conclusion

In this paper, we have shown the importance of using a thermodynamically consistent flux discretization when describing drift-diffusion processes within quantum well devices.

Using the classical Scharfetter-Gummel scheme with Fermi-Dirac statistics is an example of such an inconsistent scheme. Here we studied an (In,Ga)N multi quantum well structure as an example since it is a very important material system for optoelectronic devices. In this case, the Fermi levels show humps within the quantum wells resulting in an unphysical description of the direction of the current, e.g. assuming the usual continuous expression. This is explained by the omission of diffusion enhancement from the numerical current expression, that leads to a similar density distribution as using Boltzmann statistics. This has a knock-on effect for recombination and current-voltage behaviour, where using Fermi-Dirac statistics with a thermodynamically inconsistent scheme may incorrectly predict a Boltzmann-like behaviour.

Contrarily, for a thermodynamically consistent scheme, such as the SEDAN scheme, these unphysical humps in the Fermi levels disappear and accurate current curves and recombination
processes are predicted. Thus, thermodynamically consistent schemes are essential to address
open questions, such as the efficiency drop in modern light emitting devices and to reliably
guide their design.
# Theoretical study of the impact of alloy disorder on carrier transport and recombination processes in deep UV (Al,Ga)N light emitters

This work is an article published in Applied Physics Letters [133], in collaboration with Robert Finn, Michael O'Donovan, Patricio Farrell, Timo Streckenbach, Thomas Koprucki and Stefan Schulz.

Aluminium gallium nitride ((Al,Ga)N) has gained significant attention in recent years due to its potential for highly efficient light emitters operating in the deep ultra-violet (UV) range (< 280 nm). However, given that current devices exhibit extremely low efficiencies, understanding the fundamental properties of (Al,Ga)N-based systems is of key importance. Here, using a multi-scale simulation framework, we study the impact of alloy disorder on carrier transport, radiative and non-radiative recombination processes in a c-plane Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N quantum well embedded in a p-n junction. Our calculations reveal that alloy fluctuations can open "percolative" pathways that promote transport for the electrons and holes into the quantum well region. Such an effect is neglected in conventional, and widely used transport simulations. Moreover, we find that the resulting increased carrier density and alloy induced carrier localization effects significantly increase non-radiative Auger-Meitner recombination in comparison to the radiative process. Thus, to avoid such non-radiative process and potentially related material degradation, a careful design (wider well, multi quantum wells) of the active region is required to improve the efficiency of deep UV light emitters.

Light emitters operating in the ultraviolet (UV) spectral range have received significant attention in recent years due to their importance for a wide range of applications [182]. A region of particular interest is the deep-UV-C wavelength window (< 280 nm), which enables for instance water purification and sterilization. Aluminium gallium nitride ((Al,Ga)N) based
light emitting diodes (LEDs) are ideal candidates for such applications since their emission wavelengths can in principle be tuned across the entire UV spectrum and they do not require the use of toxic mercury. Unfortunately, the external quantum efficiency (EQE) of (Al,Ga)N-based LEDs is very low (≤1%) in the deep UV-C range and multiple factors contribute to this, e.g. low light extraction efficiencies and poor carrier injection, high defect densities and thus reduced radiative recombination in the active region [4, 182].

Understanding the fundamental material and device properties is of central importance in improving the efficiency of such devices. Theory and simulation plays an important role for guiding the device design. However, III-N systems and alloys, such as (In,Ga)N, (Al,In)N but also (Al,Ga)N, exhibit in general very different properties when compared to other III-V materials, e.g. (In,Ga)As. For instance, experimental and theoretical studies reveal strong alloy fluctuation induced carrier localization effects in III-N alloys and heterostructures [82, 243, 242, 103, 261, 134]. Accounting for such effects present a significant challenge for an accurate theoretical and numerical description, since a three-dimensional (3D) treatment of, for instance, a quantum well (QW), is required ideally on an atomistic level. Numerical challenges are further amplified when it comes to simulating LED structures, which require not only large 3D supercells but also self-consistent calculations. Significant progress has been made to address these questions [225, 103, 221]. But, in comparison to (In,Ga)N-based systems, understanding the fundamental impact of alloy fluctuations on electronic, optical and carrier transport properties in (Al,Ga)N-based QW systems is sparse.

Here, we target these questions for an Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N QW embedded in *p*- and *n*-doped regions, thus a system relevant for deep UV (Al,Ga)N-based LEDs. The calculations build on our recently established quantum corrected multi-scale simulation framework that connects atomistic electronic structure theory with a drift-diffusion scheme to analyze carrier transport. To target recombination processes, the model is coupled with the widely used ABC model [226], where empirical *A*, *B* and *C* coefficients from the literature are employed. Our results indicate that alloy fluctuations can promote carrier transport through the intrinsic barriers by opening up percolative paths that are lower in energy. These pathways appear to improve the carrier injection into the QW active region – a feature that is neglected in ‘conventional’ calculations utilizing a virtual crystal approximation (VCA). Moreover, while we find that alloy disorder and carrier localization effects enhance radiative recombination, we find also that non-radiative Auger-Meitner recombination processes are strongly increased. The resulting increase in high energy carriers via non-radiative Auger-Meitner recombination may also lead to a degradation of the material as indicated in experimental studies [148]. As a consequence, our results suggest that alloy disorder can have a detrimental effect on the device efficiency. To reduce Auger-Meitner related processes, increasing the well width or employing multi QW (MQW) active regions to reduce the carrier density may be a way to improve quantum efficiencies in UV emitters, as the above discussed percolation paths may be exploited to distribute the carriers more evenly between MQWs.

In order to investigate the impact alloy disorder has on carrier transport and recombination processes in (Al,Ga)N-based LEDs, we describe the electronic structure of the active region by means of an empirical nearest neighbour sp³ tight-binding (TB) model. We note that while we treat the TB matrix elements as free parameters to reproduce hybrid-functional DFT band structures, the TB matrix could also be constructed from maximally localized Wannier functions generated from DFT [259]. This would thus provide an avenue to extend our method towards a DFT-based TB Hamiltonian. In general, our employed TB framework treats strain and polarization field fluctuations due to random alloy fluctuations with an atomistic resolution. We have recently employed the TB model to study the influence of alloy fluctuations on the electronic and optical properties of (Al,Ga)N/AlN QWs. Our work revealed that random alloy fluctuations are sufficient
<figure><img src="image_59.png" /><figcaption>Figure D.1 – Schematic illustration of the (Al,Ga)N-based *p*-i-*n* system underlying the simulations. The intrinsic barriers plus quantum well region is denoted as the "atomistic" mesh in the main text, while the *n*- and *p*-regions are described by a sparser device mesh.</figcaption></figure>

to lead to carrier localization effects in such systems. A detailed description of the model and
electronic structure results is given in Ref. [134].

To connect our TB model to drift-diffusion (DD) simulations we have recently developed the following procedure [221]. As a first step a "local" TB Hamiltonian is constructed. By diagonalis- ing this local Hamiltonian an energy landscape that accounts for alloy induced fluctuations in the conduction and valence band profile is obtained. This energy landscape with an atomistic resolution is mapped onto a finite element mesh using the software library pdelib [141]. The underlying finite element mesh has as many nodes as the atoms in the system. To target LED device simulations, the "atomistic" mesh is connected to a sparser device mesh representing *p*- and *n*-doped regions. Using localization landscape theory (LLT) [225] in conjunction with the software package ddfermi, [104] quantum corrected DD simulations are performed. More details about the carrier transport calculations, the mesh generation and numerical aspects can be found in Refs. [221] and [128].

The above outlined quantum corrected multi-scale simulation framework is here applied to a *p*-i-*n* system where the intrinsic (*i*) region contains a 2.3 nm wide Al<sub>0.7</sub>Ga<sub>0.3</sub>N QW embedded into approximately 10 nm wide intrinsic Al<sub>0.8</sub>Ga<sub>0.2</sub>N barriers, following the UV-C device design in Ref. [235]. The atomistic mesh region describing the active region of the device has the dimensions of roughly 17.4 × 15.1 × 22.2 nm<sup>3</sup>, corresponding to ≈ 560,000 atoms, and is coupled with a sparser device mesh by adding 40 nm of *n*- and *p*-doped regions, respectively. The doping density is 5 × 10<sup>19</sup> cm<sup>-3</sup> following the set-up described in Ref. [244]. A schematic illustration of the structure is given in Fig. D.1. To connect the conduction band and valence band edges between the atomistic and the sparser mesh, we use a blend of Gaussian softening and LLT to avoid discontinuities at such an interface. Details on the Gaussian softening and the effect of LLT on the energy landscape can be found elsewhere [221]. To test the impact of the in-plane dimensions of the simulation cell and thus also the influence of the alloy microstructure on the results, the calculations have been repeated for different cell sizes (not shown). These different calculations revealed in general the same trends as reported below, so that the above detailed simulation set up gives a representative picture of the impact of alloy disorder on the carrier transport in the discussed device.

To investigate the impact alloy fluctuations have on the carrier transport and the recombina-
tion process we proceed as follows: Our starting point is a fully atomistic description, meaning
that alloy fluctuations in both the intrinsic (Al,Ga)N QW and (Al,Ga)N barriers are accounted
for (cf. Fig. D.1). To study the impact of alloy disorder on the carrier injection, we investigate
a second set-up in which alloy fluctuations in the (Al,Ga)N QW are treated on an atomistic
level, however the intrinsic (Al,Ga)N barriers are described in a VCA. Finally, we employ a
<figure><img src="image_60.png" /></figure>

Figure D.2 – (a) Current-voltage curves, on a log-scale, for an Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N quantum well embedded in a *p*-*n* junction. The data are shown for simulations that (i) account for random alloy (RA) fluctuations in the well and barriers (blue), (ii) account for RA fluctuations in the well but treat the barrier in a virtual crystal approximation (VCA) (red), and (iii) use a VCA in the entire *p*-*i*-*n* structure (pink). The inset shows the I-V curves on a linear scale.

VCA description of the entire QW plus barrier system. This last setting is basically equivalent to one-dimensional carrier transport simulations widely employed in commercially available software packages. In order to establish a comparable VCA description of the targeted system, without introducing additional free parameters (e.g. alloy dependence and or bowing parameters of material parameters), the atomistic energy landscape is averaged over the corresponding areas (e.g. well or barrier). To achieve this averaging procedure, the QW was treated initially with Dirichlet boundary conditions so that, except from alloy fluctuation induced variations, the conduction and valence band edges are flat. A VCA description of the barriers is then straightforward and achieved by averaging the atomistic energy landscape over the entire barrier volume. However, the QW VCA description is slightly more complicated as piezoelectric and spontaneous polarization induced built-in fields across the QW cause the band edges to slope. Thus, a volume average over the QW region is not possible. In this case, the atomistic energy landscape is averaged across each *c*-plane along the growth direction (*c*-axis). While this can lead to slight fluctuations in the energy landscape, the approach chosen here is better suited to compare the different calculations, since it avoids constructing a VCA that requires specific information on how to average e.g. material and TB parameters with composition, or to perform a separate VCA built-in field calculation.

Building on these different simulation set-ups, Fig. D.2 displays the calculated current-voltage (I-V) curves of the above described (Al,Ga)N *p*-*i*-*n* structure using (i) the fully atomistic approach (RA), (ii) accounting for alloy fluctuations in the well but not the barrier (RA QW/VCA Barriers) and (iii) a VCA of the entire system. One can infer from this figure that the fully atomistic treatment gives independent of the applied voltage the highest current density, while the full VCA description of the system (well and barriers) results in the lowest density. Figure D.2 also reveals that the current density, especially at lower bias values, is increased when accounting for alloy fluctuations in the barrier; this can be seen from the comparison of the I-V curves for the fully atomistic and ‘partial’ VCA (RA QW, VCA Barriers) treatment. Thus, our results indicate
<figure><img src="image_61.png" /></figure>

Figure D.3 – 3D isosurface plot of the current density in the intrinsic (Al,Ga)N quantum well and barrier regions at a bias of 5.5 V, using the fully atomistic description. The current density for the holes (blue) and electrons (red) are plotted at 700 Å cm⁻², reflecting locally high current densities due to alloy fluctuations; dashed black lines indicate well boundaries.

that alloy fluctuations are beneficial for the carrier injection into the QW region.

To visualise this aspect further, Fig. D.3 depicts an isosurface plot of the current density in the intrinsic (Al,Ga)N QW and barrier regions, cf. Fig. D.1, at a fixed bias point of 5.5 V obtained from the fully atomistic description of the system. One can infer that alloy fluctuation induced band edge energy fluctuations lead to regions of higher and lower current densities. Thus, the fluctuations in the band edge energies enable percolation paths into the active region; such pathways are absent in a VCA description of the system. We note that this percolation path effect has also be observed in (In,Ga)N QWs and (Al,Ga)N barriers [263, 45, 231, 264]. We note that both electron and hole transport are affected by these percolation paths. However, due to the lower effective electron mass, the electron charge density is more smeared out, and compared to the holes, fewer regions reach the isosurface value at which the current density is plotted in Fig. D.3. As a consequence of these percolation paths one expects higher carrier densities in the active region in comparison to a simulation that neglects alloy disorder. Also, as we will see below, the percolation paths can lead to an increased hole density in the intrinsic barrier on the *n*-side of the studied system. Such a situation may be beneficial for inter-well transport in (Al,Ga)N multi-QW systems as recently reported in (In,Ga)N-based LEDs [203].

To study the impact of alloy fluctuations on the carrier transport and recombination in (Al,Ga)N-based LED structures in further detail, Fig. D.4 depicts electron and hole carrier densities above the turn-on voltage of the device, at a bias of 5.5 V. Figure D.4 (a) displays the *averaged* carrier density, over each *c*-plane, for the random alloy system along the *c*-axis (growth direction) together with the VCA data which neglects alloy disorder. When comparing the *averaged* random alloy data with the VCA results, it becomes clear that (*on average*) both electron and hole carrier densities are increased in the well when taking alloy fluctuations into account. Figure D.4 (a) indicates that the *average* hole charge density does not leak as much into the intrinsic (Al,Ga)N barrier material on the *n*-side of system (region below 10.7 nm), when compared to the hole density profile in VCA; one may attribute this to alloy induced hole localization effects which are not accounted for by the VCA description. The *average* electron density extends further into the (Al,Ga)N barrier region on the *p*-side (region above 13 nm) when compared to the VCA data; we relate this to the effect that alloy disorder and quantum corrections effectively reduce the quantum confinement between well and barrier and that the electrons are less strongly affected by alloy induced carrier localization effects as already seen
<figure><img src="image_62.png" /></figure>

Figure D.4 – Electron (red) and hole (blue) density along the *c*-axis of a Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N quantum well embedded in a *p*-*i*-*n* junction at a bias of 5.5 V; data are shown for the intrinsic regions of the device. (a) Averaged carrier distribution of the fully atomistic calculation over each *c*-plane along the *c*-axis (solid lines), and the carrier distribution for the virtual crystal approximation (dashed lines); (b) Scatter plot of the carrier distribution in the fully atomistic calculation. The vertical dashed black lines indicate well boundaries. Note (b) is an order of magnitude larger than (a) on the *y*-axis.

in (In,Ga)N-based systems [221]. However, it must be noted that the random alloy data in Fig. D.4 (a) are averaged over each *c*-plane. Thus, while on average the hole density may be smaller in e.g. the intrinsic barrier of the *n*-region when compared to VCA, local alloy regions may have high hole densities. To shed light onto this behavior, Fig. D.4 (b) displays a scatter plot of the carrier densities along the *c*-axis in the intrinsic (Al,Ga)N regions. The symbols give the carrier density at each lattice site within the respective *c*-planes. Figure D.4 (b) reveals strong fluctuations in the *local* densities which may be an order of magnitude higher than the *average* densities displayed in Fig. D.4 (a). The *local* hole carrier density is also noticeably higher
<figure><img src="image_63.png" /></figure>

Figure D.5 – Radiative and Auger-Meitner non-radiative recombination rates along the *c*-axis of the Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N *p*-*i*-*n* system at a bias of 5.5 V; the data are shown for the intrinsic regions of the device. Solid lines: results from the atomistic calculations taking alloy fluctuations into account (averaged over each *c*-plane). Dashed lines: data obtained from a virtual crystal approximation of the same system; vertical dashed black lines indicate well boundaries.

when compared to the electrons. In general, due to the higher effective mass, carrier localization for holes is more pronounced when compared to electrons [134]. Therefore, we conclude that percolative pathways, can promote easier transport of the carriers through the barriers. Thus, holes may transport to the *n*-side of the system more easily when compared to predictions from a VCA model. This aspect is of interest when designing multi-QW-based (Al,Ga)N LEDs, since it can help with a more equal distribution of holes between the wells, and thus be important for recombination processes.

Equipped with this information we turn our attention now to the recombination process in the system. To do so we employ the widely used ABC model. In the following we mainly focus our attention on the radiative and Auger-Meitner (non-radiative) processes; with increased current densities radiative and non-radiative Auger-Meitner processes become dominant over Shockley-Read-Hall (SRH) recombination as these contributions scale with carrier density *n* as ∝ *n*² (radiative) and ∝ *n*³ (Auger-Meitner) compared to ∝ *n* (SRH). We note that defect related non-radiative processes, such as SRH and trap assisted tunnelling (TAT) processes can also affect the device performance of UV-C (Al,Ga)N-based LEDs significantly. While we take SRH into account, TAT processes have not been considered here but could be accounted for as a SRH-like process with different capture rates as detailed in Ref. [233]. Turning to radiative and Auger-Meitner (non-radiative) recombination rates, Fig. D.5 depicts these rates calculated (i) including random alloy fluctuations and thus carrier localization effects and (ii) using VCA for the entire system, consequently excluding alloy disorder. Figure D.5 displays the random alloy (RA) data averaged over each *c*-plane. We stress that the recombination rates are calculated using the *local* carrier density data, thus accounting for the effect that locally one may encounter both high electron and hole densities. The required radiative (*B*) and Auger-Meitner (*C*<sub>nnp</sub> and *C*<sub>npp</sub>) recombination coefficients are taken from Ref. [244]. We note that building on this often used carrier density dependent recombination model, the wave function character of electrons and holes is not taken into account. While the symmetry of the wave functions can affect radiative and non-radiative recombination rates, analyzing the spatial overlap of
electron and hole charge density provides first insight into the impact of alloy disorder on the recombination rates. Comparing the atomistic calculation with the VCA data it becomes clear that both radiative and non-radiative contributions are enhanced by alloy disorder since this can improve the carrier injection but also induce carrier localization in the wells. But, upon further inspection, Fig. D.5 reveals that the non-radiative rate increases more strongly when compared to the radiative. In the VCA the ratio of radiative to non-radiative peak values is $r^{\text{VCA}} \approx 19$ while in the atomistic simulation it is $r^{\text{RA}} \approx 2$ at a bias of 5.5 V. In general, and using the ABC model [226] here as a starting point, the radiative recombination contributes with $Bn^2$ while the Auger-Meitner non-radiative rate with $Cn^3$, where *n* is the carrier density. Thus, when alloy disorder increases the carrier density in local regions of the well, the Auger-Meitner rate is expected to increase more quickly than the radiative rate. Moreover, the Auger-Meitner rate is a three-particle process and has basically two contributions, namely an electron-electron-hole and a hole-hole-electron process [180]. Given that alloy induced carrier localization effects are significant for holes, locally the hole-hole wave function overlap is (strongly) increased. In turn this can result in increased Auger-Meitner rates [210]. Taking all this together, alloy disorder and related carrier localization effects are expected to lead to a significant increase of the non-radiative Auger-Meitner effect and thus have a detrimental effect on the performance of (deep) UV LEDs; a similar effect has been seen (In,Ga)N-base LEDs [264]. Moreover, strong Auger-Meitner non-radiative recombination may also lead to material degradation effects in deep-UV LEDs, as recently indicated in experiments [148, 233].

In summary, using a quantum corrected multi-scale simulation framework, we presented a detailed analysis of the impact of alloy fluctuations on carrier transport and recombination processes in an Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N single quantum well system embedded in a *p*-*n* junction. Understanding the fundamental properties of such systems provides useful insight into the physics of deep UV light emitters. Our results show that alloy disorder can lead to improved carrier injection into the active region, i.e. quantum well, through percolative pathways. Also, these pathways can give rise to higher hole densities in the intrinsic region near the *n*-side of a device, an aspect of potential interest for tailoring inter-well transport in (Al,Ga)N-based multi quantum well systems. In general, these percolation effects are not captured by standard virtual crystal approximations employed in many commercially available carrier transport simulation packages. Moreover, we find an increase in the carrier density in the well which is accompanied by alloy induced carrier localization effects. In comparison to virtual crystal simulations, both radiative and non-radiative Auger-Meitner recombination processes are increased. We attribute this to alloy induced carrier localization effects, resulting in high local carrier densities which affect the Auger-Meitner effect more strongly than the radiative process. Moreover, the high energy of such Auger-Meitner generated carriers may lead to a degradation of the material and thus an enhancement of Shockley-Read Hall recombination. Overall, our results indicate that alloy disorder opens up percolative pathways improving carrier transport, but the enhanced Auger effect can have a detrimental effect on the efficiency of the device. While we have focused here on relatively high Al contents in the well and barrier, the above drawn conclusions are also expected to hold for systems with lower Al contents, given that our recent atomistic studies reveal that even at only 10% Al in (Al,Ga)N hole localization effects are observed. Thus, in general a careful design of the active region of an UV LED, e.g. by using wider wells as discussed in recent literature [217, 262] or multi quantum wells by exploiting percolation pathways to improve the interwell transport, is necessary to reduce non-radiative recombination processes and prevent material damage.
Bibliography

[1] D. ABDEL, C. CHAINAIS-HILLAIRET, P. FARRELL, AND M. HERDA, *Numerical analysis of a finite volume scheme for charge transport in perovskite solar cells*, IMA Journal of Numerical Analysis, (2023). in press.

[2] D. ABDEL, P. FARRELL, AND J. FUHRMANN, *Assessing the quality of the excess chemical potential flux scheme for degenerate semiconductor device simulation*, Optical and Quantum Electronics, 53 (2021).

[3] T. AIKI AND A. MUNTEAN, *Existence and uniqueness of solutions to a mathematical model predicting service life of concrete structures*, Adv. Math. Sci. Appl., 19 (2009), pp. 109–129.

[4] H. AMANO, R. COLLAZO, C. DE SANTI, S. EINFELDT, M. FUNATO, J. GLAAB, S. HAGEDORN, A. HIRANO, H. HIRAYAMA, R. ISHII, Y. KASHIMA, Y. KAWAKAMI, R. KIRSTE, M. KNEISSL, R. MARTIN, F. MEHNKE, M. MENEGHINI, A. OUGAZZADEN, P. J. PARBROOK, S. RAJAN, P. REDDY, F. RÖMER, J. RUSCHEL, B. SARKAR, F. SCHOLZ, L. J. SCHOWALTER, P. SHIELDS, Z. SITAR, L. SULMONI, T. WANG, T. WERNICKE, M. WEYERS, B. WITZIGMANN, Y.-R. WU, T. WUNDERER, AND Y. ZHANG, *The 2020 UV emitter roadmap*, J. Phys. D: Appl. Phys, 53 (2020), p. 503001.

[5] D. ANDERSON AND J. DRONIOU, *An arbitrary-order scheme on generic meshes for miscible displacements in porous media*, SIAM J. Sci. Comput., 40 (2018), pp. b1020–b1054.

[6] B. ANDREIANOV, C. CANCÈS, AND A. MOUSSA, *A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs*, J. Funct. Anal., 273 (2017), pp. 3633–3670.

[7] P. F. ANTONIETTI, L. BEIRÃO DA VEIGA, AND G. MANZINI, eds., *The virtual element method and its applications*, vol. 31 of SEMA SIMAI Springer Ser., Cham: Birkhäuser, 2022.

[8] A. ARNOLD, J. A. CARRILLO, L. DESVILLETTES, J. DOLBEAULT, A. JÜNGEL, C. LEDERMAN, P. A. MARKOWICH, G. TOSCANI, AND C. VILLANI, *Entropies and equilibria of many-particle systems: An essay on recent research*, Monatsh. Math., 142 (2004), pp. 35–43.

[9] D. N. ARNOLD, *An interior penalty finite element method with discontinuous elements*, SIAM J. Numer. Anal., 19 (1982), pp. 742–760.

[10] D. N. ARNOLD, F. BREZZI, B. COCKBURN, AND L. D. MARINI, *Unified analysis of discontinuous Galerkin methods for elliptic problems*, SIAM J. Numer. Anal., 39 (2002), pp. 1749–1779.

[11] D. N. ARNOLD, G. DAVID, D. JERISON, S. MAYBORODA, AND M. FILOCHE, *Effective Confining Potential of Quantum States in Disordered Media*, Phys. Rev. Lett., 116 (2016), p. 056602.

[12] M. AUF DER MAUR, A. PECCHIA, G. PENAZZI, W. RODRIGUES, AND A. DI CARLO, *Efficiency Drop in Green InGaN/GaN Light Emitting Diodes: The Role of Random Alloy Fluctuations*, Phys. Rev. Lett., 116 (2016), p. 027401.

[13] I. BABUŠKA, *The finite element method with penalty*, Math. Comput., 27 (1973), pp. 221–228.
[14] U. BANDELOW, H. GAJЕWSKI, AND R. HÜNLICH, Fabry-perot lasers: *Thermodynamics-based modeling*, Optoelectronic Devices: Advanced Simulation and Analysis, (2005), pp. 63–85.

[15] R. BANK AND D. ROSE, *Some error estimates for the box method*, SIAM J Numer Anal, 24 (1987), pp. 777–787.

[16] C. BATAILLON, F. BOUCHON, C. CHAINAIS-HILLAIRET, C. DESGRANGES, E. HOARAU, F. MARTIN, S. PERRIN, M. TUPIN, AND J. TALANDIER, *Corrosion modelling of iron based alloy in nuclear waste repository*, Electrochimica Acta, 55 (2010), pp. 4451–4467.

[17] J. BEAR, *Dynamics of Fluids in Porous Media*, Courier Corporation, 1988.

[18] L. BEAUDE AND S. LEMAIRE, *ParaSkel++: a C++ platform for the high-performance, arbitrary-order, 2/3D numerical approximation of PDEs on general polytopal meshes using skeletal Galerkin methods*, Aug. 2021.

[19] L. BEIRÃO DA VEIGA, J. DRONIOU, AND G. MANZINI, *A unified approach for handling convection terms in finite volumes and mimetic discretization methods for elliptic problems*, IMA J. Numer. Anal., 31 (2011), pp. 1357–1401.

[20] L. BEIRÃO DA VEIGA, F. BREZZI, A. CANGIANI, G. MANZINI, L. D. MARINI, AND A. RUSSO, *Basic principles of virtual element methods*, Math. Models Methods Appl. Sci., 23 (2013), pp. 199–214.

[21] L. BEIRÃO DA VEIGA, F. BREZZI, L. D. MARINI, AND A. RUSSO, *The Hitchhiker's guide to the virtual element method*, Math. Models Methods Appl. Sci., 24 (2014), pp. 1541–1573.

[22] L. BEIRÃO DA VEIGA, F. DASSI, AND A. RUSSO, *High-order virtual element method on polyhedral meshes*, Comput. Math. Appl., 74 (2017), pp. 1110–1122.

[23] N. BEN ABDALLAH AND R. EL HAJJ, *Diffusion and guiding center approximation for particle transport in strong magnetic fields*, Kinet. Relat. Models, 1 (2008), pp. 331–354.

[24] I. BEN GHARBIA, C. CANCÈS, T. FANEY, M. JONVAL, AND Q. H. TRAN, *Robust resolution of single-phase chemical equilibrium using parametrization and Cartesian representation techniques*. working paper or preprint, 2023.

[25] M. BENDAHMANE, Z. KHALIL, AND M. SAAD, *Convergence of a finite volume scheme for gas-water flow in a multi-dimensional porous medium*, Math. Models Methods Appl. Sci., 24 (2014), pp. 145–185.

[26] M. BESSEMOULIN-CHATARD, *A finite volume scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme*, Numer. Math., 121 (2012), pp. 637–670.

[27] M. BESSEMOULIN-CHATARD AND C. CHAINAIS-HILLAIRET, *Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems*, J. Numer. Math., 25 (2017), pp. 147–168.

[28] ——, *Uniform-in-time bounds for approximate solutions of the drift-diffusion system*, Numer. Math., 141 (2019), pp. 881–916.

[29] M. BESSEMOULIN-CHATARD, C. CHAINAIS-HILLAIRET, AND M.-H. VIGNAL, *Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit*, SIAM J. Numer. Anal., 52 (2014), pp. 1666–1691.

[30] S. BIRNER, T. ZIBOLD, T. ANDLAUER, T. KUBIS, M. SABATHIL, A. TRELLAKIS, AND P. VOGL, *nextnano: General Purpose 3-D Simulations*, IEEE Transactions on Electron Devices, 54 (2007), pp. 2137–2142.

[31] J. BLAKEMORE, *Approximations for Fermi-Dirac integrals, especially the function $F_{1/2}(\eta)$ used to describe electron density in a semiconductor*, Solid-State Electronics, 25 (1982), pp. 1067–1076.
[32] X. BLANC AND E. LABOURASSE, *A positive scheme for diffusion problems on deformed meshes*, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 96 (2016), pp. 660–680.

[33] T. BODINEAU, J. LEBOWITZ, C. MOUHOT, AND C. VILLANI, *Lyapunov functionals for boundary-driven nonlinear drift-diffusion equations*, Nonlinearity, 27 (2014), pp. 2111–2132.

[34] L. BOLTZMANN, *Weitere Studien über das Wärmegleichgewicht unter Gasmolekülen*, k. und k. Hof- und Staatsdr., 1872.

[35] J. BONELLE, D. A. DI PIETRO, AND A. ERN, *Low-order reconstruction operators on polyhedral meshes: Application to Compatible Discrete Operator schemes*, Computer Aided Geometric Design, 35–36 (2015), pp. 27–41.

[36] F. BONIZZONI, M. BRAUKHOFF, A. JÜNGEL, AND I. PERUGIA, *A structure-preserving discontinuous Galerkin scheme for the Fisher-KPP equation*, Numer. Math., 146 (2020), pp. 119–157.

[37] M. BRAUKHOFF, I. PERUGIA, AND P. STOCKER, *An entropy structure preserving space-time formulation for cross-diffusion systems: Analysis and galerkin discretization*, SIAM Journal on Numerical Analysis, 60 (2022), pp. 364–395.

[38] K. BRENNER AND C. CANCÈS, *Improving Newton's method performance by parametrization: the case of the Richards equation*, SIAM J. Numer. Anal., 55 (2017), pp. 1760–1785.

[39] F. BREZZI, K. LIPNIKOV, AND M. SHASHKOV, *Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes*, SIAM J. Numer. Anal., 43 (2005), pp. 1872–1896.

[40] F. BREZZI, K. LIPNIKOV, AND V. SIMONCINI, *A family of mimetic finite difference methods on polygonal and polyhedral meshes*, Math. Models Methods Appl. Sci., 15 (2005), pp. 1533–1551.

[41] F. BREZZI, L. D. MARINI, S. MICHELETTI, P. PIETRA, R. SACCO, AND S. WANG, *Discretization of semiconductor device problems*. I, in Handbook of numerical analysis. Vol XIII. Special volume: Numerical methods in electromagnetics., Amsterdam: Elsevier/North Holland, 2005, pp. 317–441.

[42] F. BREZZI, L. D. MARINI, AND P. PIETRA, *Numerical simulation of semiconductor devices*, Comput. Methods Appl. Mech. Eng., 75 (1989), pp. 493–514.

[43] ——, *Two-dimensional exponential fitting and applications to drift-diffusion models*, SIAM J. Numer. Anal., 26 (1989), pp. 1342–1355.

[44] D. BROWNE, B. MAZUMDER, Y.-R. WU, AND J. SPECK, *Electron transport in unipolar InGaN/GaN multiple quantum well structures grown by NH₃ molecular beam epitaxy*, J. Appl. Phys., 117 (2015), p. 185703.

[45] D. A. BROWNE, M. N. FIREMAN, B. MAZUMDER, L. Y. KURITZKY, Y.-R. WU, AND J. S. SPECK, *Vertical transport through AlGaN barriers in heterostructures grown by ammonia molecular beam epitaxy and metalorganic chemical vapor deposition*, Semiconductor Science and Technology, 32 (2017), p. 025010.

[46] C. BUET AND S. DELLACHERIE, *On the Chang and Cooper scheme applied to a linear Fokker-Planck equation*, Commun. Math. Sci., 8 (2010), pp. 1079–1090.

[47] M. BURGER, M. DI FRANCESCO, J.-F. PIETSCHMANN, AND B. SCHLAKE, *Nonlinear cross-diffusion with size exclusion*, SIAM J. Math. Anal., 42 (2010), pp. 2842–2871.

[48] J.-S. CAMIER AND F. HERMELINE, *A monotone nonlinear finite volume method for approximating diffusion operators on general meshes*, Int. J. Numer. Methods Eng., 107 (2016), pp. 496–519.

[49] C. CANCÈS, *Energy stable numerical methods for porous media flow type problems*, Oil & Gas Science and Technology - Rev. IFP Énergies nouvelles, 73 (2018).
[50] C. CANCEès, C. CHAINAIS-HILLAIRET, J. FUHRMANN, AND B. GAUDEUL, *A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model*, IMA J. Numer. Anal., **41** (2021), pp. 271–314.

[51] C. CANCEès, C. CHAINAIS-HILLAIRET, M. HERDA, AND S. KRELL, *Large time behavior of nonlinear finite volume schemes for convection-diffusion equations*, SIAM J. Numer. Anal., **58** (2020), pp. 2544–2571.

[52] C. CANCEès, C. CHAINAIS-HILLAIRET, AND S. KRELL, *Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations*, Comput. Methods Appl. Math., **18** (2018), pp. 407–432.

[53] C. CANCEès, C. CHAINAIS-HILLAIRET, B. MERLET, F. RAIMONDI, AND J. VENEL, *Mathematical analysis of a thermodynamically consistent reduced model for iron corrosion*, Z. Angew. Math. Phys., **74** (2023), p. 31. Id/No 96.

[54] C. CANCEès AND B. GAUDEUL, *A convergent entropy diminishing finite volume scheme for a cross-diffusion system*, SIAM J. Numer. Anal., **58** (2020), pp. 2684–2710.

[55] C. CANCEès AND C. GUICHARD, *Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations*, Math. Comp., **85** (2016), pp. 549–580.

[56] ——, *Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure*, Found. Comput. Math., **17** (2017), pp. 1525–1584.

[57] C. CANCEès, F. NABET, AND M. VOHRALÍK, *Convergence and a posteriori error analysis for energy-stable finite element approximations of degenerate parabolic equations*, Math. Comput., **90** (2021), pp. 517–563.

[58] X. CAO AND H. HUANG, *An adaptive conservative finite volume method for Poisson-Nernst-Planck equations on a moving mesh*, Commun. Comput. Phys., **26** (2019), pp. 389–412.

[59] M. A. CARO, S. SCHULZ, S. B. HEALY, AND E. P. O'ReILLY, *Built-in field control in alloyed c-plane III-N quantum dots and wells*, J. Appl. Phys., **109** (2011), p. 084110.

[60] M. A. CARO, S. SCHULZ, AND E. P. O'ReILLY, *Theory of local electric polarization and its relation to internal strain: impact on the polarization potential and electronic properties of group-III nitrides*, Phys. Rev. B, **88** (2013), p. 214103.

[61] J. A. CARRILLO, A. JÜNGEL, P. A. MARKOWICH, G. TOSCANI, AND A. UNTERREITER, *Entropy dissipation methods for degenerate parabolic problems and generalized Sobolev inequalities*, Monatsh. Math., **133** (2001), pp. 1–82.

[62] J. A. CARRILLO AND G. TOSCANI, *Exponential convergence toward equilibrium for homogeneous Fokker-Planck-type equations*, Math. Methods Appl. Sci., **21** (1998), pp. 1269–1286.

[63] ——, *Asymptotic $L^1$-decay of solutions of the porous medium equation to self-similarity*, Indiana Univ. Math. J., **49** (2000), pp. 113–142.

[64] C. CHAINAIS-HILLAIRET, *Discrete duality finite volume schemes for two-dimensional drift-diffusion and energy-transport models*, Int. J. Numer. Methods Fluids, **59** (2009), pp. 239–257.

[65] C. CHAINAIS-HILLAIRET AND J. DRONIOU, *Convergence analysis of a mixed finite volume scheme for an elliptic-parabolic system modeling miscible fluid flows in porous media*, SIAM J. Numer. Anal., **45** (2007), pp. 2228–2258.

[66] ——, *Finite-volume schemes for noncoercive elliptic problems with Neumann boundary conditions*, IMA J. Numer. Anal., **31** (2011), pp. 61–85.
[67] C. CHAINAIS-HILLAIRET, R. EYMARD, AND J. FUHRMANN, *A monotone numerical flux for quasilinear convection-diffusion equation*, arXiv preprint arXiv:2209.11657, (2022).

[68] C. CHAINAIS-HILLAIRET AND F. FILBET, *Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model*, IMA J. Numer. Anal., 27 (2007), pp. 689–716.

[69] C. CHAINAIS-HILLAIRET AND M. HERDA, *Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations*, IMA J. Numer. Anal., 40 (2020), pp. 2473–2504.

[70] C. CHAINAIS-HILLAIRET, M. HERDA, S. LEMAIRE, AND J. MOATTI, *Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches*, Numer. Math., 151 (2022), pp. 963–1016.

[71] C. CHAINAIS-HILLAIRET, S. KRELL, AND A. MOUTON, *Study of discrete duality finite volume schemes for the Peaceman model*, SIAM J. Sci. Comput., 35 (2013), pp. a2928–a2952.

[72] ——, *Convergence analysis of a DDFV scheme for a system describing miscible fluid flows in porous media*, Numer. Methods Partial Differ. Equations, 31 (2015), pp. 723–760.

[73] J. CHANG AND G. COOPER, *A practical difference scheme for fokker-planck equations*, Journal of Computational Physics, 6 (1970), pp. 1–16.

[74] D. CHAUDHURI, J. C. KELLEHER, M. R. O'BRIEN, E. P. O'Reilly, AND S. SCHULZ, *Electronic structure of semiconductor nanostructures: A modified localization landscape theory*, Physical Review B, 101 (2020), p. 035430.

[75] D. CHAUDHURI, M. O'DONOVAN, T. STRECKENBACH, O. MARQUARDT, P. FARRELL, S. K. PATRA, T. KOPRUCKI, AND S. SCHULZ, *Multiscale simulations of the electronic structure of III-nitride quantum wells with varied indium content: Connecting atomistic and continuum-based models*, J. Appl. Phys., 129 (2021), p. 073104.

[76] F. CHAVE, D. A. DI PIETRO, AND S. LEMAIRE, *A discrete Weber inequality on three-dimensional hybrid spaces with application to the HHO approximation of magnetostatics*, Math. Models Methods Appl. Sci., 32 (2022), pp. 175–207.

[77] G. CHEN, P. MONK, AND Y. ZHANG, *An HDG method for the time-dependent drift-diffusion model of semiconductor devices*, J. Sci. Comput., 80 (2019), pp. 420–443.

[78] Y. CHEN AND B. COCKBURN, *Analysis of variable-degree HDG methods for convection-diffusion equations. I: General nonconforming meshes*, IMA J. Numer. Anal., 32 (2012), pp. 1267–1293.

[79] ——, *Analysis of variable-degree HDG methods for convection-diffusion equations. II: Semi-matching nonconforming meshes*, Math. Comput., 83 (2014), pp. 87–111.

[80] H. M. CHENG AND J. T. T. BOONKKAMP, *A generalised complete flux scheme for anisotropic advection-diffusion equations*, Adv. Comput. Math., 47 (2021).

[81] H. M. CHENG, J. TEN THIJE BOONKKAMP, J. JANSSEN, D. MIHAILOVA, AND J. VAN DIJK, *Combining the hybrid mimetic mixed method with the scharfetter-gummel scheme for magnetised transport in plasmas*, 2022.

[82] S. F. CHICHIBU, A. UEDONO, T. ONUMA, B. A. HASKELL, A. CHAKRABORTY, T. KOYAMA, P. T. FINI, S. KELLER, S. P. DENBAARS, J. S. SPECK, U. K. MISHRA, S. NAKAMURA, S. YAMAGUCHI, S. KAMIYAMA, H. AMANO, I. AKASAKI, J. HAN, AND T. SOTA, *Origin of defect-insensitive emission probability in In-containing (Al,In,Ga)N alloy semiconductors*, Nature Mater., 5 (2006), p. 810.

[83] M. CICUTTIN, A. ERN, AND N. PIGNET, *Hybrid high-order methods. A primer with applications to solid mechanics*, SpringerBriefs Math., Cham: Springer, 2021.
[84] B. COCKBURN, D. A. DI PIETRO, AND A. ERN, *Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods*, ESAIM, Math. Model. Numer. Anal., 50 (2016), pp. 635–650.

[85] B. COCKBURN, B. DONG, AND J. GUZMÁN, *A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems*, Math. Comput., 77 (2008), pp. 1887–1916.

[86] B. COCKBURN, B. DONG, J. GUZMÁN, M. RESTELLI, AND R. SACCO, *A hybridizable discontinuous Galerkin method for steady-state convection-diffusion-reaction problems*, SIAM J. Sci. Comput., 31 (2009), pp. 3827–3846.

[87] B. COCKBURN, J. GOPALAKRISHNAN, AND R. LAZAROV, *Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems*, SIAM J. Numer. Anal., 47 (2009), pp. 1319–1365.

[88] L. CODECASA, R. SPECOGNA, AND F. TREVISAN, *A new set of basis functions for the discrete geometric approach*, J. Comput. Phys., 229 (2010), pp. 7401–7410.

[89] P. COLIN, *Numerical analysis of drift-diffusion models: convergence and asymptotic behaviors*, PhD thesis, Université de Lille 1, June 2016.

[90] Y. COUDIÈRE AND F. HUBERT, *A 3d discrete duality finite volume method for nonlinear elliptic equations*, SIAM J. Sci. Comput., 33 (2011), pp. 1739–1764.

[91] M. DARWISH AND A. GAGLIARDI, *A drift-diffusion simulation model for organic field effect transistors: on the importance of the gaussian density of states and traps*, J. Phys. D: Appl. Phys., (2020), p. 105102.

[92] A. DAVID, M. J. GRUNDMANN, J. F. KAEDING, N. F. GARDNER, T. G. MIHOPOULOS, AND M. R. KRAMES, *Carrier distribution in (0001)InGaN/GaN multiple quantum well light-emitting diodes*, Applied Physics Letters, 92 (2008), p. 053502.

[93] A. DAVID, N. G. YOUNG, C. A. HURNI, AND M. D. CRAVEN, *Quantum Efficiency of III-Nitride Emitters: Evidence for Defect-Assisted Nonradiative Recombination and its Effect on the Green Gap*, Phys. Rev. Applied, 11 (2019), p. 031001.

[94] A. DAVID, N. G. YOUNG, C. LUND, AND M. D. CRAVEN, *Review—The Physics of Recombinations in III-Nitride Emitters*, ECS Journal of Solid State Science and Technology, 9 (2019), p. 016021.

[95] P. DAWSON, S. SCHULZ, R. A. OLIVER, M. J. KAPPERS, AND C. J. HUMPHREYS, *The nature of carrier localisation in polar and nonpolar InGaN/GaN quantum wells*, Journal of Applied Physics, 119 (2016), p. 181505.

[96] J.-D. DEUSCHEL AND D. W. STROOCK, *Large deviations*, vol. 137 of Pure and Applied Mathematics, Academic Press, Inc., Boston, MA, 1989.

[97] D. A. DI PIETRO AND J. DRONIOU, *The hybrid high-order method for polytopal meshes. Design, analysis, and applications*, vol. 19 of MS&A, Model. Simul. Appl., Cham: Springer, 2020.

[98] D. A. DI PIETRO, J. DRONIOU, AND A. ERN, *A discontinuous-skeletal method for advection-diffusion-reaction on general meshes*, SIAM J. Numer. Anal., 53 (2015), pp. 2135–2157.

[99] D. A. DI PIETRO AND A. ERN, *A hybrid high-order locking-free method for linear elasticity on general meshes*, Comput. Methods Appl. Mech. Eng., 283 (2015), pp. 1–21.

[100] D. A. DI PIETRO, A. ERN, AND S. LEMAIRE, *An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators*, Comput. Methods Appl. Math., 14 (2014), pp. 461–472.

[101] D. A. DI PIETRO AND S. LEMAIRE, *An extension of the Crouzeix–Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow*, Math. Comp., 84 (2015), pp. 1–31.
[102] A. DI VITO, A. PECCHIA, A. DI CARLO, AND M. AUF DER MAUR, Characterization of non-uniform InGaN alloys: spatial localization of carriers and optical properties, Japanese Journal of Applied Physics, 58 (2019), p. SCCC03.

[103] ——, *Simulating random alloy effects in III-nitride light emitting diodes*, J. Appl. Phys., 128 (2020), p. 041102.

[104] D. H. DOAN, P. FARRELL, J. FUHRMANN, M. KANTNER, T. KOPRUCKI, AND N. ROTUNDO, ddfermi – a drift-diffusion simulation tool, ddfermi – a drift-diffusion simulation tool, Weierstrass Institute (WIAS), doi: http://doi.org/10.20347/WIAS.SOftware.DDFERMI, 2020.

[105] K. DOMELEVO AND P. OMNES, A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids, ESAIM, Math. Model. Numer. Anal., 39 (2005), pp. 1203–1249.

[106] J. DRONIOU, *Non-coercive linear elliptic problems*, Potential Analysis, 17 (2002), pp. 181–203.

[107] ——, *Degrés topologiques et applications*, 2006.

[108] ——, *Remarks on discretizations of convection terms in Hybrid Mimetic Mixed methods*, Networks & Heterogeneous Media, 5 (2010), pp. 545–563.

[109] ——, *Finite volume schemes for diffusion equations: Introduction to and review of modern methods*, Math. Models Methods Appl. Sci., 24 (2014), pp. 1575–1619.

[110] J. DRONIOU AND R. EYMARD, *A mixed finite volume scheme for anisotropic diffusion problems on any grid*, Numer. Math., 105 (2006), pp. 35–71.

[111] J. DRONIOU, R. EYMARD, T. GALLOUËT, C. GUICHARD, AND R. HERBIN, *The Gradient Discretisation Method*, vol. 82 of Mathématiques & Applications, Springer International Publishing, Cham, Switzerland, 2018.

[112] J. DRONIOU, R. EYMARD, T. GALLOUËT, AND R. HERBIN, *A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods*, Math. Models Methods Appl. Sci., 20 (2010), pp. 265–295.

[113] J. DRONIOU AND C. LE POTIER, *Construction and convergence study of schemes preserving the elliptic local maximum principle*, SIAM J. Numer. Anal., 49 (2011), pp. 459–490.

[114] J. DRONIOU AND J.-L. VÁZQUEZ, *Noncoercive convection-diffusion elliptic problems with Neumann boundary conditions*, Calc. Var., 34 (2009), pp. 413–434.

[115] D. A. DUNAVANT, *High degree efficient symmetrical Gaussian quadrature rules for the triangle*, Int. J. Numer. Methods Eng., 21 (1985), pp. 1129–1148.

[116] A. DUSSAIGNE, F. BARBIER, B. DAMILANO, S. CHENOT, A. GRENIER, A. M. PAPON, B. SAMUEL, B. BEN BAKIR, D. VAUFREY, J. C. PILLET, A. GASSE, O. LEDOUX, M. ROZHAVSKAYA, AND D. SOTTA, *Full InGaN red light emitting diodes*, Journal of Applied Physics, 128 (2020), p. 135704.

[117] C. ERIGNOUX, *Hydrodynamic limit of boundary driven exclusion processes with nonreversible boundary dynamics*, J. Stat. Phys., 172 (2018), pp. 1327–1357.

[118] A. ERN AND J.-L. GUERMOND, *Finite Elements I: Approximation and Interpolation*, Springer, Feb. 2021.

[119] L. C. EVANS, *Partial Differential Equations: Second Edition*, vol. 19 of Graduate Studies in Mathematics, American Mathematical Society, Providence, R.I., 2010.

[120] R. EYMARD, J. FUHRMANN, AND K. GÄRTNER, *A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local Dirichlet problems*, Numer. Math., 102 (2006), pp. 463–495.
[121] R. EYMARD, T. GALLOUËT, AND R. HERBIN, *Finite volume methods*, in Techniques of Scientific Computing (Part 3), Handb. Numer. Anal., VII, North-Holland, Amsterdam, 2000, pp. 713–1020.

[122] —, *Benchmark on anisotropic problems. SUSHI: a scheme using stabilization and hybrid interfaces for anisotropic heterogeneous diffusion problems*, in Finite Volumes for Complex Applications V - Problems & Perspectives, R. Eymard and J.-M. Hérard, eds., ISTE, London, 2008, pp. 801–814.

[123] —, *Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces*, IMA J. Numer. Anal., 30 (2010), pp. 1009–1043.

[124] R. EYMARD, C. GUICHARD, AND R. HERBIN, *Small-stencil 3d schemes for diffusive flows in porous media*, ESAIM, Math. Model. Numer. Anal., 46 (2012), pp. 265–290.

[125] P. FARRELL, T. KOPRUCKI, AND J. FUHRMANN, *Computational and analytical comparison of flux discretizations for the semiconductor device equations beyond Boltzmann statistics*, J. Comput. Phys., 346 (2017), pp. 497–513.

[126] P. FARRELL, J. MOATTI, M. O'DONOVAN, S. SCHULZ, AND T. KOPRUCKI, *Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities*, Optical and Quantum Electronics, 55 (2023), p. 978.

[127] P. FARRELL, M. PATRIARCA, J. FUHRMANN, AND T. KOPRUCKI, *Comparison of thermodynamically consistent charge carrier flux discretizations for Fermi–Dirac and Gauss–Fermi statistics*, Optical and Quantum Electronics, 50 (2018), pp. 1–10.

[128] P. FARRELL, N. ROTUNDO, D. H. DOAN, M. KANTNER, J. FUHRMANN, AND T. KOPRUCKI, *Mathematical Methods: Drift-Diffusion Models*, in Handbook of Optoelectronic Device Modeling and Simulation, J. Piprek, ed., vol. 2, CRC Press, 2017, ch. 50, pp. 733–771.

[129] D. FERRY, S. RAMEY, L. SHIFREN, AND R. AKIS, *The effective potential in device modeling: The good, the bad and the ugly*, Journal of Computational Electronics, 1 (2002), pp. 59–65.

[130] D. A. FICK, *On liquid diffusion*, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 10 (1855), pp. 30–39.

[131] F. FILBET AND M. HERDA, *A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure*, Numer. Math., 137 (2017), pp. 535–577.

[132] M. FILOCHE, M. PICCARDO, Y.-R. WU, C.-K. LI, C. WEISBUCH, AND S. MAYBORODA, *Localization landscape theory of disorder in semiconductors. I. Theory and modeling*, Physical Review B, 95 (2017), p. 144204.

[133] R. FINN, M. O'DONOVAN, P. FARRELL, J. MOATTI, T. STRECKENBACH, T. KOPRUCKI, AND S. SCHULZ, *Theoretical study of the impact of alloy disorder on carrier transport and recombination processes in deep UV (Al,Ga)N light emitters*, Applied Physics Letters, 122 (2023). 241104.

[134] R. FINN AND S. SCHULZ, *Impact of random alloy fluctuations on the electronic and optical properties of (Al,Ga)N quantum wells: Insights from tight-binding calculations*, J. Chem. Phys., 157 (2022), p. 244705.

[135] R. A. FISHER, *The wave of advance of advantageous genes*, Annals of Eugenics, 7 (1937), pp. 355–369.

[136] F. FLANDOLI, F. RUSSO, AND J. WOLF, *Some SDEs with distributional drift. I: General calculus*, Osaka J. Math., 40 (2003), pp. 493–542.
[137] —, *Some SDEs with distributional drift. II: Lyons-Zheng structure, Itô's formula and semimartingale characterization*, Random Oper. Stoch. Equ., 12 (2004), pp. 145–184.

[138] J. FOURIER, *Théorie analytique de la chaleur*, Manuscripta; History of science, 18th and 19th century, Chez Firmin Didot, père et fils, 1822.

[139] G. FU, W. QIU, AND W. ZHANG, *An analysis of HDG methods for convection-dominated diffusion problems*, ESAIM, Math. Model. Numer. Anal., 49 (2015), pp. 225–256.

[140] J. FUHRMANN, *Comparison and numerical treatment of generalised Nernst-Planck models*, Comput. Phys. Commun., 196 (2015), pp. 166–178.

[141] J. FUHRMANN, T. STRECKENBACH, ET AL., *pdelib: A finite volume and finite element toolbox for pdes. [software]*., pdelib: a finite volume and finite element toolbox for pdes. [software]. version: 2.4.20190405, Weierstrass Institute (WIAS), http://pdelib.org, 2019.

[142] H. GAJEWSKI, *On existence, uniqueness and asymptotic behavior of solutions of the basic equations for carrier transport in semiconductors*, Z. Angew. Math. Mech., 65 (1985), pp. 101–108.

[143] H. GAJEWSKI AND K. GÄRTNER, *On the discretization of van Roosbroeck's equations with magnetic field*, Z. Angew. Math. Mech., 76 (1996), pp. 247–264.

[144] H. GAJEWSKI AND K. GRÖGER, *Semiconductor equations for variable mobilities based on Boltzmann statistics or Fermi-Dirac statistics*, Math. Nachr., 140 (1989), pp. 7–36.

[145] B. GALLER, A. LAUBSCH, A. WOJCIK, H. LUGAUER, A. GOMEZ-IGLESIAS, M. SABATHIL, AND B. HAHN, *Investigation of the carrier distribution in InGaN-based multi-quantum-well structures*, physica status solidi c, 8 (2011), pp. 2372–2374.

[146] T. GALLOUËT, A. NATALE, AND G. TODESCHI, *From geodesic extrapolation to a variational BDF2 scheme for Wasserstein gradient flows*. working paper or preprint, May 2023.

[147] B. GAUDEUL AND J. FUHRMANN, *Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints*, Numer. Math., 151 (2022), pp. 99–149.

[148] J. GLAD, J. RUSCHEL, N. L. BLOCH, H. K. CHO, F. MAHNKE, L. SULMONI, M. GUTTMANN, T. WERNICKE, M. WEYERS, S. EINFELDT, AND M. KNEISSL, *Impact of operation parameters on the degradation of 233 nm AlGaN-based far UVC LEDs*, J. Appl. Phys., 131 (2022), p. 014501.

[149] A. GLITZKY, M. LIERO, AND G. NIKA, *An existence result for a class of electrothermal drift-diffusion models with Gauss-Fermi statistics for organic semiconductors*, Analysis and Applications, 19 (2021), pp. 275–304.

[150] S. K. GODUNOV AND I. BOHACHEVSKY, *Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics*, Matematičeskij sbornik, 47(89) (1959), pp. 271–306.

[151] S. GOTTLIEB, *On high order strong stability preserving Runge-Kutta and multi step time discretizations*, J. Sci. Comput., 25 (2005), pp. 105–128.

[152] S. GOTTLIEB, D. KETCHESON, AND C.-W. SHU, *Strong stability preserving Runge-Kutta and multistep time discretizations.*, Hackensack, NJ: World Scientific, 2011.

[153] A. GUIONNET AND B. ZEGARLINSKI, *Lectures on Logarithmic Sobolev Inequalities*, Séminaire de probabilités de Strasbourg, 36 (2002), pp. 1–134.

[154] R. HERBIN AND F. HUBERT, *Benchmark on discretization schemes for anisotropic diffusion problems on general grids*, in Finite Volumes for Complex Applications V - Problems & Perspectives, R. Eymard and J.-M. Hérard, eds., ISTE, London, 2008, pp. 659–692.
[155] M. HERDA, *On massless electron limit for a multispecies kinetic system with external magnetic field*, J. Differ. Equations, 260 (2016), pp. 7861–7891.

[156] M. HERDA AND A. ZUREK, *Study of an entropy dissipating finite volume scheme for a nonlocal cross-diffusion system*, ESAIM: M2AN, 57 (2023), pp. 1589–1617.

[157] F. HERMELINE, *A finite volume method for the approximation of diffusion operators on distorted meshes*, Journal of Computational Physics, 160 (2000), pp. 481–499.

[158] C. J. HUMPHREYS, *Solid-State Lighting*, MRS Bulletin, 33 (2008), p. 459.

[159] J.-I. HWANG, R. HASHIMOTO, S. SAITO, AND S. NUNOUE, *Development of InGaN-based red LED grown on (0001) polar surface*, Applied Physics Express, 7 (2014), p. 071003.

[160] A. M. IL'IN, *Differencing scheme for a differential equation with a small parameter affecting the highest derivative*, Mathematical Notes of the Academy of Sciences of the USSR, 6 (1969), pp. 596–602.

[161] E. ISSOGLIO, *Transport equations with fractal noise-existence, uniqueness and regularity of the solution*, Z. Anal. Anwend., 32 (2013), pp. 37–53.

[162] ——, *A non-linear parabolic PDE with a distributional coefficient and its applications to stochastic analysis*, J. Differ. Equations, 267 (2019), pp. 5976–6003.

[163] F. JOCHMANN, *Existence of weak solutions of the drift diffusion model coupled with Maxwell's equations*, J. Math. Anal. Appl., 204 (1996), pp. 655–676.

[164] ——, *Galerkin approximation of weak solutions of the drift diffusion model for semiconductors coupled with Maxwell's equations*, Math. Methods Appl. Sci., 19 (1996), pp. 1471–1488.

[165] ——, *A singular limit in the drift diffusion model for semiconductors coupled with Maxwell's equations*, Appl. Anal., 67 (1997).

[166] ——, *Convergence to stationary states of solutions of the transient drift diffusion equations for semiconductor devices with prescribed currents*, Asymptotic Anal., 18 (1998), pp. 67–109.

[167] ——, *Uniqueness and regularity for the two-dimensional drift-diffusion model for semiconductors coupled with Maxwell's equations*, J. Differ. Equations, 147 (1998), pp. 242–270.

[168] ——, *Some analytical results concerning the drift diffusion equations for semiconductor devices coupled with Maxwell's equations*, in Hyperbolic problems: Theory, numerics, applications. Proceedings of the 7th international conference, Zürich, Switzerland, February 1998. Vol. II, Basel: Birkhäuser, 1999, pp. 525–534.

[169] C. M. JONES, C.-H. TENG, Q. YAN, P.-C. KU, AND E. KIOUPAKIS, *Impact of carrier localization on recombination in InGaN quantum wells and the efficiency of nitride light-emitting diodes: Insights from theory and numerical simulations*, Applied Physics Letters, 111 (2017), p. 113501.

[170] C. JOURDANA, A. JÜNGEL, AND N. ZAMPONI, *Three-species drift-diffusion models for memristors*. working paper or preprint, Apr. 2022.

[171] A. JÜNGEL, *Numerical approximation of a drift-diffusion model for semiconductors with nonlinear diffusion*, ZAMM, 75 (1995), pp. 783–799.

[172] ——, *Quasi-hydrodynamic semiconductor equations*, vol. 41 of Prog. Nonlinear Differ. Equ. Appl., Basel: Birkhäuser, 2001.

[173] ——, *The boundedness-by-entropy method for cross-diffusion systems*, Nonlinearity, 28 (2015), pp. 1963–2001.

[174] ——, *Entropy Methods for Diffusive Partial Differential Equations*, SpringerBriefs in Mathematics, Springer International Publishing, Cham, Switzerland, 2016.
[175] A. JÜNGEL AND P. PIETRA, *A discretization scheme for a quasi-hydrodynamic semiconductor model*, Math. Models Methods Appl. Sci., 7 (1997), pp. 935-955.

[176] A. JÜNGEL AND S. SCHUCHNIGG, *Entropy-dissipating semi-discrete Runge-Kutta schemes for nonlinear diffusion equations*, Commun. Math. Sci., 15 (2017), pp. 27-53.

[177] A. JÜNGEL AND A. ZUREK, *A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system*, SIAM J. Numer. Anal., 59 (2021), pp. 2286-2309.

[178] S. KARPOV, *ABC-model for interpretation of internal quantum efficiency and its droop in III-nitride LEDs: a review.*, Optical and Quantum Electronics, 47 (2015), pp. 1293-1303.

[179] E. F. KELLER AND L. A. SEGEL, *Initiation of slime mold aggregation viewed as an instability*, Journal of Theoretical Biology, 26 (1970), pp. 399-415.

[180] E. KIOUPAKIS, D. STEIAUF, P. RINKE, K. T. DELANEY, AND C. G. V. DE WALLE, *First-principles calculations of indirect Auger recombination in nitride semiconductors*, Phys. Rev. B, 92 (2015), p. 035207.

[181] E. KIOUPAKIS, Q. YAN, D. STEIAUF, AND C. G. VAN DE WALLE, *Temperature and carrier-density dependence of Auger and radiative recombination in nitride optoelectronic devices*, New Journal of Physics, 15 (2013), p. 125006.

[182] M. KNEISSL, T.-Y. SEONG, J. HAN, AND H. AMANO, *The emergence and prospects of deep-ultraviolet light-emitting diode technologies*, Nat. Photonics, 13 (2019), p. 233.

[183] A. KOLMOGOROV, I. PETROVSKII, AND N. PISCOUNOV, *A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem.*, Moscow University Bulletin of Mathematics, 1 (1937).

[184] T. KOPRUCKI AND K. GÄRTNER, *Discretization scheme for drift-diffusion equations with strong diffusion enhancement*, Opt Quant Electron, 45 (2013), pp. 791-796.

[185] ——, *Generalization of the scharfetter-gummel scheme*, in Numerical Simulation of Optoelectronic Devices (NUSOD), 2013 13th International Conference on, 2013, pp. 85-86.

[186] T. KOPRUCKI, N. ROTUNDO, P. FARRELL, D. DOAN, AND J. FUHRMANN, *On thermodynamic consistency of a scharfetter-gummel scheme based on a modified thermal voltage for drift-diffusion equations with diffusion enhancement*, Opt Quant Electron, 47 (2015), pp. 1327-1332.

[187] S. KRELL AND J. MOATTI, *Structure-preserving schemes for drift-diffusion systems on general meshes: DDFV vs HFV*, in Finite Volumes for Complex Applications X - elliptic and parabolic problems, 2023. in press.

[188] P. D. LAX, *Weak solutions of nonlinear hyperbolic equations and their numerical computation*, Commun. Pure Appl. Math., 7 (1954), pp. 159-193.

[189] C. LE BRIS, F. LEGOLL, AND F. MADIOT, *Stabilisation de problèmes non coercifs via une méthode numérique utilisant la mesure invariante (Stabilization of non-coercive problems using the invariant measure)*, Comptes Rendus Mathématique, 354 (2016), pp. 799-803.

[190] C. LEHRENFELD, *Hybrid Discontinuous Galerkin methods for solving incompressible flow problems*, PhD thesis, Rheinisch-Westfälische Technische Hochschule (RWTH) Aachen, 05 2010.

[191] C. LEHRENFELD AND J. SCHÖBERL, *High order exactly divergence-free Hybrid Discontinuous Galerkin methods for unsteady incompressible flows*, Comput. Methods Appl. Mech. Eng., 307 (2016), pp. 339-361.
[192] S. LEMAIRE, *Bridging the hybrid high-order and virtual element methods*, IMA J. Numer. Anal., 41 (2021), pp. 549–593.

[193] G. LHEUREUX, C. LYNSKY, Y. WU, J. SPECK, AND C. WEISBUCH, *A 3d simulation comparison of carrier transport in green and blue c-plane multi-quantum well nitride light emitting diodes*, J Appl Phys, 128 (2020), p. 235703.

[194] C.-K. LI, M. PICCARDO, L.-S. LU, S. MAYBORODA, L. MARTINELLI, J. PERETTI, J. S. SPECK, C. WEISBUCH, M. FILOCHE, AND Y.-R. WU, *Localization landscape theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes*, Phys. Rev. B, 95 (2017), p. 144206.

[195] L. LI AND J.-G. LIU, *Large time behaviors of upwind schemes and B-schemes for Fokker-Planck equations on $\mathbb{R}$ by jump processes*, Mathematics of Computation, (2020).

[196] H. LIU AND Z. WANG, *An entropy satisfying discontinuous Galerkin method for nonlinear Fokker-Planck equations*, J. Sci. Comput., 68 (2016), pp. 1217–1240.

[197] H. LIU, Z. WANG, P. YIN, AND H. YU, *Positivity-preserving third order DG schemes for Poisson-Nernst-Planck equations*, J. Comput. Phys., 452 (2022), p. 22. Id/No 110777.

[198] H. LIU AND H. YU, *Maximum-principle-satisfying third order discontinuous Galerkin schemes for Fokker-Planck equations*, SIAM J. Sci. Comput., 36 (2014), pp. a2296–a2325.

[199] J. P. LIU, J.-H. RYOU, R. D. DUPUIS, J. HAN, G. D. SHEN, AND H. B. WANG, *Barrier effect on hole transport and carrier distribution in InGaN/GaN multiple quantum well visible light-emitting diodes*, Applied Physics Letters, 93 (2008), p. 021102.

[200] Y. LIU AND C.-W. SHU, *Analysis of the local discontinuous Galerkin method for the drift-diffusion model of semiconductor devices*, Sci. China, Math., 59 (2016), pp. 115–140.

[201] Q. LV, J. LIU, C. MO, J. ZHANG, X. WU, Q. WU, AND F. JIANG, *Realization of Highly Efficient InGaN Green LEDs with Sandwich-like Multiple Quantum Well Structure: Role of Enhanced Interwell Carrier Transport*, ACS Photonics, 6 (2019), pp. 130–138.

[202] C. LYNSKY, A. ALHASSAN, G. LHEUREUX, B. BONEF, S. DENBAARS, S. NAKAMURA, Y. WU, C. WEISBUCH, AND J. SPECK, *Barriers to carrier transport in multiple quantum well nitride-based c-plane green light emitting diodes*, Phys Rev Materials, 4 (2020), p. 054604.

[203] C. LYNSKY, G. LHEUREUX, B. BONEF, K. S. QWAH, R. C. WHITE, S. P. DENBAARS, S. NAKAMURA, Y.-R. WU, C. WEISBUCH, AND J. S. SPECK, *Improved Vertical Carrier Transport for Green III-Nitride LEDs Using (In,Ga)N Alloy Quantum Barriers*, Phys. Rev. Applied, 17 (2022), p. 054048.

[204] R. MACNEAL, *An asymmetrical finite difference network*, Quart Math Appl, 11 (1953), pp. 295–310.

[205] F. MADIOT, *Multiscale finite element methods for advection-diffusion problems*, PhD thesis, Université Paris-Est, December 2016.

[206] A. J. MAJDA AND A. L. BERTOZZI, *Vorticity and incompressible flow*, Camb. Texts Appl. Math., Cambridge: Cambridge University Press, 2002.

[207] P. A. MARKOWICH, *The stationary semiconductor device equations*, Computational Microelectronics, Springer-Verlag, Vienna, 1986.

[208] P. A. MARKOWICH, C. A. RINGHOFER, AND C. SCHMEISER, *Semiconductor equations*, Wien: Springer-Verlag, 1990.

[209] P. A. MARKOWICH AND A. UNTERREITER, *Vacuum solutions of a stationary drift-diffusion model*, Annali della Scuola Normale Superiore di Pisa - Classe di Scienze, Ser. 4, 20 (1993), pp. 371–386.
[210] J. M. McMAHON, E. KIOUPAKIS, AND S. SCHULZ, *Atomistic analysis of Auger recombination in c-plane (In,Ga)N/GaN quantum wells: Temperature-dependent competition between radiative and nonradiative recombination*, Physical Review B, 105 (2022), p. 195307.

[211] J. M. McMAHON, D. S. P. TANNER, E. KIOUPAKIS, AND S. SCHULZ, *Atomistic analysis of radiative recombination rate, Stokes shift and density of states in c-plane InGaN/GaN qunatum wells*, Appl. Phys. Lett., 116 (2020), p. 181104.

[212] J. MOATTI, *A skeletal high-order structure preserving scheme for advection-diffusion equations*, in Finite Volumes for Complex Applications X - elliptic and parabolic problems, 2023. in press.

[213] ——, *A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes*, ESAIM: M2AN, 57 (2023), pp. 2557-2593.

[214] M. S. MOCK, *An initial value problem from semiconductor device theory*, SIAM J. Math. Anal., 5 (1974), pp. 597-612.

[215] ——, *An example of nonuniqueness of stationary solutions in semiconductor device models*, COMPEL-The international journal for computation and mathematics in electrical and electronic engineering, 1 (1982), pp. 165-174.

[216] ——, *Analysis of mathematical models of semiconductor devices*, vol. 3 of Adv. Numer. Comput. Ser., Boole Press, Dublin, 1983.

[217] G. MUZIOL, H. TURSKI, M. SIEKACZ, K. SZKUDLAREK, L. JANICKI, M. BARANOWSKI, S. ZOLUD, R. KUDRAWIEC, T. SUSKI, AND C. SKIERBISZEWSKI, *Beyond Quantum Efficiency Limitations Originating from the Piezoelectric Polarization in Light-Emitting Devices*, ACS Photonics, 6 (2019), pp. 1963-1971.

[218] W. NERNST, *Die elektromotorische wirksamkeit der jonen*, Zeitschrift für physikalische Chemie, 4 (1889), pp. 129-181.

[219] M. O'DONOVAN, P. FARRELL, J. MOATTI, T. STRECKENBACH, T. KOPRUCKI, AND S. SCHULZ, *Impact of random alloy fluctuations on the carrier distribution in multi-color (In,Ga)N/GaN quantum well systems*. working paper or preprint, Sept. 2022.

[220] M. O'DONOVAN, M. LUISIER, E. P. O'Reilly, and S. SCHULZ, *Impact of random alloy fluctuations on inter-well transport in InGaN/GaN multi-quantum well systems: an atomistic non-equilibrium Green's function study*, J. Phys.: Condens. Matter, 33 (2021), p. 045302.

[221] M. O'DONOVAN, D. CHAUDHURI, T. STRECKENBACH, P. FARRELL, S. SCHULZ, AND T. KOPRUCKI, *From atomistic tight-binding theory to macroscale drift-diffusion: Multiscale modeling and numerical simulation of uni-polar charge transport in (In,Ga)N devices with random fluctuations*, Journal of Applied Physics, 130 (2021), p. 065702.

[222] M. O'DONOVAN, P. FARRELL, T. STRECKENBACH, T. KOPRUCKI, AND S. SCHULZ, *Multiscale simulations of uni-polar hole transport in (In,Ga)N quantum well systems*, Optical and Quantum Electronics, 54 (2022), p. 405.

[223] G. PAASCH AND S. SCHEINERT, *Charge carrier density of organics with Gaussian density of states: Analytical approximation for the Gauss-Fermi integral*, Journal of Applied Physics, 107 (2010), 104501.

[224] P. PELLETIER, D. DELANDE, V. JOSSE, A. ASPECT, S. MAYBORODA, D. N. ARNOLD, AND M. FILOCHE, *Spectral functions and localization-landscape theory in speckle potentials*, Phys. Rev. A, 105 (2022), p. 023314.

[225] M. PICCARDO, C.-K. LI, Y.-R. WU, J. S. SPECK, B. BONEF, R. M. FARRELL, M. FILOCHE, L. MARTINELLI, J. PERETTI, AND C. WEISBUCH, *Localization landscape theory of disorder in semiconductors. II. Urbach tails of disordered quantum well layers*, Phys. Rev. B, 95 (2017), p. 144205.
[226] J. PIPREK, *Efficiency droop in nitride-based light-emitting diodes*, Phys. Status Solidi A, 207 (2010), pp. No. 10, 2217-2225.

[227] M. PLANCK, *Ueber die erregung von electricität und wärme in electrolyten*, Annalen der Physik, 275 (1890), pp. 161-186.

[228] O. PURBO, D. CASSIDY, AND S. CHISHOLM, *Numerical model for degenerate and heterostructure semiconductor devices*, J Appl Phys, 66 (1989), pp. 5078-5082.

[229] W. QIU AND K. SHI, *An HDG method for convection diffusion equation*, J. Sci. Comput., 66 (2016), pp. 346-357.

[230] J. QUASTEL, *Diffusion of color in the simple exclusion process*, Commun. Pure Appl. Math., 45 (1992), pp. 623-679.

[231] K. S. QWAH, M. MONAVARIAN, G. LHEUREUX, J. WANG, Y.-R. WU, AND J. S. SPECK, *Theoretical and experimental investigations of vertical hole transport through unipolar AlGaN structures: Impacts of random alloy disorder*, Applied Physics Letters, 117 (2020), p. 022107.

[232] X. REN, Z. WANG, W. SHA, AND W. CHOY, *Exploring the way to approach the efficiency limit of perovskite solar cells by drift-diffusion model*, ACS Photonics, 4 (2017), pp. 934-942.

[233] N. ROCCATO, F. PIVA, C. D. SANTI, M. BUFFOLO, M. FREGOLENT, M. PILATI, N. SUSILO, D. H. VIDAL, A. MUHIN, L. SULMONI, T. WERNICKE, M. KNEISSL, G. MENEGHESSO, E. ZANONI, AND M. MENEGHINI, *Modeling the electrical degradation of AlGaN-based UV-C LEDs by combined deep-level optical spectroscopy and TCAD simulations*, Appl. Phys. Lett., 122 (2023), p. 161105.

[234] R. RODRIGUEZ-TORRES, E. GUTIERREZ-DOMINGUEZ, R. KLIMA, AND S. SELBERHERR, *Analysis of split-drain magfets*, IEEE Transactions on Electron Devices, 51 (2004), pp. 2237-2245.

[235] F. RÖMER, M. GUTTMANN, T. WERNICKE, M. KNEISSL, AND B. WITZIGMANN, *Effect of inhomogeneous Broadening in Ultraviolet III-Nitride Light-Emitting Diodes*, Materials, 14 (2021), p. 7890.

[236] F. RÖMER AND B. WITZIGMANN, *Luminescence distribution in the multi-quantum well region of iii-nitride light emitting diodes*, Proc. of SPIE, 10124 (2017), pp. 101240Y-1.

[237] ——, *Signature of the ideality factor in iii-nitride multi quantum well light emitting diodes*, Opt Quant Electron, 50 (2018), pp. 1-10.

[238] F. RUSSO AND G. TRUTNAU, *Some parabolic PDEs whose drift is an irregular random noise in space*, Ann. Probab., 35 (2007), pp. 2213-2262.

[239] D. L. SCHARFETTER AND H. K. GUMMEL, *Large-signal analysis of a silicon Read diode oscillator*, IEEE Transactions on Electron Devices, 16 (1969), pp. 64-77.

[240] D. SCHIAVON, M. BINDER, M. PETER, B. GALLER, P. DRECHSEL, AND F. SCHOLZ, *Wavelength-dependent determination of the recombination rate coefficients in single-quantum-well GaInN/GaN light emitting diodes*, physica status solidi (b), 250 (2013), pp. 283-290.

[241] M. SCHNEIDER, L. AGÉLAS, G. ENCHÉRY, AND B. FLEMISCH, *Convergence of nonlinear finite volume schemes for heterogeneous anisotropic diffusion on general meshes*, J. Comput. Phys., 351 (2017), pp. 80-107.

[242] S. SCHULZ, M. A. CARO, C. COUGHLAN, AND E. P. O'Reilly, *Atomistic analysis of the impact of alloy and well width fluctuations on the electronic and optical properties of InGaN/GaN quantum wells*, Phys. Rev. B, 91 (2015), p. 035439.

[243] S. SCHULZ, M. A. CARO, AND E. P. O'Reilly, *Impact of cation-based localized electronic states on the conduction and valence band structure of Al₁-xInxN alloys*, Appl. Phys. Lett., 104 (2014), p. 172102.
[244] H.-T. SHEN, Y.-C. CHANG, AND Y.-R. WU, *Analysis of Light-Emission Polarization Ratio in Deep-Ultraviolet Light-Emitting Diodes by Considering Random Alloy Fluctuations with the 3D k.p Method*, Phys. Status Solidi, 16 (2022), p. 2100498.

[245] Z. SHENG, J. YUE, AND G. YUAN, *Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes*, SIAM J. Sci. Comput., 31 (2009), pp. 2915-2934.

[246] N. SHIGESADA, K. KAWASAKI, AND E. TERAMOTO, *Spatial segregation of interacting species*, Journal of Theoretical Biology, 79 (1979), pp. 83-99.

[247] H. SI, *TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator*, ACM Transactions on Mathematical Software, 41 (2015), pp. 1-36.

[248] H. SI, K. GÄRTNER, AND J. FUHRMANN, *Boundary conforming delaunay mesh generation*, Comput Math Math Phys, 50 (2010), pp. 38-53.

[249] S. STODTMANN, R. LEE, C. WEILER, AND A. BADINSKI, *Numerical simulation of organic semiconductor devices with high carrier densities*, J Appl Phys, 112 (2012), p. 114909.

[250] S. SU AND H. TANG, *A positivity-preserving and free energy dissipative hybrid scheme for the Poisson-Nernst-Planck equations on polygonal and polyhedral meshes*, Comput. Math. Appl., 108 (2022), pp. 33-48.

[251] M. SZYMAŃSKI, D. TU, AND R. FORCHHEIMER, *2-d drift-diffusion simulation of organic electrochemical transistors*, IEEE Transactions on Electron Devices, 64 (2017), pp. 5114-5120.

[252] D. S. P. TANNER, M. A. CARO, E. P. O'ReILLY, AND S. SCHULZ, *Random alloy fluctuations and structural inhomogeneities in c-plane $In_xGa_{1-x}N$ quantum wells: theory of ground and excited electron and hole states*, RSC Adv., 6 (2016), pp. 64513-64530.

[253] D. S. P. TANNER, P. DAWSON, M. J. KAPPERS, R. A. OLIVER, AND S. SCHULZ, *Polar InGaN/GaN quantum wells: Revisiting the impact of carrier localization on the "green gap" problem*, Phys. Rev. Applied, 13 (2020), p. 044068.

[254] G. TOSCANI, *Entropy production and the rate of convergence to equilibrium for the Fokker-Planck equation*, Quart. Appl. Math., 57 (1999), pp. 521-541.

[255] W. TRESS, K. LEO, AND M. RIEDE, *Optimum mobility, contact properties, and open-circuit voltage of organic solar cells: A drift-diffusion simulation study*, Phys Rev B, (2012), p. 155201.

[256] T.-Y. TSAI, K. MICHALCZEWSKI, P. MARTYNIUK, C.-H. WU, AND Y. WU, *Application of localization landscape theory and the k-p model for direct modeling of carrier transport in a type ii superlattice inas/inassb photoconductor system*, J Appl Phys, 127 (2020), p. 033104.

[257] S. L. M. VAN MENSFOORT AND R. COEHOORN, *Effect of gaussian disorder on the voltage dependence of the current density in sandwich-type devices based on organic semiconductors*, Phys. Rev. B, 78 (2008), p. 085207.

[258] W. VAN ROOSBROECK, *Theory of the flow of electrons and holes in germanium and other semiconductors*, The Bell System Technical Journal, 29 (1950), pp. 560-607.

[259] M. WAHN AND J. NEUGEBAUER, *Generalized Wannier functions: An efficient way to construct ab-initio tight-binding parameters for group-III nitrides*, phys. stat. sol. (b), 243 (2006), pp. No. 7, 1583-1587.

[260] D. WATSON-PARRIS, M. J. GODFREY, P. DAWSON, R. A. OLIVER, M. J. GALTREY, M. J. KAPPERS, AND C. J. HUMPHREYS, *Carrier localization mechanisms in $In_xGa_{1-x}N/GaN$*, Phys. Rev. B, 83 (2011), p. 115321.

[261] C. WEISBUCH, S. NAKAMURA, Y.-R. WU, AND J. S. SPECK, *Disorder effects in nitride semiconductors: impact on fundamental and device properties*, Nanophotonics, 10(1) (2021), pp. 3-21.
[262] B. WITZIGMANN, F. RÖMER, M. MARTENS, C. KUHN, T. WERNICKE, AND M. KNEISSL, *Calculation of optical gain in AlGaN quantum wells for ultraviolet emission*, AIP Advances, 10 (2020), p. 095307.

[263] C.-K. WU, C.-K. LI, AND Y.-R. WU, *Percolation transport study in nitride based LED by considering the random alloy fluctuation*, J. Comp. Elec., 14 (2015), pp. 416–424.

[264] T.-J. YANG, R. SHIVARAMAN, J. S. SPECK, AND Y.-R. WU, *The influence of random indium alloy fluctuations in indium gallium nitride quantum wells on the device behavior*, J. Appl. Phys., 116 (2014), p. 113104.

[265] Z. YU AND R. DUTTON, *SEDAN III – A one-dimensional device simulator*. www-tcad.stanford.edu/tcad/programs/sedan3.html, 1988.

[266] J. H. ZHU, S. M. ZHANG, H. WANG, D. G. ZHAO, J. J. ZHU, Z. S. LIU, D. S. JIANG, Y. X. QIU, AND H. YANG, *The investigation on carrier distribution in InGaN/GaN multiple quantum well layers*, Journal of Applied Physics, 109 (2011), p. 093117.
Table des matières

<table><tr><td>Résumé</td><td>xi</td></tr><tr><td>Sommaire</td><td>xiii</td></tr><tr><td>Liste des tableaux</td><td>xvii</td></tr><tr><td>Table des figures</td><td>xix</td></tr><tr><td>Avant-propos</td><td>1</td></tr><tr><td>Introduction générale</td><td>3</td></tr><tr><td>1 Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches</td><td>31</td></tr><tr><td>2 A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes</td><td>77</td></tr><tr><td>3 A comparison of structure-preserving schemes for drift-diffusion systems on general meshes: DDFV vs HFV</td><td>117</td></tr><tr><td>4 High-order polytopal schemes for advection-diffusion equations: linear and nonlinear approaches</td><td>127</td></tr><tr><td>Conclusion and perspectives</td><td>161</td></tr><tr><td>A Semiconductor models with varying band edge energies: an overview</td><td>167</td></tr><tr><td>B Impact of random alloy fluctuations on the carrier distribution in multi-color (In,Ga)N/GaN quantum well systems</td><td>173</td></tr><tr><td>C Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities</td><td>189</td></tr><tr><td>D Theoretical study of the impact of alloy disorder on carrier transport and recombination processes in deep UV (Al,Ga)N light emitters</td><td>199</td></tr><tr><td>Bibliography</td><td>207</td></tr><tr><td>Table des matières</td><td>223</td></tr></table>
<figure><img src="image_64.png" /></figure>
Résumé

Dans cette thèse, nous nous intéressons à l'approximation numérique de problèmes de convection-diffusion, potentiellement anisotropes, par des schémas d'ordre élevé sur maillages généraux. Notre objectif est de proposer des méthodes fiables, précises et efficaces : les solutions numériques doivent préserver les propriétés physiques des solutions continues (conservation de la masse, positivité des densités, comportement en temps long) tout en autorisant une large gamme de paramètres de discrétisation (pas de temps grands, maillages spatiaux généraux) et en optimisant la précision de calcul à coût donné. Les problèmes considérés sont des équations d'advection-diffusion ainsi que des systèmes couplés de dérive-diffusion qui modélisent les composants semi-conducteurs.

On se concentre d'abord sur une équation d'advection-diffusion seule, pour laquelle nous proposons et analysons trois méthodes d'ordre bas de type volumes finis hybrides (HFV). Cette comparaison met en avant la nécessité d'utiliser un schéma non-linéaire afin de préserver la positivité des solutions, tant d'un point de vue théorique que numérique. On s'intéresse alors à l'approximation d'un système de dérive-diffusion, constitué de deux équations d'advection-diffusion couplées avec une équation de Poisson. Pour ce problème, on introduit un schéma non-linéaire basé sur la méthode précédente, qui préserve les bornes des densités calculées (et en particulier leur positivité) et le comportement en temps long de la solution. Ce schéma HFV pour la dérive-diffusion est ensuite comparé numériquement avec un schéma présentant des propriétés similaires basé sur la méthode volumes finis en dualité discrète (DDFV).

Nous nous intéressons alors à des schémas d'ordre élevé (en espace). Ces schémas sont basés sur des méthodes de type hybride d'ordre élevé (HHO) qui peuvent être interprétées comme des extensions à l'ordre arbitraire des méthodes HFV. Nous introduisons deux méthodes pour les équations d'advection-diffusion linéaires. La première est linéaire, tandis que la seconde est non-linéaire, et permet de préserver la positivité. Pour ces deux schémas, nous prouvons l'existence de solutions discrètes et établissons des résultats de comportement en temps long. Nous confirmons également ces résultats numériquement, et mettons en avant la nécessité d'utiliser une méthode non-linéaire pour préserver la positivité. Par ailleurs, on observe que la méthode non-linéaire converge à l'ordre attendu. De plus, la montée en ordre permet un gain d'efficacité (précision de l'approximation pour un coût de calcul donné) conséquent par rapport aux méthodes d'ordre bas des premiers chapitres.

Ces travaux sont complétés par l'étude de problèmes de convection-diffusion avec convection très irrégulière, effectuée en collaboration avec des physiciens durant la thèse. Les recherches menées visent à comprendre comment concevoir des diodes électroluminescentes efficaces émettant dans l'ultraviolet profond, et soulèvent divers enjeux relatifs à la modélisation, l'analyse numérique et la simulation de ces problèmes.

**Mots clés :** analyse numérique, méthode d'entropie, méthodes numériques sur maillages généraux, méthodes numériques d'ordre élevé, semi-conducteurs, edp paraboliques
Abstract

In this thesis, we are interested in the numerical approximation of anisotropic convection-diffusion problems using high-order schemes on general meshes. Our objective is to develop reliable, accurate, and efficient methods: the numerical solutions must preserve the physical properties of the continuous solutions (mass conservation, positivity of densities, long-time behaviour) while allowing for a wide range of discretisation parameters (large time steps, general spatial meshes) and optimising the accuracy at a given computational cost. The problems under study are advection-diffusion equations as well as coupled drift-diffusion systems that model semiconductor devices.

We first focus on a single advection-diffusion equation for which we propose and analyse three different hybrid finite volume (HFV) methods to approximate the solution. Their comparison highlights the necessity of using a nonlinear scheme to ensure the preservation of positivity, both from a theoretical and a numerical level. We then consider the numerical approximation of a drift-diffusion system, consisting of two coupled advection-diffusion equations with a Poisson equation. For this problem, we introduce a nonlinear scheme based on the previous method, which preserves the bounds on the computed densities (including their positivity) and the long-time behaviour of the solution. This HFV scheme for drift-diffusion is then numerically compared with another scheme with similar properties based on the discrete duality finite volume (DDFV) method.

We then focus on high-order (in space) schemes. These schemes are based on hybrid high-order (HHO) methods which can be interpreted as an arbitrary-order extension of HFV methods. We introduce two schemes for linear advection-diffusion: the first method is linear, while the second one is nonlinear and preserves the positivity of the solution. For both of these schemes, we prove the existence of discrete solutions and establish results about their long-time behaviour. We also confirm these results numerically and emphasise the need to use a nonlinear method to preserve positivity. We observe that the nonlinear method converges at the expected order. Furthermore, increasing the order leads to a significant gain in efficiency (accuracy of the approximation for a given computational cost) compared to the low-order methods discussed in the first chapters.

These works are complemented with the study of convection-diffusion problems with highly irregular convection, carried out in collaboration with physicists during the thesis. These investigations aim to understand how to design efficient deep ultraviolet-emitting electroluminescent diodes, and raise various issues related to the modeling, numerical analysis and simulation of these problems.

**Keywords:** numerical analysis, entropy method, numerical methods on general meshes, high-order numerical methods, semiconductors, parabolic pdes