2 The SHE Equations

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 [23, 110]. 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. [37, 38] 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 [107, 106]. 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 [108, 109], 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 [82, 81, 83] 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. [26, 35] 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. [69, 45] 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. [75, 74, 76] 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 [18, 41, 20], 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 Ω Ω

(2.2)

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

|k|2∂|k| 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) B

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.

Term $∂f∕∂t$ : Since the time derivative can be pulled in front of the projection integral, one immediately obtains $∂[f ]l,m ∕∂t$ , where $∫ [f ]l,m := Yl,mf Zd Ω (2.7) Ω$ is in general different from the expansion coefficient $fl,m$ as outlined above.
Term $v ⋅∇xf$ : Similar to the previous term, the gradient with respect to the spatial coordinate $x$ can be pulled in front of the integral, thus leading to $∫ ∇x ⋅ Yl,mvf Zd Ω . (2.8) Ω$ It is convenient to associate the integral with the generalized current density $∫ jl,m (x, ε,t) = Yl,mvf Zd Ω . (2.9) Ω$

Term $F ∕ℏ⋅∇kf$ : In contrast to the other terms, the derivative cannot be pulled out of the projection integral. Since the gradient with respect to the wavevector $k$ of the distribution function $f$ would lead to problems for a separation of $fl′,m ′$ and $l′,m′ Y$ after inserting the expansion (2.1), an integration by parts is carried out. To avoid formal difficulties with the delta distribution, the projection is in addition first multiplied with a smooth test function $ψ (ε)$ and integrated over energy:

∫ ₀^∞ψ(ε)	∫ _Bδ(ε − ε(k))Y ^l,m ⋅∇_kfdkdε
	= ⋅∫ _Bψ(ε(k))Y ^l,m∇_kfdk
	= − ⋅∫ _B∇_k(ψ(ε(k))Y ^l,m)fdk
	= − ⋅∫ _Bfdk

The gradient in spherical coordinates $(|k |,θ,φ)$ in $k$ -space is given by

∂Y-l,m- 1--∂Y-l,m-- ---1---∂Y-l,m- ∇kYl,m = ∂|k| e|k| + |k | ∂θ eθ + |k|sin θ ∂φ eφ , (2.10)

where $e|k|$ , $eθ$ and $eφ$ denote the unit vectors in $|k|$ , $θ$ and $φ$ -direction, respectively. The invariance of spherical harmonics with respect to the radial direction and the relation $v = ∇k ε∕ℏ$ leads to

∫ ₀^∞ψ(ε)	∫ _Bδ(ε − ε(k))Y ^l,m ⋅∇_xfdkdε
	= −F ⋅∫ ₀^∞∫ _ΩY ^l,mvfZdΩdε
	−F ⋅∫ ₀^∞ψ(ε)∫ _ΩfZdΩdε
	= −F ⋅∫ ₀^∞j_l,m + ψΓ_l,mdε
	= F ⋅∫ ₀^∞ψdε ,

where the angular coupling term is given by

∫ ( l,m l,m ) Γ l,m = --1- ∂Y----eθ + -1--∂Y----eφ fZd Ω . (2.11) Ω ℏ|k| ∂θ sin θ ∂ φ

Since the test function $ψ$ can be taken ‘arbitrarily’, one obtains

∫ ( ∂j ) --2-- δ(ε − ε(k))Yl,m F- ⋅∇xf dk = F ⋅ --l,m- − Γ l,m . (2.12) (2π )3 B ℏ ∂ ε

Term $Q{f }$ : The scattering operator in the low-density approximation is considered, which neglects the nonlinearity introduced by Pauli’s exclusion principle: $∫ Q {f} = --1-- s(x,k ′,k)f(x, k′,t)− s(x, k,k′)f (x, k,t)dk′ . (2.13) (2π)3 B$ Note that it is common to multiply the integral with a small sample volume $Vs$ in front, which is then included in the denominator of the scattering rates. In this work, the sample volume $V s$ is not written explicitly. To allow for several different scattering processes like acoustical and optical phonon scattering, the index $η$ is used for writing the scattering rate as $′ ∑ ′ ′ s(x,k,k ) = ση(x,k, k)δ(ε(k)− ε(k )± ℏωη) , (2.14) η$ where the second and third arguments of $s$ denote the initial and the final state respectively. The minus sign stands for emission of energy and the plus sign for absorption. In the case of multiple energy bands, summation over all energy bands has to be added to (2.14), cf. [53]. In the following only a single band is considered, a generalization to multiple bands mainly consists of a summation over all energy bands involved.
The scattering integral is split into an in-scattering term
$in 1 ∫ ′ ′ ′ Q {f} = (2π-)3 s(x, k,k )f (x, k ,t)dk (2.15) B$ and an out-scattering term $∫ out -1--- ′ ′ Q {f } = (2π )3 f(x,k,t) s(x,k,k )dk , (2.16) B$ which are projected onto spherical harmonics one after another. Projection of the in-scattering term (2.15) yields
$∫ in --2-- l,m in ∑ in Q l,m {f} = (2π)3 δ(ε − ε(k ))Y Q {f}dk =: Q η,l,m {f} , B η$ (2.17)

