Chapter 2
The SHE Equations

This chapter gives a brief introduction to the spherical harmonics expansion (SHE) method for the deterministic numerical solution of the BTE. First, the developments since the introduction of the method in the early 1990s are discussed. Connections with the new contributions presented in this work are established as much as possible. Then, the formulation of the SHE equations as proposed by Jungemann et al.  [53] and later refined by Hong et al.  [42] is presented, which is also the foundation for the contributions in subsequent chapters. The chapter closes with a comparison of the requirements on modern TCAD given in Sec. 1.2 and the current state-of-the-art of the SHE method.

2.1 Historical Overview

First journal publications on the SHE method date back to Gnudi et al.  [21] as well as Goldsman et al.  [25] in 1991. However, a precursor of the method has already been used in the PhD thesis of Goldsman in 1989  [24]. Even though the method formally relies on an expansion of the distribution function f(x,k,t)  with spatial location x , wave vector k and time t  into spherical harmonics Y l,m  of the form

           ∑∞  ∑ l
f (x, k,t) =        fl,m (x,ε,t)Y l,m (θ,φ ) ,                (2.1)
           l=0 m=− l
where the three-dimensional wave vector k is written in spherical coordinates ε  , θ  and φ  on equi-energy surfaces, the derivation of equations for the expansion coefficients was rather based on perturbations of the ground state at l = 0  than on a clean mathematical approach such as a Galerkin scheme. Consequently, the derivation of the equations required a lot of bookkeeping. Nevertheless, the first-order expansion was extended to the two- and multi-dimensional case soon after  [23110]. It is worthwhile to mention the H  -transform, which has been proposed already at a very early stage in the development of the SHE method  [23].

In the mid-1990s a large number of publications dealt with the SHE method. Lin et al. derived a Scharfetter-Gummel-type stabilization  [65] for the first-order SHE method, and coupled the SHE method for electrons with the Poisson equation and the hole continuity equation  [66]. Hennacy et al.  [3738] extended the method to arbitrary order expansions. A comparison of different boundary conditions has been carried out by Schroeder et al.  [92]. Vecchi et al.  [105] proposed an efficient solution scheme based on a multigrid-like refinement of the grid near the conduction band edge. Moreover, a decoupling of the system of linear equations after discretization is discussed for the case that only one inelastic scattering mechanism with constant energy transfer is considered. Even though the techniques presented therein differ substantially from those discussed in Chap. 6 and Chap. 7, the crucial role of inelastic scattering for the coupling structure and methods for the resolution enhancement near the band edge have been discussed already. The same group proposed a scheme for incorporating full-band effects for both the conduction and the valence band  [107106]. The publication by Rahmat et al.  [78] observed that the spatial terms of the SHE equations can be handled separately from the angular coupling, which is the key observation used throughout Chap. 4. However, this separation was given as an implementation hint for the assembly of the Jacobian matrix only, while in Chap. 4 the system matrix is kept in a compressed form and is never set up explicitly. The same publication also covers aspects of numerical stability and proposes an upwinding scheme based on trial and error. The attractiveness of the SHE method for the investigation of hot-carrier effects is reflected in publications on the inclusion of electron-electron scattering for first-order SHE  [108109], and by investigations of impact ionization  [22] as well as hot-electron injections  [77]. Singh  [97] discussed the advantage of using inelastic scattering in order to avoid spurious oscillations in the distribution function and proved that the resulting system matrix for first-order SHE is an M  -matrix, which ensures the positivity of the distribution function at the discrete level.

The number of publications on the SHE method by the engineering community declined towards the end of the 20th century. Instead, publications from mathematicians increased. Ben Abdallah  [2] put the first-order SHE method in context with established moment-based methods and derived high-field approximations  [3]. Ringhofer proposed various expansion approaches for the BTE based on entropy functionals  [828183] and proposed an expansion of the energy coordinate into Hermite polynomials  [80]. Hansen et al.  [36] analyzed the SHE method for plasma physics.

