Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities

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.


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) (Roemer et al. 2018;Li et al. 2017;O'Donovan et al. 2022), transistors (Szymaǹski et al. 2017) (Darwish and Gagliardi 2020) and solar cells (Ren et al. 2017;Tress et al. 2012).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 (Bessemoulin-Chatard 2012).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 (Humphreys 2008) and for which considerable effort has been undertaken to develop advanced carrier transport models (O'Donovan et al. 2021).
For the Boltzmann distribution, the classical Scharfetter-Gummel scheme (Scharfetter and Gummel 1969) 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 (Purbo et al. 1989;Jüngel 1995;Stodtmann et al. 2012).Unfortunately, these schemes are not thermodynamically consistent.In Koprucki and Gärtner (2013a), Koprucki and Gärtner (2013b) a thermodynamically consistent generalization for Blakemore statistics (which is itself a special case of Eymard et al. 2006) is presented in the spirit of (Scharfetter and Gummel 1969) 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 Chainais-Hillairet et al. (2022).On the other hand, in Bessemoulin-Chatard (2012) 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 Koprucki et al. (2015a).Finally, the SEDAN scheme (Yu and Dutton 1988;Cancès et al. 2020;Abdel 2021) 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 Sect.2, we describe the bipolar drift-diffusion model for charge transport in semiconductors.Its finite volume discretization including the flux discretizations is described in Sect.3. The formal definition of discrete thermodynamic consistency is also presented in this section.We compare thermodynamically consistent and inconsistent schemes in Sect. 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 Sect. 5.

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 Farrell et al. (2017a).The dependence of the carrier densities n and p on the chemical potentials for electrons and holes n and p are described by a statistical distribution function F as well as conduction and valence band densities of states N c and N v via the state equations n = N c F( n ) and p = N v F( p ) .Typical choices for the distribution function are F( ) = exp( ) , the so- called Boltzmann approximation, or e E− +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 n and p via Here q denotes the elementary charge, 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 .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 (Ferry et al. 2002;Li et al. 2017).One of these approaches is based on the so-called localization landscape theory (LLT), which allows to extract a (nonlocal) 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 (Chaudhuri et al. 2020).The resulting (effective) confining potentials, E eff c,v , exhibit band edges that 978 Page 4 of 12 are softened and approximate the finite extent of carrier wavefunctions. 1 In the framework of LLT a linear system of equations is solved to determine the effective potential without introducing extra free parameters (Filoche et al. 2017) We model a bipolar semiconductor device as a domain Ω ⊂ ℝ 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 and continuity equations for electrons and holes: Here, 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 or for any strictly monotonic Fermi-like distribution function F( ) in drift-diffusion form where n and 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(F −1 ) � (x) , for x ∈ ℝ .This factor is the so-called diffusion enhancement appearing as a density- dependent modification factor in the generalized Einstein relation, see (van Mensfoort and Coehoorn 2008), leading in general to a non-linear diffusion coefficient.For the Boltzmann distribution, F( ) = exp( ) , we have g ≡ 1 and the current expressions (4) reduce to the usual ones with linear diffusion. (1) (3) j n = −q n n∇ n , j p = −q p p∇ p , (4) 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, .

Finite volume space discretization
We discretize the domain Ω using the Voronoï box based finite volume method introduced in Macneal (1953), also known as "box method" due to Bank and Rose (1987).It uses a simplical boundary conforming Delaunay grid (Si et al. 2010) which allows to obtain control volumes surrounding each given collocation point x K by joining the circumcenters of the simplices containing it, see (Farrell et al. 2017a) for details.
Let K denote the boundary of the control volume K, and | | the measure of a geometrical object .For each control volume K, we integrate the continuity Eq. 2) and apply Gauss's theorem to the integral of the flux divergence.Restricting our considerations to the electron transport equation, we obtain where n is the internal unit normal to K and n KL is the internal unit normal to the interface K ∩ 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 K ∩ 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 Farrell et al. (2017a).

Discrete thermodynamic consistency
One property which holds on a continuous level to avoid unphysical state dissipation is the preservation of thermodynamic equilibrium (Farrell et al. 2017a).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 Koprucki et al. 2015a;Farrell et al. 2017a): a numerical flux j = j KL is said to be thermodynamically consistent if it satisfies an analogous discrete relation, i.e. where . We point out that the condition (6) holds in equilibrium.Here, we introduce a stronger notion of thermodynamic consistency, which holds outside of equilibrium, namely 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 (3).Thermodynamic consistency is also important, when coupling the van Roosbroeck system to heat transport models (Farrell et al. 2017a).We discuss now different numerical fluxes that may be used within a Voronoï finite volume framework. (5) ) j = 0 implies KL = 0, (7) j ≤ 0 implies KL ≥ 0 and j ≥ 0 implies KL ≤ 0.