where
$∫ ∫ Qiηn,l,m {f} = Y l,mZ ση(x,k(ε± ℏω η,θ′,φ′),k(ε,θ,φ)) Ω Ω × f(x, ε± ℏω η,t)Z (ε ± ℏωη,θ′,φ′)d Ω′dΩ .$ (2.18)

Similarly, the out-scattering term is evaluated as
$Qout{f } = ∑ Qout {f } , (2.19) l,m η η,l,m$ where
$∫ ∫ out l,m ′ ′ Q η,l,m{f } = Y f Z Z (ε ± ℏωη,θ ,φ ) Ω Ω ′ ′ ′ × ση(x,k(ε,θ,φ),k(ε ± ℏωη,θ,φ ))dΩ dΩ .$ (2.20)

A considerable simplification can be achieved if the transition rate is assumed to be velocity randomizing, i.e. the coefficient $ση$ in (2.14) depends only on the initial and final energy, but not on the angles. This allows for rewriting (2.18) as
$∫ ∫ in,VR l,m ′ ′ ′ Q η,l,m {f} = ση(x,ε ± ℏωη,ε) Y Zd Ω f(x,ε± ℏωη,t)Z (ε± ℏω η,θ,φ )dΩ Ω Ω = ση(x,ε ± ℏωη,ε)Zl,m (ε)[f-(x,-ε±-ℏω-η)]0,0-, Y0,0$ (2.21)

where (2.7) was used and $Zl,m (ε)$ denotes the orthonormal projection of the generalized density of states onto the spherical harmonic $Yl,m$ using the standard inner product on the sphere. In the case of spherical bands, $[f] = f Z 0,0 0,0$ and thus only $f 0,0$ is coupled into the balance equation for $l = m = 0$ .
With the assumption of velocity-randomization, the out-scattering term can be simplified to
$∫ ∫ out,VR l,m ′ Q η,l,m {f} = ση(x,ε,ε ± ℏωη) × Y fZd Ω × Z (ε± ℏω η,θ,φ)dΩ ∑ Ω Ω = ση(x,ε,ε± ℏω η)[f]l,m Z0,0(ε±-ℏ-ωη)-. η Y0,0$ (2.22)

Thus, the out-scattering term is proportional to $fl,m$ in the case of a spherically symmetric density of states $Z$ . If an expansion of the generalized distribution function $g = f Z$ is carried out instead of an expansion of $f$ , then the out-scattering term is proportional to $gl,m$ irrespective of any spherical symmetry of $Z$ .
Summing up, the full projected scattering operator using velocity randomization (VR) is thus given by
$[ QVR {f} = --1- ∑ Z (ε)σ (x, ε± ℏω ,ε)[f(x,ε ± ℏω ,t)] l,m Y 0,0 η l,m η η η 0,0 ] − [f]l,m(x,ε,t)ση(x,ε,ε∓ ℏωη)Z0,0(ε ∓ ℏωη) .$ (2.23)

Even though the scattering operator after projection is not an integral operator any longer, shifted arguments on the right hand side show up whenever inelastic collisions characterized by $ℏω η ⁄= 0$ are considered.

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(ε∓ ℏ ωη)

(2.24)

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 ∓ ℏωη) .

(2.26)

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) q

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 ∓ ℏωη) ,

(2.29)

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 ∓ ℏω η) ,

(2.30)

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 ∓ ℏωη) ,

(2.34)

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 ∓ ℏωη) ,

(2.35)

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) [53, 42].

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 η

(2.38)

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

|q| − Δ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$ .

Solve the Poisson equation (2.43) for $ψ(k+1)$ using $n(k)$ and $p(k)$ .
Solve the SHE equations (2.44) for $n(k+1)$ using $ψ (k+1)$ .
Solve the hole continuity equation (2.45) for $p(k+1)$ using $ψ(k+1 )$ .
Stop if $(ψ(k+1),n(k+1 ),p(k+1))$ is sufficiently accurate.
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$