In the early 2000s, Goldsman et al.  [2635] applied a first-order SHE method to the modified Boltzmann equation taking contributions from the Wigner equation into account. The publication of Jungemann et al.  [53] in 2006 applied the maximum entropy dissipation scheme developed by Ringhofer to a discretization in (x,ε)  -space, where ε  denotes kinetic energy, and put arbitrary order expansions on solid grounds by using a Galerkin scheme for the angular components. Furthermore, higher-order expansions were shown to be required for devices in the nanometer regime, a box-integration scheme suitable also for a certain class of unstructured grids was proposed, and the need for good preconditioners in order to obtain convergence of iterative solvers was discussed. Hong et al.  [44] presented the first arbitrary-order implementation of the SHE method in two spatial dimensions in 2008, refined the numerical scheme in 2009 to cover magnetic fields, and re-introduced the H  -transform in order to preserve numerical stability in the deca-nanometer regime  [42]. A full-band SHE of the valence band was presented by Pham et al.  [73]. Recently, Matz et al.  [6945] presented a fitted band structure for the conduction band for use with the SHE method, which further increases the accuracy of the SHE method, but at the expense of considerably increased computational costs. In 2011, Jin et al.  [48] proposed a simplified method for the inclusion of full-band effects, which preserves the accuracy of the fitted band structure without notably increasing the computational effort. Hong et al.  [43] also extended the SHE method to obey Pauli’s exclusion principle, while Pham et al.  [757476] coupled the Schrödinger equation with a SHE method of reduced dimensionality for multi-subband solutions of PMOSFETs.

2.2 Derivation of the SHE Equations

The SHE equations are derived in the following. The calculations mostly follow those given by Jungemann et al.  [53]. For reasons of clarity, function arguments are omitted whenever appropriate. Furthermore, the explicit inclusion of the Herring-Vogt transform as well as the inclusion of a magnetic field, which are detailed by Hong et al.  [42], are avoided. For details on spherical harmonics the reader is referred to the literature  [184120], and for the particular phase factors used within this thesis to  [86].

The derivation is carried out for an expansion of the distribution function f (x,k,t)  into spherical harmonics as in (2.1) on equi-energy surfaces, which requires that the mapping (ε,θ,φ) ↦→ k is a bijection. For a given distribution function f  , the expansion coefficients are obtained from projections:

             --2--     l,m                                3
fl,m (x,ε,t) = (2π)3 BY    (θ(k),φ (k ))f (x, k,t)δ(ε − ε(k))dk
           =    Y l,m (θ,φ)f(x,k (ε,θ,φ ),t)Z (ε,θ,φ )d Ω

Here, Ω  denotes the unit sphere, and the generalized density of states including spin degeneracy in the absence of magnetic fields is given by

Z(ε,θ,φ) = 4π3  ∂ε  ,                          (2.3)
where a possible spatial dependence due to the use of different materials is ignored for the sake of conciseness. One should be careful when comparing generalized densities of states from different authors, because additional prefactors are in use.

It is important to keep in mind that the spherical harmonics are not orthogonal with respect to the bilinear form

(f,g) :=    fgZd Ω ,                           (2.4)
     e    Ω
unless the generalized density of states Z  does not depend on the angles. For this reason, one may alternatively expand the generalized distribution function g := fZ  and use the standard inner product on the unit sphere, which is the approach taken by Jungemann et al.  [53]. In the following, however, an expansion of f  is carried out, but differences to an expansion with respect to g  are pointed out on a regular basis.

One pitfall in the derivation of the SHE equations for an expansion of f  is related to the representation (2.2) for a given distribution function f  . The BTE is to be solved for an unknown distribution function f  , hence f
 l,m  is not obtained directly from a given function. Instead, the unknown distribution function f  is expanded into spherical harmonics, which ultimately results in equations for the expansion coefficients fl,m  . Keeping in mind that the spherical harmonics are not necessarily orthogonal with respect to the bilinear form (2.4), there holds

