Assessing the quality of the excess chemical potential flux scheme for degenerate semiconductor device simulation

The van Roosbroeck system models current flows in (non-)degenerate semiconductor devices. Focusing on the stationary model, we compare the excess chemical potential discretization scheme, a flux approximation which is based on a modification of the drift term in the current densities, with another state-of-the-art Scharfetter–Gummel scheme, namely the diffusion-enhanced scheme. Physically, the diffusion-enhanced scheme can be interpreted as a flux approximation which modifies the thermal voltage. As a reference solution we consider an implicitly defined integral flux, using Blakemore statistics. The integral flux refers to the exact solution of a local two point boundary value problem for the continuous current density and can be interpreted as a generalized Scharfetter–Gummel scheme. All numerical discretization schemes can be used within a Voronoi finite volume method to simulate charge transport in (non-)degenerate semiconductor devices. The investigation includes the analysis of Taylor expansions, a derivation of error estimates and a visualization of errors in local flux approximations to extend previous discussions. Additionally, drift-diffusion simulations of a p–i–n device are performed.


Introduction
The standard drift-diffusion model for semi-classical charge transport of free electrons and holes due to a self-consistent electric field in a semiconductor device is the van Roosbroeck system. We consider Voronoi finite volume schemes for the discretization of the semiconductor device equations. We are interested in the numerically more challenging degenerate case and pay particular attention to the choice of flux approximations. One of the most classical numerical flux scheme is arguably the Scharfetter-Gummel scheme (Scharfetter and Gummel 1969), which yields a numerical stable and thermodynamically consistent numerical flux, but cannot be used for general charge carrier statistics. A generalization of the Scharfetter-Gummel scheme is available (Eymard et al. 2006), but computationally expensive. Hence, several thermodynamically consistent numerical flux schemes, which modify this generalization to lower the computational costs, are proposed in the literature (Farrell et al. 2017bPatriarca et al. 2018). We focus on the excess chemical potential scheme (Yu and Dutton 1988) which appears to be used in parts of the device simulation community (Silvaco International 2016; Synopsys, Inc 2010). However, unfortunately, there seem to be no direct comparisons of this scheme with other recent modified Scharfetter-Gummel schemes. This paper aims to fill this gap by comparing it to the diffusion enhanced scheme (Bessemoulin-Chatard 2012), which was already compared to other modified Scharfetter-Gummel schemes and seems to be the most promising modified Scharfetter-Gummel flux (Farrell et al. 2017b). The integral flux by Eymard et al. is used as a reference flux. It supplements previously made temperaturedependent observations (Kantner 2020).