Solve (2.46) for $δψ (k)$ using $n(k)$ and $p(k)$ .
Set $ψ(k+1) = ψ(k) + αδψ (k)$
Solve the SHE equations (2.44) for $n(k+1)$ using $ψ (k+1)$ .
Solve the hole continuity equation (2.45) for $p(k+1)$ using $ψ(k+1 )$ .
Stop if $(ψ(k+1),n(k+1 ),p(k+1))$ is sufficiently accurate.
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:

$∂n ∕∂f$ : The electron density is obtained from the distribution function as $∫ ∫ n = ∞ fZd Ωd ε , (2.54) 0 Ω$ hence the derivatives with respect to the expansion coefficients are given by $∫ ∞ ∫ -∂n--= Yl,mZd Ωd ε . (2.55) ∂fl,m 0 Ω$ In the case that $Z$ does not show an angular dependence, a nonzero derivative is obtained for $f0,0$ only.
$∂(L − Q )∕∂ψ l,m l,m$ : Here, the derivative of the SHE equations (2.34) and (2.35) with respect to the potential $ψ$ need to be taken. Since most terms depend on the kinetic energy $ε$ , which for a fixed total energy $H$ depends on the potential $ψ$ , the derivative is very elaborate. For the sake of brevity, only the derivative of the generalized current density term $jl,m$ is outlined explicitly: $∂jl,m-(ε)-= ∂jl,m ∂ε-= − q∂jl,m , (2.56) ∂ψ ∂ε ∂ψ ∂ε$ where it should be noted that $q$ refers to the signed charge of the carrier. In the same way derivatives of the terms $f l,m$ , $Γ l,m$ and $Z$ are computed. The force term is handled in the same way as for the hole continuity equation. For the case that no analytical expression for $jl,m$ , $Γ l,m$ or $Z$ is available, a discrete differential quotient is taken.

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.

Accuracy: The main rationale for the introduction of the deterministic SHE method was the ability to compute solutions of the full BTE instead of simpler macroscopic equations obtained from moments of the BTE. Moreover, the SHE method has been shown e.g. by Hennacy et al. [37] or Jungemann et al. [53] to yield results in agreement with those obtained from the Monte Carlo method, which is commonly considered to be the reference for all macroscopic models. Therefore, SHE easily fulfills the requirement for accuracy in modern TCAD.
Charge Conservation: The box integration scheme proposed by Jungemann et al. [53] and later refined by Hong et al. [42] ensures current conservation by construction.
Self-Consistency: As discussed in Sec. 2.4, the SHE equations are solved self-consistently with the Poisson equation. The convergence behavior is similar to that of the drift-diffusion model.
Resolution of Complicated Domains: Except for the use of spatially triangular grids for first-order SHE [105], publications on the SHE method have relied on structured grids so far. For one-dimensional simulations, this is clearly not a concern. However, the first two-dimensional higher-order SHE simulations have been reported just recently in 2008 by Hong et al. [44] on structured grids. Even though the discretization scheme is in principle suitable for unstructured grids, no results have been reported by then. Nevertheless, a commercial implementation of a SHE simulator using triangular grids became available from Synopsys Inc. in 2011.
Computational Resources: The SHE method imposes high memory requirements due to the additional energy coordinate. Memory consumptions in the range of 25 Gigabytes were reported for a third-order SHE simulation of a spatially two-dimensional simulation of a SiGe HBT [42], leading to single-threaded execution times of several hours. While execution times are still well below those of the Monte Carlo method, they prohibit parametric studies within a reasonable amount of time, as well as time-dependent simulations. Moreover, the parallelization of the SHE method has not been investigated so far, while fully parallel Monte Carlo simulators are already available [114]. Consequently, improvements from the numerics point of view as well as from the scientific computing point of view are required to make the SHE method attractive of modern TCAD.
Extendibility: Due to the formulation by means for partial differential equations, the SHE method can be extended in a similar way as the drift-diffusion method. For example, the inclusion of magnetic fields has been proposed by Hong et al. [42]. Except for carrier-carrier scattering, for which a scheme is proposed in Sec. 3.4, the different scattering operators used with the Monte Carlo method have also been employed for the SHE method [46]. Thus, the SHE method fully complies to the requirements of modern TCAD in terms of extendibility.

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.

[next] [prev] [prev-tail] [front] [up]

Chapter 2The SHE Equations