∫              ∫                                    ∫
     l,m             l,m ∑         l′,m ′      ∑            l,m  l′,m′
 Ω Y   fZd Ω ≃   ΩY    ′  ′fl′,m′Y    Zd Ω =  ′ ′fl′,m ′ Ω Y   Y    Zd Ω ,    (2.5)
                       l,m                  l,m
which results in the expansion coefficient fl,m  only in the case of a generalized density of states independent of the angles, i.e. Z = Z (ε)  , but not for the general case of an angular dependency.

In order to derive a set of equations for the expansion coefficients fl,m  for an unknown function f  known to fulfill the BTE (1.1), one can proceed in two ways. The first possibility is to insert the expansion into the BTE and then project the resulting equation onto spherical harmonics, while the second possibility is to first project the BTE and then insert the expansion of f  . Both operations are linear and the BTE itself is linear when using linear scattering operators only, thus both methods yield the same result provided that all required exchanges of summation and integration are valid. Since the latter method leads to less notational clutter, first a projection of the BTE onto each of the spherical harmonics Y l,m  of the form

       2  ∫      l,m             3
X ↦→  (2π)3   XY    δ(ε− ε(k))dk                      (2.6)
is carried out. These projections onto Y l,m  are detailed in the following term-by-term. Function arguments are usually suppressed to increase the legibility of the expressions.

Collecting all individual terms, we obtain the system of equations for the BTE upon projection onto spherical harmonics on equi-energy surfaces under the assumption of velocity randomization:

∂[f]                         (∂j          )
---l,m--+ ∇x ⋅jl,m (x, ε,t) + F ⋅ ---l,m − Γ l,m   =
  ∂t         ∑   [              ∂ε
       = -1--    Zl,m(ε)ση(x,ε ± ℏωη,ε)[f ]0,0(x, ε± ℏωη,t)
         Y 0,0  η
                   − [f]l,m (x,ε,t)ση(x,ε,ε∓ ℏω η)Z0,0(ε∓ ℏ ωη)

It is worthwhile to note that an expansion of g = fZ  instead of f  leads to the same system of equations when replacing [f]  with g  .

2.3 H  -Transform and MEDS

A discretization of the system (2.24) in (x, ε)  -space leads to spurious oscillations for very small devices, because the grid is not aligned with the trajectories of carriers in free flight. The numerical stability can be improved substantially by a transformation from kinetic energy ε  to total energy H  , because trajectories are then given by H = const  . The price to pay for the increased numerical stability is the additional effort for handling the band-edge given by ε = 0  , because the simulation domain needs to be adjusted with every update of the electrostatic potential.


Figure 2.1: Trajectories of carriers in free flight within the device are given by constant total energy H  .

The H  -transform was suggested already soon after the SHE method had been introduced  [23]. In the following, the application of the H  -transform to (2.24) in the way proposed by  [42] is shown. Consider a variable transformation from (x, ε)  to (˜x,H )  by

˜x = x ,                    H =  ε+ qΨ (x) ,              (2.25)
where Ψ (x)  can be an arbitrary function of x and q  is the signed charge of the carrier. The system (2.24) is then transformed to
∂[f]l,m-+ ∇ ˜x ⋅jl,m + [F + ∇xH ]⋅ ∂jl,m − F ⋅Γ l,m
  ∂t            [               ∂H
       = --1-∑   Z   σ (˜x, H ± ℏω ,H )[f]  (x˜, H ± ℏω  ,t)
         Y0,0 η    l,m η          η      0,0          η
                   − [f]l,mσ η(˜x, H,H  ∓ ℏωη)Z0,0(H  ∓ ℏωη)  .

Since the force term is given by

F  = − ∇x [Ec (x )+ qψ (x )] ,                      (2.27)
where Ec   is the position-dependent valley minimum and ψ  the electrostatic potential, the derivative of jl,m  with respect to H  is eliminated by the choice
Ψ (x) = ψ(x) + Ec(x)-.                         (2.28)
This particular choice of Ψ  implies that H  refers to the total energy. The H  -transformed system is thus given by
∂[f]l,m-+ ∇  ⋅j   − F  ⋅Γ   =  -1--∑   Z  σ  (x, H ± ℏω  ,H )[f]  (x,H ± ℏω  ,t)
  ∂t       x  l,m        l,m    Y0,0      l,m  η          η      0,0          η
                                   η                                      ]
                                     −  [f]l,mσ η(x, H,H  ∓ ℏωη)Z0,0(H  ∓ ℏωη)  ,

where x is written instead of ˜x .

Inserting the expansion (2.1) into (2.29) yields

   l′,m ′  ∂fl′,m ′       l′,m′            l′,m ′
[Y    ]l,m--∂t-- + ∇x ⋅jl,m  fl′,m ′ − F ⋅Γ l,m fl′,m′
                    1 ∑  [                        ′ ′
                = ----    Zl,mση(x, H ± ℏωη,H )[Yl,m ]0,0fl′,m′(x,H  ± ℏωη,t)
                  Y0,0 η
                               l′,m ′                                    ]
                          − [Y    ]l,mfl′,m′ση(x,H, H ∓ ℏω η)Z0,0(H ∓ ℏω η) ,

where Einstein’s summation convention over pairs of upper and lower indices is employed. The coefficients

   l′,m ′          l,m  l′,m′
[Y    ]l,m :=    Y    Y    ZdΩ  ,                                 (2.31)
      ′  ′   ∫Ω        ′ ′
     jll,,mm  :=    Y l,mvY  l,m Zd Ω ,                                (2.32)
             ∫Ω     (                       )
      l′,m ′      -1--  ∂Yl,m--    -1--∂Y-l,m-      l′,m ′
    Γ l,m  :=   Ω ℏ|k |   ∂θ  eθ + sin θ  ∂φ  eφ  Y    Zd Ω         (2.33)
determine the coupling between the different equations and will be thoroughly investigated in Chap. 4.

The even part of the carrier distribution function can be associated with densities like the carrier density or the energy density, while the odd part of the carrier distribution function is viewed as fluxes. To improve numerical stability, these fluxes need to be stabilized  [82]. Based on entropy dissipation principles, Ringhofer proposed to multiply the projected SHE equations with the entropy function exp(H ∕(kBTL))  , where kB   is the Boltzmann constant and TL   denotes lattice temperature, and to take the negative adjoint form of the projected equations for odd l   [82]. While the multiplication with the entropy function is crucial for a formulation based on kinetic energy as in (2.24) and  [53], it is just a constant factor in a formulation based on total energy such as (2.30). Therefore, it is sufficient to simply take the negative adjoint operator  [42]:

          ′ ′   ∂fl′,m′       (  ′ ′    )        ′ ′
l even : [Yl,m ]l,m----- + ∇x ⋅ jll,,mm fl′,m′  − F ⋅Γ ll,,mm fl′,m ′
                  ∂t         ∑  [
                       = -1--    Zl,m ση(x,H  ± ℏωη,H )[Y l′,m ′]0,0fl′,m′(x,H ± ℏω η,t)
                         Y0,0 η
                                     l′,m′                                     ]
                                 − [Y    ]l,mfl′,m ′σ η(x, H,H  ∓ ℏωη)Z0,0(H  ∓ ℏωη)  ,

l odd : [Y l′,m ′]l,m ∂fl′,m-′+ jl′,m′⋅∇xfl ′,m′ + F ⋅ ˆΓ l′,m′fl′,m ′
                  ∂t     l,m     [            l,m
                         -1--∑                          l′,m′    ′  ′
                      =  Y0,0     Zl,m ση(x,H ± ℏω η,H )[Y    ]0,0fl,m (x,H  ± ℏωη,t)
                              η                                               ]
                                 − [Y l′,m ′]l,mfl′,m′ση(x,H, H ∓ ℏ ωη)Z0,0(H ∓  ℏωη) ,

where Einstein’s summation convention over pairs of upper and lower indices is employed. The coupling coefficients [Yl′,m′]l,m  and   ′ ′
jll,,mm  are self-adjoint and thus unchanged, while  l′,m′
ˆΓl,m  is given by

  ′ ′           ∫      (    l′,m′            l′,m′  )
ˆΓ ll,,mm := Γ l,m′ ′ =  -1--  ∂Y----eθ + --1- ∂Y----eφ  Y l,mZd Ω  .       (2.36)
         l,m     Ω ℏ|k |    ∂θ       sinθ  ∂ φ
This stabilization scheme is commonly referred to as maximum entropy dissipation scheme (MEDS)  [5342].

It is worthwhile to compare (2.34) and (2.35) with the system obtained from an expansion of

                          ∑∞   ∑l              l,m
g(x,k, t) = f(x,k, t)Z (k) =         gl,m(x,H, t)Y    (θ,φ ) ,       (2.37)
                           l=0 m=− l
which is obtained for even l  after repeating the previous steps as
∂gl,m-      ˜ l′,m′           ˜l′,m ′
 ∂t  + ∇x ⋅j l,m gl′,m′ − F ⋅Γ l,m gl′,m′
        -1-∑  [
     =  Y0,0    Zl,mσ η(x, H ± ℏω η,H )g0,0(x,H  ± ℏωη,t)
             η                                     ]
               − g   σ (x,H, H ∓ ℏ ω )Z  (H ∓ ℏ ω )
                  l,m  η             η  0,0        η

with coupling coefficients

˜l′,m′       l,m   l′,m′
jl,m  =   ΩY   vY     dΩ ,                                    (2.39)
  ′ ′   ∫      (   l,m             l,m   )
˜Γ ll,,mm =    -1-- ∂Y----eθ + -1--∂Y----eφ  Y l′,m ′dΩ  .          (2.40)
         Ω ℏ|k |   ∂θ       sin θ  ∂φ
The equations for odd l  are obtained in the same way by taking the negative adjoint form. If the density of states Z  does not show an angular dependence, the two approaches are equivalent, otherwise the truncated expansions will yield different results in general. To the knowledge of the author, a systematic comparison of the formulations (2.34) and (2.35) with (2.38) for nonspherical energy bands has not been carried out yet.

The system (2.34) and (2.35) consists of an infinite number of equations due to an expansion of f  of the form (2.1). The conforming Galerkin procedure for obtaining a finite set of equations is to consider a truncated expansion

           ∑L  ∑l
f(x,k,t) =         fl,m(x, H,t)Y l,m (θ,φ) ,               (2.41)
           l=0m= −l
for a reasonable maximum (finite) expansion order L  , and to consider only a finite subset of the projected equations. The typical choice l = 0,...,L  and m  = − l,...,l  for the projected equations is used in the remainder of this thesis, hence the number of unknown expansion coefficients and the number of equations agree. Otherwise, approximations in the least-squares sense for the over- or under-determined system would have to be considered. However, least-squares problems lead to additional computational effort compared to the solution of a linear system, thus no investigations have been carried out in that direction yet.

2.4 Self-Consistency with Poisson’s Equation

So far, the force term F has been considered to be a given quantity. In a homogeneous material it is linked to the electrostatic potential ψ  by F  = qE =  − q∇x ψ  , where q  denotes the signed charge of the particle and E the electric field. With Gauss’ Law, a scalar permittivity ϵ  and the charge density |q|(p− n + C )  , where p  and n  denote the density of holes and electrons respectively, and C  refers to the density of fixed charges in the material, the Poisson equation

− Δx ψ = -ϵ (p − n + C)                        (2.42)
is obtained.

Since the carrier density depends on the solution of the BTE, which in turn depends on the solution of the Poisson equation, a self-consistent solution of both equations has to be found. Even though both equations are linear in their unknowns, the coupling is nonlinear due to the inner product of the force F with the gradient ∇kf  of the distribution function with respect to the wave vector k . Consequently, a nonlinear iteration scheme has to be employed for the solution of the coupled system

                  − Δx ψ = |q|(p − n + C) ,                           (2.43)
                 Ll,m {f} = Ql,m {f} ,  l = 0,...L, m = − l,...l ,       (2.44)
∇x ⋅(− Dp ∇xp + pμp∇x ψ) = 0 ,                                       (2.45)
where electrons are treated by the SHE method, and holes are considered by a continuity equation without recombination. For the simulation of devices based on unipolar operation, the hole continuity equation is often neglected. A similar system is obtained if the SHE method is employed for holes and a continuity equation is used for electrons. Note that in the cases where a unipolar simulation is sufficient, the continuity equation is often neglected.

2.4.1 Gummel Iteration

Given the iterates   (k)
ψ   ,  (k)
n   and  (k)
p   for the potential, electron and hole density respectively, updates can be obtained by solving the equations (2.43), (2.44), and (2.45) sequentially in a Gauss-Seidel-type manner  [32]:

Algorithm 1 (Simplest Form of Gummel Iteration). Input: Initial guesses  (0)
ψ   ,  (0)
n   and  (0)
p   , k = 0  .

  1. Solve the Poisson equation (2.43) for ψ(k+1)   using n(k)   and p(k)   .
  2. Solve the SHE equations (2.44) for n(k+1)   using ψ (k+1)   .
  3. Solve the hole continuity equation (2.45) for p(k+1)   using ψ(k+1 )   .
  4. Stop if (ψ(k+1),n(k+1 ),p(k+1))  is sufficiently accurate.
  5. Set k ← k + 1  and repeat.

In its simplest form, the Gummel iteration diverges in most cases due to the exponential dependence of the carrier concentrations on the potential.

The convergence behavior can be improved substantially by considering the dependencies n = n (ψ ) = ni exp ((ψ − φn )∕VT )  and p = p(ψ) = niexp(− (ψ − φp )∕VT )  , where φn   and φp   denote the quasi-Fermi potentials of electrons and holes respectively and VT =  |q|∕(kBT)  is the thermal voltage. A linearization leads to the damped Poisson equation for δψ (k) = ψ(k+1) − ψ(k)   :

− Δ δψ(k) + n-+-pδψ (k) = Δψ (k) + |q|(p − n + C).            (2.46)
            VT                   ϵ
Consequently, the modified Gummel iteration is obtained as follows:

Algorithm 2 (Modified Gummel Iteration). Input: Initial guesses ψ(0)   , n (0)   and p(0)   , k = 0  , α > 0

  1. Solve (2.46) for δψ (k)   using n(k)   and p(k)   .
  2. Set ψ(k+1) = ψ(k) + αδψ (k)
  3. Solve the SHE equations (2.44) for n(k+1)   using ψ (k+1)   .
  4. Solve the hole continuity equation (2.45) for p(k+1)   using ψ(k+1 )   .
  5. Stop if (ψ(k+1),n(k+1 ),p(k+1))  is sufficiently accurate.
  6. Set k ← k + 1  and repeat.

The damping parameter α  allows for additional control of the damping. Typical values are about 0.5  and may be chosen differently for each iteration. In order to ensure a reduction of the residual, an additional control of α  can be employed, where e.g. subsequent trials       −i
αi = 2  for i = 1,2,...  are chosen until the residual in the current step is reduced  [16].

A typical convergence plot for an undamped modified Gummel iteration (α = 1  ) is depicted in Fig. 2.2. Typical features are the almost constant potential update during the first iterations, and the nonuniform reduction of the potential correction in later iteration steps. These features are more pronounced for more complicated devices under higher bias.


Figure 2.2: Comparison of the typical convergence of an undamped modified Gummel scheme and an undamped Newton scheme for the simulation of a one-dimensional   +  +
n  nn   -diode.

2.4.2 Newton’s Method

The standard method for the solution of nonlinear systems of equation is Newton’s method. The big advantage over many other methods is the quadratic convergence sufficiently close to the true solution. For a system of equations of the form

        (  g1(x1,x2,...,xN)  )
        |  g (x ,x ,...,x )  |
G (x ) = ||  2  1  2.     N   || = 0 ,                  (2.47)
        (          ..         )
           gN(x1,x2,...,xN )
the update of a guess  (k)
x   to the true solution is determined by
x(k+1) = x (k) − J −G1G (x(k)) .                    (2.48)
The Jacobian matrix JG  of G is given by
      (                     )
         ∂g1   ∂g1-  ...  ∂g1-
      |  ∂∂xg12   ∂∂x2g2-      ∂∂xgN2-|
J G = ||  ∂x1   ∂x2   ...  ∂xN ||  .                    (2.49)
      |(   ...     ...   ...   ...  |)
         ∂gN- ∂gN-  ...  ∂gN-
         ∂x1   ∂x2       ∂xN
Consequently, every Newton step requires the solution of a system of linear equations and the evaluation of the residual G (x(k))  . The method fails if the Jacobian matrix does not have full rank, which is, however, unlikely in practice due to the regularizing effect of round-off errors.

A Newton method for the system (2.43)-(2.45) can be derived by writing the system in the form

  − Δ  ψ − |q|(p − n + C) = 0 =: g (ψ,f ,p) ,                            (2.50)
      x     ϵ                   1
       Ll,m {f} − Ql,m {f} = 0 =: g2 (ψ, f,p), l = 0,...L,m = − l,...l ,    (2.51)

∇x ⋅(− Dp∇xp  + pμp∇x ψ) = 0 =: g3(ψ,f ,p) ,                            (2.52)
with                            T
f = (f0,0,f1,−1,f1,0,...,fL,L)   for some maximum expansion order L  . The Jacobian matrix is then obtained as
(                                                                 )
         − Δx        (1∕VT + |q|∕ε)∂n ∕∂f        1∕VT − |q|∕ε
( ∂ (Ll,m − Ql,m )∕∂ψ    Ll,m {⋅} − Ql,m{⋅}              0            )  .   (2.53)
       ∇x ⋅p∇x                0           ∇x  ⋅(− Dp ∇x + μp ∇x ψ)
Here, two terms need additional considerations:

In practice, the Newton scheme will be used together with a damping scheme similar to the one discussed for the Gummel iteration. The existence of a damping parameter αi  leading to a reduction of the overall residual norm is ensured under reasonable assumptions about the system of equations.

An examplary convergence plot for the Newton method is given in Fig. 2.2. The quadratic convergence is readily visible and leads to a much smaller number of iterations compared to the modified Gummel method.

2.5 Compatibility with Modern TCAD

Now as the SHE equations are derived, the state-of-the-art for the SHE method – excluding the author’s own contributions presented in the remainder of this thesis – is compared to the requirements for modern TCAD established in Sec. 1.2.

The issues raised in terms of the resolution of complicated domains as well as in the use of computational resources will be addressed in the remainder of this thesis. Chap. 5 discusses the use of unstructured triangular and tetrahedral meshes for the SHE method. Chap. 6 addresses the quadratically increasing computational costs of the SHE method with the expansion order L  . In Chap. 7 a parallel preconditioner is developed, which enables the efficient use of modern many- and multi-core computing architectures.