Van Roosbroeck model
The stationary variant of the van Roosbroeck system is given by where q denotes the elementary charge, s the dielectric permittivity, N A and N D describe the density of singly ionized acceptor and donor atoms, and R the recombination. The set of unknowns is expressed by the electrostatic potential and the quasi Fermi potentials for electrons n and holes p . The current densities in the continuity Eqs. (1b) and (1c) are given by where the electron and hole densities n and p are defined by (Sze and Ng 2006) (2) n = −q n n∇ n , p = −q p p∇ p , The strictly monotonously increasing statistics function F will be discussed later. The conduction and valence band density of states are given by N c and N v , the mobilities by n and p and the Boltzmann constant by k B . The conduction and valence band-edge energies are denoted by E c and E v and T refers to the temperature. With help of the generalized Einstein relation it is possible to model the diffusion coefficients D n and D p via the the nonlinear diffusion enhancement by (introducing the thermal voltage U T = k B T∕q) With this relation, we can rewrite the electric fluxes (2) into a drift-diffusion form In general, inorganic semiconductor devices can be modeled by choosing the Fermi-Dirac integral of order one-half (Sze and Ng 2006) for the statistics function F . Non-degenerate semiconductors are modeled with the Boltzmann approximation F( ) = exp( ) . In this case, the diffusion enhancement (4) is equal to one. In this work, we focus on degenerate semiconductors, i.e. nonlinear diffusive problems. To compare our flux approximations, we choose the Blakemore statistics function F( ) = (exp(− ) + ) −1 with = 0.27 which is a valid approximation of the Fermi-Dirac integral of order one-half in the low density limit < 1.3 with a relative error of ≤ 0.03 (Blakemore 1952(Blakemore , 1982. Additionally, an expensive but accurate numerical flux is known for the Blakemore case (Koprucki and Gärtner 2013). The different statistics with the corresponding diffusion enhancements are depicted in Fig. 1. For brevity, we consider only the current density for electrons from now on and will partially omit the index n.

Scharfetter-Gummel type fluxes
The open, bounded domain , on which the model (1) is defined, is partitioned into N control volumes K such that = ⋃ N K=1 K , where each K is associated with a collocation point K ∈ K . We are interested in a numerical flux j along the edge connecting the collocation points of two neighboring control volumes K and L which is assumed to be aligned with the normal direction with respect to the interface K ∩ L . In the following, a subindex K corresponds to an evaluation of a physical quantity at node K and a subindex L to an evaluation at node L , respectively. Integrating (1b) over K , using the Gauss divergence theorem and one point quadrature rules yields the discrete counterpart where N( K ) denotes the set of all control volumes neighboring K . The nonlinear flux function j n;KL = j n;KL ( K , L , K , L ) approximates the projected flux ⋅ KL locally along the edge K ∩ L , where KL is the corresponding normal vector along K ∩ L . For details concerning the finite volume method see Farrell et al. (2017a).
Furthermore, 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. A numerical flux j = j KL is now said to be thermodynamically consistent, if it satisfies an analogous discrete relation, i.e.
where KL = L − K and KL = ( L − K )∕U T . Thermodynamic consistency is also important, when coupling the van Roosbroeck system to heat transport models (Farrell et al. 2017b). We discuss now different thermodynamically consistent numerical fluxes that may be used within a Voronoi finite volume framework.

Generalized Scharfetter-Gummel scheme
Under the assumption that the flux n and the electric field −∇ are constant along each face of a Voronoi cell, the flux can be projected onto the shared edge between two neighboring control volumes. Then, an integral equation can be derived (Eymard et al. 2006), which shall be satisfied by the unknown local numerical flux j where is defined in (3). The integration limits are given by K = n K , K and L = n L , L and h KL denotes the Euclidean distance between two neighboring nodes K and L . The existence of a solution to (6) was proven by Gärtner (2015), even though the integral equation is in general not explicitly solvable. We refer to the solution of (6) as generalized Scharfetter-Gummel flux. Note that for non-degenerate semiconductor devices the generalized scheme reduces to the classical Scharfetter-Gummel scheme (Scharfetter and Gummel 1969) for a non-dimensionalized edge current j sg = j∕j 0 with B as the Bernoulli function which is defined by B(x) ∶= x∕(e x − 1) , B(0) = 1 . Additionally, it was shown by Koprucki and Gärtner (2013) that for degenerate semiconductors based on Blakemore statistics the integral Eq. (6) can be reduced to a fixed point equation, namely The implicit Eq. (8) can be solved within a few Newton steps, but the efficiency of this flux is highly dependent on the choice of initial value. Hence, computationally less expensive flux discretization schemes are needed as an alternative. Still, we will use this scheme as a reference flux for the case of degenerated semiconductors, modeled by Blakemore statistics.

Diffusion enhanced scheme
Recently, another modified Scharfetter-Gummel discretization scheme was introduced (Bessemoulin-Chatard 2012). There, a logarithmic average for the nonlinear diffusion enhancement g in (4) is considered, resulting in the local flux approximation We stress that, in case of a denominator in (9) near zero, i.e. K ≈ L , regularization strategies need to be developed to handle the removable singularity.

"Sedan" scheme
The earliest reference we could find for the excess chemical potential scheme is the source code of the SEDAN III simulator (Yu and Dutton 1988), therefore in the following, we will likewise refer to this scheme as the Sedan scheme. There are benchmarks computed by the device simulator SEDAN III itself available in literature, but to the best of our knowledge there are barely any comparisons of this flux discretization scheme with other schemes known. A numerical analysis focused comparison of this flux approximation is given by Cancès et al. (2021). The scheme is motivated by rearranging the drift part in (5) to include the excess chemical potential, ex = log F( ) − , yielding with (7)

Comparison of flux discretizations
This paper aims to extend a previous discussion (Farrell et al. 2017b) by examining similar aspects for the excess chemical potential flux approximation introduced in Sect. 3.3.

Taylor expansions
Taylor expansions of the following form dependent on KL and KL can be derived for the flux approximations introduced in Sect. 3 for a sufficiently smooth statistics F , when expanding in ̄K L = ( L + K )∕2 , see Abdel (2020) and Farrell et al. (2017b). Here, the fluxes j g , j s and j d correspond to one of the flux discretization schemes introduced in Sect. 3. For the prefactors j , j ∈ {1, 2, 3} , we have The absolute error in these prefactors j between the Taylor expansions of the generalized Scharfetter-Gummel scheme and the two modified ones is depicted in Fig. 2. For the depicted figures we choose the statistics function F as the Blakemore statistics. For large negative arguments of the function F the Boltzmann and the Blakemore statistics nearly coincide, corresponding to the non-degenerate case. Hence, the modified Scharfetter-Gummel schemes converge towards the classical scheme (7) and we observe nearly vanishing errors in Fig. 2. For large positive arguments we observe that the errors in the prefactors corresponding to the diffusion enhanced scheme increase with the exponents of KL , whereas the error corresponding to the excess chemical potential scheme nearly vanishes. Due to this observation it gives rise to think that in case of no electrical field the excess chemical potential flux performs better than the diffusion enhanced scheme. However, neither the diffusion enhanced nor the excess chemical potential scheme is third-order accurate. To measure the quality of the fluxes in a different manner, second-order error estimates for the local flux errors are considered next.
When neglecting third-order terms, the following error bounds between the modified and the generalized flux dependent on the diffusion enhancement can be derived (Abdel 2020;Farrell et al. 2017b) The error bounds indicate a better performance of the diffusion enhanced scheme for large values of the diffusion enhancement g, i.e. for statistics strongly deviating from the Boltzmann regime.

Error between local flux approximations
We study the logarithmic error between the modified flux schemes and the generalized scheme for two fixed averages ̄K L . The errors for the simulation of a degenerate semiconductor can be seen in Fig. 3. The black dashed lines correspond to thermodynamic consistency, as well as pure drift currents, i.e. K = L . In both cases, the modified schemes agree exactly with the generalized scheme.
Since F(̄K L )∕g(̄K L ) = F � (̄K L ) , the derivative of the statistics function appears in the error estimates (14) and (15). The derivative of the Blakemore statistics decreases for large positive arguments. Hence, we observe in Fig. 3 that increasing the average ̄K L results in a comparatively smaller error. Both, the error estimates (14), (15) and Fig. 3 indicate a larger area, where the diffusion enhanced and the generalized scheme agree well for small values of KL and large values of the potential difference KL . Further, the red dashed line in Fig. 3 indicates agreement of the excess chemical potential scheme and the exact solution of (6) for a purely diffusive flux , i.e. a vanishing electrical field KL = 0 . This can be proven analytically, see Abdel (2020). In this specific case, the excess chemical potential scheme is the best possible flux approximation.

Numerical example
Finally, we study the impact of the different flux approximations on the simulation of degenerate semiconductor devices for a 6 m long GaAs p-i-n diode with a width of 0.5 m and a depth of 1.0 ⋅ 10 −4 cm . A motivation for this example is given in Farrell et al. (2017b). Since the 3D device only varies along one axis, it suffices to simulate the device along one spatial dimension. On each 2 m long layer (p/n−doped or intrinsic) we choose N = 3 ⋅ 2 n ref −1 uniform nodes. The donor and acceptor doping densities are given by N D = 4.35 × 10 17 cm −3 and N A = 4.20 × 10 18 cm −3 . An open source solver, based on VoronoiFVM (Fuhrmann 2019(Fuhrmann -2020, was used which allows to use automatic differentiation. The stationary van Roosbroeck system (1) with zero recombination supplemented with Dirichlet-Neumann boundary conditions is considered. The resulting current voltage curves for a refinement level n ref = 3 and the L ∞ errors in the computed total currents based on the different flux approximations for the first nine refinement level are depicted in Fig. 4. It can be observed that eventually the errors in the computed total currents based on the flux schemes converge with order O(h 2 ) . Furthermore, it suggests that on coarse meshes, which are hard to avoid for expensive 3D simulations, the excess chemical potential flux performs better than the diffusion enhanced scheme, yet multidimensional device set-ups were not simulated.

Conclusion
Our goal was to assess the quality of the excess chemical potential flux (11) which has received surprisingly little attention in the literature. To this end, we compared it to another modified Scharfetter-Gummel scheme (10) by studying its error with respect to the more accurate but expensive integral flux (6). For this, we analyzed Taylor expansions of the flux discretization schemes, the errors in the local flux approximations and simulated a p-i-n benchmark. Further applications and the impact of this scheme on realistic multidimensional device settings will be part of future research.