The Scharfetter-Gummel scheme
First, we introduce the well known, classical Scharfetter-Gummel flux approximation (Scharfetter and Gummel 1969) given by where the constant j 0 is given by j It is important to point out that Scharfetter and Gummel introduced this numerical flux only in the Boltzmann regime, i.e.F = exp .In this case, the flux is thermodynamically consistent in the sense of (6).However, once we leave the Boltzmann regime, i.e.F ≠ exp , and continue using (8) this numerical flux will no longer be thermodynamically consistent.

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 a excess chemical potential scheme is the source code of the SEDAN III simulator (Yu and Dutton 1988), which explains the reason we use this name.A numerical analysis focused comparison of this flux approximation is given in Cancès et al. (2020) and simulation results are presented in Abdel (2021).The scheme is motivated by rearranging the drift part to include the excess chemical potential, ex = ln F( ) − , yielding with 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 Therefore, it is easy to see that the SEDAN flux satisfies (7), since both B(Q KL ) and F( L ) are positive: it is a thermodynamicaly consistent numerical flux.
Note that when applying Boltzmann statistics 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.

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 O' Donovan et al. (2022).For large negative (8) values of ( ≤ −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. O'Donovan 2022).
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 1a displays the conduction band edge and quasi Fermi energies of the MQW system when treated using Boltzmann statistics at a bias of 3.3 V. Applying the Scharfetter-Gummel scheme (8) to Fermi-Dirac statistics for a bias of 3.3 V leads to humps in the quasi Fermi energy within each QW, see Fig. 1b.
Recalling the condition (7), for thermodynamically consistent schemes, the discrete gradient of the discrete quasi Fermi level should indicates the direction of the discrete electron flow (as it is in the continuous case, according to (3)).However in Fig. 1b, 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 (7).Moreover, as there is no generation of carriers in the system this change in direction of the electron flux is highly unphysical.
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 (Fig. 1a).Similarly, using Fermi-Dirac statistics with the SEDAN scheme does not exhibit the unphysical humps in the Fermi level (Fig. 1c).
Another perspective to interpret thermodynamic inconsistency is to note the incorrect interplay between quasi Fermi energies and local current fluxes.We see in Fig. 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 (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 Fig. 1a, c.However, for the inconsistent case the Previous studies of the numerical flux approximations have shown that a thermodynamically inconsistent scheme can result in the incorrect sign of the particle flux (Koprucki et al. 2015a).This is also reflected by the fact that (6) holds only for consistent schemes, which guarantee the physically correct sign of the current also far from equilibrium.
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.3 V.The Fermi-Dirac function grows like a polynomial while the Boltzmann approximation grows for large densities exponentially, see for example Figure 50.9 in Farrell et al. (2017a).This different behaviour leads to nonlinear diffusion, the so-called diffusion enhancement, for non-Boltzmann statistics, see (4), or Koprucki et al. (2015a).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. 3a) and using Fermi-Dirac statistics with the Scharfetter-Gummel scheme (Fig. 3b) 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 Fig. 1b.Comparing these densities with the correctly calculated density using Fermi-Dirac statistics in combination with the thermodynamically consistent SEDAN scheme (Fig. 3c) 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.3 V.
Because it influences the carrier density, 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 using Fermi-Dirac statistics with the SEDAN scheme (as the thermodynamically inconsistent SG scheme with Fermi-Dirac statistics results in the same densities as the Boltzmann case this is equivalent to comparing the SG and SEDAN schemes with Fermi-Dirac statistics).This is highlighted in Fig. 4a, where the differences between Fermi-Dirac and Boltzmann-like behaviour are shown for the three recombination rates, calculated via where R s i is the recombination rate associated with the process i (i ∈ {Shockley-Read-Hall, Radiative, Auger}) calculated with the scheme s (s ∈ {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.3 V.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 Fig. 4b, 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 order of magnitude at 3.6 V.
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 (Bandelow et al. 2005), expressed by the so-called Fermi voltage.In this case the unphysical humps seen in Fig. 1b will lead to an incorrect prediction of the transparency density.

Conclusion
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.

Fig. 1
Fig.1Conduction band edge (black) and quasi Fermi energy (red) at a bias of 3.3 V 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

Fig. 2 Fig. 3
Fig. 2 Numerical electron flux averaged over each atomic plane at a bias of 3.3 V shown for Boltzmann statisitcs (black), Fermi-Dirac statistics incorrectly using the Scharfetter-Gummel (SG) flux discretization (red) and Fermi-Dirac statistics correctly using the SEDAN flux discretization (blue).(Colorfigure online)

Fig. 4 a
Fig. 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.3 V, 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.(Colorfigure online) (Tsai et al. 2020)on Lheureux et al. 2020;l.(2021withadrift-diffusionbasedcarriertransport solver is discussed in more detail inO'Donovan et al. (2022) and O'Donovan et al. (2021).In O'Donovan et al. (2022)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 semiconductor structures including LEDs(Lynskey et al. 2020;Lheureux et al. 2020; Lynskey et al. 2022; Romer  and Witzigmann 2017), blocking layers(Qwah et al. 2020) and superlattices(Tsai et al. 2020); moreover, these methods have been used successfully alongside experimental studies to gain insight into device behaviour(Lynskey et al. 2022; Romer and Witzigmann  2017).