6 Adaptive Variable-Order SHE

Chapter 6
Adaptive Variable-Order SHE

The structural properties discussed in Chap. 4 show that in the case of spherical energy bands the coupling between the SHE equations is sparse. This allows for a reduction of memory requirements for the system matrix of the discretized equations from $O (N L4)$ to $O (N L2)$ . The proposed matrix compression scheme further reduces the memory requirements for the system matrix to $O(N + L2)$ , which finally results in the total memory requirements being dominated by the $N (L+ 1 )2 = O (NL2 )$ unknowns already at an expansion order of five. For a three-dimensional device simulation using $100$ discrete total energies and $10000$ grid points ( $6 N = 10$ ), the number of unknowns for a first-order expansion are up to $6 4× 10$ , where $6 10$ even-order expansion coefficients enter the linear solver. On the same grid, a ninth-order expansion leads to a total of $8 10$ unknowns, of which $7 4.5× 10$ enter the linear solver. Even though the matrix compression scheme ensures that the total memory requirements stay within a reasonable amount of a few gigabytes, the high computational effort due to the large number of unknowns leads to long simulation times.

With the use of unstructured grids for the SHE method as described in Chap. 5, the number of grid points in the $(x,H )$ -space can be reduced to lower numbers than for structured grids without sacrificing accuracy in critical device regions. Therefore, the total number of unknowns in the linear system is reduced from $N (L + 1)2$ to $N ′(L+ 1 )2$ , where $N$ and $N ′$ refer to the number of unknowns in $(x,H )$ -space and $N ′$ can be a factor two to five smaller than $N$ , cf. Sec. 5.4.

Similar to unstructured grids, which allow for a high resolution in critical regions and a lower resolution in less important regions, a higher expansion order can be chosen at locations and energies with high influence on a target quantity such as current, carrier density, or average carrier velocity. Consequently, instead of a fixed-order expansion as in (2.41), variable-order expansions of the form

L(x∑i,Hn )∑ l f (xi,k (Hn, θ,φ ),t) = fl,m(xi,Hn, t)Yl,m (θ,φ) (6.1) l=0 m=− l

are considered in this chapter. The aim is to further reduce the total number of unknowns to $N ′(L ′ + 1)2$ , where the average expansion order $L ′$ for a grid point $xi$ obtained by averaging all positive expansion orders over all positive discrete kinetic energies can be well below an equivalent uniform expansion order $L$ without sacrificing accuracy.

6.1 Variable-Order SHE

The advantage of the SHE method is that the equilibrium distribution is described exactly by a zeroth-order expansion. This is in contrast to macroscopic transport models derived from moments of the BTE, where higher-order moments do not vanish. Therefore, it is reasonable to expect that for the resolution of small deviations from the equilibrium distribution, only a low expansion order is required for the SHE method in order to obtain good accuracy.

In this section the discretization of the SHE equations using variable-order expansions is presented. For a maximum expansion order $Lmax$ in the simulation domain, a variable-order expansion can be obtained by assembling a system matrix for maximum expansion order $Lmax$ and then imposing homogeneous Dirichlet conditions on all expansion coefficients $fl,m(xi,Hn )$ with $l > L(xi,Hn )$ . However, such an approach is impractical due to the unnecessary effort of setting up a much larger system matrix than actually required. In addition, the benefit of reduced memory required for variable expansion orders is eliminated when setting up a system matrix for uniform expansion order $L max$ first.

Figure 6.1: System matrix structure for uniform expansion orders with spherical energy bands at a constant energy. If two boxes $Bi$ and $Bj$ do not have a common interface, the respective blocks are zero. If boxes $Bi$ and $˜ Bj$ do not overlap, the respective blocks are also zero. The odd-to-odd coupling implies a diagonal block due to the structure of the scattering operators.

It is assumed that the unknowns are enumerated in such a way that the expansion coefficients for a particular box $B ∈ B$ or a dual box $B˜ ∈ ˜B$ are consecutive. In addition, the unknowns associated with boxes $B ∈ B$ are enumerated first, then the unknowns associated with the dual boxes $˜ ˜ B ∈ B$ , which leads to the system matrix structure discussed in Sec. 4.3. The resulting structure of the system matrix for a uniform expansion order is depicted in Fig. 6.1. The block sizes for each block coupling two boxes are:

	Even	Odd

Even	$even even NL × N L$	$even odd N L × N L$
Odd	$odd even N L × N L$	$odd odd NL × NL$

Here, $Neven L$ and $N odd L$ denote the number of even and odd expansion coefficients up to and including order $L$ respectively. For instance, consider the block in the rows for box $Bi$ and in the columns for the dual box $B˜i$ , which is located at the interface between $Bi$ and $Bj$ . From the discrete equations for the even-order projections in steady-state (5.9) one obtains the system matrix block for a certain discrete energy $Hn$ directly from rewriting the discrete equations in matrix form:

⌊ ( ∫H+ 1,− 1 ∫H+ L,L ) Hn−n j0,0;ddH ... Hn−n j0,0;ddH ||∑3 || ∫ H+n 1,−1 ∫H+n L,L || || Ai,jni,j;d || Hn− j2,−2;ddH ... H−n j2,−2;ddH || |⌈d=1 |( ... ... ... |) ∫ H+n 1,−1 ∫ H+n L,L H −n jL−1,L− 1;ddH ... H −n jL−1,L− 1;ddH ( ∫H+n 1,−1 ∫H+n L,L ) ⌋ ( ) | ∫ H−n+ Γ0,0;ddH ... ∫ H−n+ Γ 0,0;ddH | | βi,n;1,− 1 ∑3 || H−nΓ 1,−1 dH ... Hn− Γ L,L dH || || || .. || − Vi,jFi,j;d | Hn 2,−.2;d . H n 2,−. 2;d | | | . | , d=1 |( .. .. .. |) |⌉ ( βi,n;L,L−2 ) ∫ H+n Γ 1,− 1 dH ... ∫H+n Γ L,L dH βi,n;L,L H −n L−1,L−1;d H−n L− 1,L−1;d

(6.2)

where the subscript index $d$ denotes the $d$ -th component of the vector and the expansion order $L$ is odd. The matrix block $(Bi,Bi)$ results from the scattering operator and is diagonal up to couplings induced by energy displacements $ℏω η$ due to inelastic scattering mechanisms, cf. Thm. 3. A similar structure is obtained for the even-to-even, the odd-to-even, and the odd-to-odd block.

Figure 6.2: System matrix structure for variable-order expansions with spherical energy bands at a constant energy. The coupling blocks have varying sizes.

When allowing variable-order expansions, the dimensions of the coupling blocks can differ for each box $B i$ and each dual box $B˜ i$ . Reconsidering the example of the matrix block for the box $Bi$ and the dual box $˜ Bi$ from above, and assuming an even expansion order $Li$ for $Bi$ and an odd expansion order $˜ Li$ for the dual box $˜ Bi$ , the discrete system of equations is given by

⌊ ( ∫ H+n 1,−1 ∫H+n ˜Li,L˜i ) | | ∫Hn+− j0,0;ddH ... ∫ H−n+ j0,0;d dH | ||∑3 || Hn− j12,−,−21;ddH ... Hn− j˜L2i,−,˜L2i;ddH || | Ai,jni,j;d| Hn . . H n . | |⌈d=1 |( .. .. .. |) ∫ H+n− j1,−1 dH ... ∫ H+n− j˜LiL˜i dH H n( Li−1,Li− 1;d H n Li,Li;d ) ⌋ ∫H+n− Γ 1,−1dH ... ∫ H+n− Γ ˜Li,˜LidH ( ) 3 | ∫ HHn+ 10,−,0;d1 ∫HHn+ L0˜,0,;d˜L | | βi,n;1,−1 ∑ || H n−n Γ2,− 2;ddH ... Hn−n Γ2,i−2i;ddH || || || ... || − Vi,jFi,jni,j;d|| .. .. .. || || |( β ˜ ˜ |) . d=1 ( ∫ + . . ∫ + ˜. ˜ ) ⌉ i,n;Li,Li− 2 HHn− Γ 1L,−i1,Li;ddH ... HHn− ΓLLii,L,Lii;ddH βi,n;L˜i,˜Li n n

(6.3)

This leads to the system matrix structure shown in Fig. 6.2, where the individual small blocks are now of different dimension according to the respective expansion order.

The expansion orders for the boxes $Bi$ and the dual boxes $B˜i$ should not be chosen arbitrarily. Consider the case of spherical energy bands, where only a single box $B$ carries the maximum even expansion order $Lmax$ . If all surrounding dual boxes carry an odd expansion order of at most $Lmax − 3$ , then the highest-order expansion coefficients $fL ,m max$ for the box $B$ do not couple due to Thm. 2 with expansion coefficients of lower order. Since the right hand side of the linear equation is zero except for Dirichlet boundary values, the expansion coefficients $fLmax,m$ are computed as zero. Moreover, since the size of the final linear system of equations is determined by the expansion orders of the boxes $Bi ∈ B$ only, Thm. 2 suggests that the odd expansion order of the dual box $B˜$ should be larger than the even expansion orders of the two overlapped boxes $Bi$ and $Bj$ , otherwise the linear system of equations yields a less accurate solution at the same computational effort. Summing up, one obtains:

Guideline 1. The odd expansion order of each of each dual box $˜ ˜ B ∈ B$ should be larger than the even expansion orders of the boxes $Bi$ and $Bj$ overlapped by $˜B$ .

Conversely, since odd expansion orders larger than required do not lead to higher accuracy because of the same reasoning with the roles of $Bi$ and $B˜i$ exchanged, the odd expansion order of a dual box is set to the maximum even expansion order of the overlapped boxes plus one, i.e.

˜ order(Bi,j) = max (order(Bi),order(Bj))+ 1 . (6.4)

Consider two neighboring boxes $Bi$ and $Bj$ with expansion orders $Li$ and $Lj$ , $Li < Lj$ . Denote the dual box overlapping $Bi$ and $Bj$ with $B˜i,j$ . From Thm. 2 and Guideline 2 it follows that the expansion coefficients defined for $Bi$ couple with expansion coefficients in $B˜ i,j$ up to order $L + 1 i$ . These expansion coefficients in $B˜ i,j$ couple according to Thm. 2 with expansion coefficients in $Bj$ up to order $Li + 2$ , thus there is no gain in accuracy for $Lj$ larger than $Li + 2$ . This leads to the second guideline:

Guideline 2. The maximum even expansion order of neighboring boxes $Bi$ and $Bj$ should not differ by more than two.

6.2 Adaptive Control of the SHE Order

In the previous section the efficient assembly of SHE equations using variable expansion orders has been discussed. However, the variable expansion orders were assumed to be given for each box $B$ and each dual box $B˜$ . For every-day TCAD purposes, such a selection should be based on an adaptive strategy until a prescribed accuracy for a certain target quantity is reached. This section thus discusses strategies to automatically pick suitable expansion orders within the simulation domain. Due to Guideline 1, it suffices to select expansion orders for the the boxes in $B$ .

What is considered to be a suitable or good expansion order for a box $B$ clearly depends on the target quantity. Macroscopic quantities are typically obtained by moments of the distribution function such as

∫ X (x,t) = Φ(k)f (x, k,t)dk3 , (6.5) ℝ3

where $Φ (k)$ is a suitable polynomial of $k$ . Due to the asymptotically exponential decay of the distribution function $f$ with the modulus of $k$ , the distribution function needs to be computed with high accuracy at low energies, while lower accuracy is sufficient at higher energies for the computation of $X$ . However, for the investigation of high-energy phenomena such as hot carrier degradation [101, 9], the distribution function needs to be fairly accurate at high energies, while there is less emphasis on accuracy at lower energies. Consequently, suitable expansion orders depend on the quantities of interest.

In the following, three different strategies for the adaptive choice of expansion orders are presented. The first approach is based on an analytical result for the SHE of a function $g(θ,φ)$ in dependence of the smoothness of $g$ . The second approach increases the expansion order particularly in regions with high weight on one or more target quantities. The third approach is a rather classical residual-based technique, which is common in finite element and finite volume methods [1, 8, 56]. Even though the three strategies are presented separately, in practice they should be combined in order to obtain best results.

6.2.1 Rate of Decay of Expansion Coefficients

The SHE can be seen as the three-dimensional extension of Fourier series. Since the latter is more widely known than the former, the motivation for the adaptive strategy for SHE is first given for Fourier series. The transition from Fourier series to SHE does not impose additional complications then.

(a) Plot of $f1(φ )$

(b) Plot of expansion coefficients of $f1(φ)$

(d) Plot of expansion coefficients of $f(φ) 2$

Figure 6.3: Comparison of the Fourier series of a differentiable function $f1$ and a discontinuous function $f2$ .

A $2π$ -periodic function $f (φ )$ can be expanded into trigonometric functions as

∞∑ f (φ ) = am cos(m φ)+ bm sin(m φ ) , (6.6) m=0

where $b0 = 0$ . The function $f$ does not need to be continuous, it is sufficient for $f$ to be integrable. The expansion coefficients are computed from $f$ as

1 ∫ π a0 = --- f(φ)dφ , (6.7) 2π∫ −π 1- π am = π −π f(φ)cos(m φ)dφ , (6.8) 1 ∫ π bm = -- f(φ)sin(mφ )dφ . (6.9) π −π

Consider the two $2 π$ -periodic functions depicted in Fig. 6.3

{ 2 4 4 4 2 2 1, |φ| ≤ π∕2 , f1(φ ) = [(x− π )(x + π)] ∕π = x ∕π − 2x ∕π + 1 , f2(φ ) = 0, π ∕2 < |φ| ≤ π . (6.10)

While $f1$ is continuously differentiable, $f2$ is discontinuous. The Fourier series of the functions are given by

-8- ∞∑ (−-1)m 1- ∑∞ -2(−-1)m--- f1(φ) = 15 − 48 π4m4 cos(m φ ) , f2(φ) = 2 − π(2m − 1) cos((2m − 1)φ ) . m=1 m=1 (6.11)

For the smooth function $f1$ , the Fourier coefficients decay with the fourth power of $m$ , while the Fourier coefficients of $f2$ only decay inversely proportional to $m$ . More generally, it can be shown that for a function with integrable derivatives up to order $k$ , the expansion coefficients decay at least with rate $m −k$ . Conversely, the rate of decay of the Fourier coefficients are a measure for the smoothness of the expanded function [28, 13].

For SHE, a similar result holds for a function $f$ defined on the unit sphere $Ω$ [31, 18]:

Theorem 4. Let $f$ be $k$ -times continuously differentiable on the unit sphere $Ω$ , where $k$ is a nonnegative integer. For the spherical harmonics expansion

∞ l ∞ ∑ ∑ l,m ∑ f = fl,mY =: Ql , (6.12 ) l=0m= −l l=0

where $Ql$ is a spherical harmonic of degree $l$ , there holds

|Ql(θ,φ)| ≤ Cl1∕2− k , (6.13 )

for a constant $C = C (f ) > 0$ independent of $k$ and all angles $θ$ and $φ$ .

Since the spherically symmetric distribution function $f$ in equilibrium is described exactly by a zeroth-order expansion, higher-order expansion coefficients provide a measure for the distortion of $f$ from equilibrium. This is in contrast to macroscopic models derived from moments of the BTE such as those described in Sec. 1.1, where even moments of the distribution function in equilibrium do not vanish.

Given a numerical solution of the SHE equations for certain expansion orders $L(xi,Hi)$ , an adaptive adjustment of the expansion order can be based on (6.13). Division by $f0,0 ≈ C$ and taking the logarithm leads for $l > 2$ to

log|Ql|− logf0,0 ≲ (1∕2 − k)logl , (6.14)

where here and in the following $log |Ql|$ refers to $max θ,φ log |Ql(θ,φ )|$ . Taking the left hand side as an estimate for $k$ and dropping the constant $1∕2$ , which acts only as an offset to the estimate on the rate of decay $k$ , results in

log |Q |− log f η := − k ∼-----l-------0,0-, (6.15) log l

with $η$ taking the role of an error estimator. The larger the value of $η$ at a certain point $(xi,Hn )$ , the slower is the decay of the distribution function, thus the higher is the expected gain in accuracy from an increased expansion order.

Since due to Guideline 1 the highest available expansion order is associated with dual boxes in $B˜$ , the rate of decay for a box $B i$ can be based either on the highest expansion order for the box $Bi$ , or on the dual boxes $˜ Bi,j$ overlapping $Bi$ . In the latter case, one obtains for the error estimator $ηi$ for the box $Bi$

[∑ ] B˜i,j log |Ql+1 | ηi ∼ ------N˜-------− log f0,0 ∕log(l + 1) , (6.16) i

where the maximum expansion order in $Bi$ is $l$ , the sum extends over the dual boxes overlapping $Bi$ and $˜Ni$ denotes the number of dual boxes overlapping $Bi$ .

The exact evaluation of the the term $log|Ql| = max θ,φlog|Ql(θ,φ)|$ is a maximization problem over all angles and leads to high overall numerical effort, because such a maximization problem needs to be solved at every grid point $(xi,Hn)$ . A numerically very cheap approximation to the maximum is

l+1 l+1 ∑ ∑ |Ql+1 | ≈ |fl+1,m||Yl,m | ≈ |fl+1,m| , (6.17) m= −l−1 m= −l−1

since the expansion coefficients $fl,m$ are readily available. The approximation $|Yl,m| ≈ 1$ leads to an overestimation of the maximum, but is not a concern because it leads to a constant offset in the error indicator $ηi$ only.

The case $l = 0$ in (6.16) requires additional treatment since $log(l + 1)$ evaluates to zero. Replacing $l$ by $l + 1$ resolves these problems and partly accounts for the overestimation of $|Q | l+1$ by (6.17). For larger values of $l$ , the replacement of $l$ by $l + 1$ does not cause a significant difference anyway. This leads to the final form of the estimator:

[∑ ˜ log[∑Li+1 |fL+1,m(B˜i,j)|] ] ηi ∼ ---Bi,j------m=−-Li−-1---i---------- − log |f0,0| ∕log(Li + 2) , (6.18) N˜i

where $Li$ denotes the maximum even expansion order of the box $Bi$ . After an evaluation of the estimator for all boxes $Bi$ , the expansion order can then be increased for the boxes with largest values of $η i$ . An additional smoothing step then ensures conformity with respect to Guidelines 1 and 2.

The estimator (6.18) is based on the rate of decay of the expansion coefficients only and does not consider the exponential decay of the distribution function at higher energies. Therefore, the estimator treats regions with large values of the distribution function as being equally important as regions with very low values of the distribution function. This is a disadvantage if the target quantity is a macroscopic quantity such as the average carrier velocity, for which only large values of the distribution function lead to significant contributions. In such situations, an additional term $α logf 0,0$ may be added to (6.18) in order to penalize low-probability regions of the distribution function. In the special case $α = 1∕ log(Li + 2)$ , this modification is equivalent to dropping the term $logf0,0$ from (6.18). It should be noted that this modification towards higher weight at regions with high probability is similar in spirit to a target quantity driven control of the expansion order discussed next.

6.2.2 Target Quantity Driven Adaptive Control

Another approach to the adaptive control of the expansion order is based on one or more dedicated target quantities, for which the expansion orders are chosen such that a prescribed accuracy is obtained. For instance, consider as target quantity the average particle energy

∫ ∞ ⟨ε⟩(x) = εf(x,ε)Z (ε)dε . (6.19) 0

The average energy of a discrete solution $f(xi,Hn )$ of the SHE equations is computed using e.g. a simple midpoint quadrature rule as

N∑H ⟨ε⟩(xi) ≈ ε(Hn )f0,0(xi,Hn)Z (Hn)(H+n − H −n ) . (6.20) n=0

Hence, this allows for an extraction of the individual weight $w(xi,Hn )$ of the grid node $(xi,Hn )$ on the target quantity as

ε(Hn )f0,0(xi,Hn )Z(Hn )(H+n − H−n ) w (xi,Hn ) = -------------⟨ε⟩(x-)------------- . (6.21) i

Weights for other target quantities such as the average carrier velocity or the carrier density can be defined in a similar manner. Therefore, $w (xi,Hn )$ can be seen as a generic weight function for an arbitrary target quantity $q(x, H)$ , where typically macroscopic target quantities of the form $q(x)$ are of interest.

An increase of the expansion order is most appropriate in regions where the weights $w$ are high. On the contrary, regions with a weight of several orders of magnitude below the boxes of highest weight can be kept at lowest expansion order, because the contribution to the target quantity is very small anyway.

(a) Plot of the two distribution functions.

(b) Plot of the normalized weight for the computation of average carrier energies from the two distribution functions.

Figure 6.4: Comparison of weights for the computation of average energies for a Maxwellian distribution $f1$ and a shifted, heated Maxwellian $f2$ . The two weights are normalized to a peak value $1$ .

A comparison of the weight functions for the calculation of the average carrier energy of a Maxwellian distribution $f 1$ and a shifted and heated Maxwellian distribution $f 2$ is depicted in Fig. 6.4. For the Maxwellian distribution it is sufficient to consider only energies up to about $0.3$ eV for the computation of the average energy at an accuracy of about one percent, because the weights at higher energies are below $10−3$ already. However, for the shifted and heated Maxwellian distribution, contributions up to an energy of about $0.6$ eV need to be considered for comparable accuracy.

6.2.3 Residual-Based Adaptive Control

An expansion order control based solely on the weight function $w$ does not account appropriately for long-range influences of regions of the distribution function with low probability on regions with higher probability. This is a particular concern in short devices at high electric fields, where strong variations of the distribution function with respect to energy are observed. Similar to the estimation of the rate of decay in Sec. 6.2.1, a residual-based approach provides a better monitor for the need for higher expansion orders than the target quantity driven approach in the previous section.

Let $L L f = (fl,m (xi,Hn))l,m;i,n$ denote the numerical solution of the SHE equations for a possibly variable expansion order $L = L(xi,Hn )$ obtained from the solution of the system

SLf L = bL . (6.22)

Now, let $L`→L ′ L`→L ′ f = (fl,m (xi,Hn ))l,m;i,n$ denote the prolongated solution $L f$ for a possibly variable expansion order $′ ′ L = L (xi,Hn ) > L (xi,Hn )$ obtained by

{ L`→L′ fLl,m(xi,Hn ) , l ≤ L(xi,Hn ) fl,m (xi,Hn) = 0 , otherwise . (6.23)

Denoting the system matrix for the expansion orders $′ L$ with $L′ S$ and the right hand side with $L′ b$ , the residual

eL`→L′ := SL ′fL`→L ′ − bL′ (6.24)

provides an indication on where to increase the expansion order. In order to make the residual invariant with respect to the scale of the values $fL (xi,Hn ) l,m$ , the linear system should be first symmetrized and rescaled as discussed in Sec. 7.2. Note that the matrix compression scheme derived and discussed in Chap. 4 allows for a convenient means to store the larger system matrix $′ S$ and to compute the matrix-vector product $L′ L`→L′ S f$ efficiently.

It is important to keep in mind that the residual is not invariant with respect to transformations of the linear system. In particular, if the $j$ -th equation is multiplied with a certain factor, then also the residual is multiplied with the same factor, even though the solution of the linear system is unchanged. Consequently, it is advisable to first apply normalizations, for instance by multiplication of each equation in the system with the reciprocal of the respective diagonal entry.

Similar to other expansion order control strategies outlined above, an increase of the expansion order can finally be carried out for boxes with highest residuals. An additional expansion order smoothing step after the increase ensures that Guidelines 1 and 2 are obeyed.

6.3 Results

Average carrier velocities along a $100$ nm $n+nn+$ -diode with a bias of $0.7$ Volt and intrinsic region between $x = 20$ nm and $x = 60$ nm are compared for different uniform expansion orders as well as for the three adaptive expansion order schemes after one (maximum SHE order 3), two (maximum SHE order 5) and three (maximum SHE order 7) adaption steps. The total energy range spans $2$ eV using an energy spacing of $12.5$ meV. A first-order expansion is kept directly at the band edge, because it has been observed that it improves numerical stability at high expansion orders. At the contacts, a Maxwell distribution is imposed as a Dirichlet boundary condition, thus the expansion is kept at first order there. The resulting velocity curves are depicted in Fig. 6.5 and show rather small differences between the different expansion orders. Nevertheless, the logarithmic plots of the relative errors in Fig. 6.9 provide full insight.

Figure 6.5: Comparison of carrier velocities for different uniform and adaptive SHE orders. The differences between the adaptive schemes are below the image resolution and analyzed in more detail in Fig. 6.9.

The adaptive strategy based on the decay of expansion coefficients as described in Sec. 6.2.1 is illustrated in Fig. 6.6. In order to emphasize refinement in high-probability regions, the term $0.25× log(f0,0)$ is added to (6.18) as discussed above. The resulting indiciator values are shifted such that the highest value is zero. The threshold for an expansion order increase is set to $− 6.0$ . The indicator is not computed at the contact, but evaluated to $0$ in Fig. 6.6. During the adaption procedure, the adaptive control stays at a first-order expansion in the left $n+$ region, where the distribution function is still close to equilibrium. The expansion order is then increased in the intrinsic region and away from the band edge at higher energies after the intrinsic regions. In summary, the adaptive strategy based on the decay of the expansion coefficients increases the expansion order mostly along the high energy tail and close to the band edge inside the intrinsic region. The increase of the expansion order near the right contact stems from the use of Dirichlet boundary conditions and can be considered to be an artifact, because the distribution function is forced from a heated distribution to a Maxwell distribution at the contact, leading to an unphysical boundary layer as discussed by Schroeder et al. [92].

A different behavior is observed for the target quantity based control shown in Fig. 6.7. The indicator is obtained by taking the logarithm of the contribution of the respective box in the simulation domain and shifting the values such that the highest contribution leads to an indicator value of zero. A threshold value of $− 2.0$ has been used for increasing the expansion order. Since the distribution function changes its shape only mildly at higher energies, the indicator is essentially unchanged during the adaption, leading to higher expansion orders near the band edge only. Note that the high-energy tail of the distribution function right after the intrinsic region is resolved by the increased expansion order. The error plot in Fig. 6.9(b) shows that virtually the same accuracy as for uniform expansions is obtained. However, the expansion order at later adaption steps abruptly changes from first-order to highest order as quickly as possible without violating Guidelines 1 and 2, thus it is reasonable to expect that a less abrupt change of the expansion order can preserve the accuracy using a lower number of unknowns.

Fig. 6.8 depicts the effect of the residual based expansion order control. The residual has been normalized first by scaling the uniform system such that a unit diagonal is obtained. The indicator is then obtained by taking the logarithm of the residual values and applying a suitable shift such that the highest residual value leads to an indicator value of zero. The expansion order is increased if the indicator is larger than $− 2.0$ . In principle, one may assign less weight to higher energies such as for the adaptive strategy based on the rate of decay of expansion coefficients. For illustration purposes, the unmodified indicator is shown. It should be noted that the strategy based on the decay leads to qualitatively similar indicator values if no modification at higher energies is applied. The pure residual based strategy increases the expansion order to third-order up to high energies in the right half of the device. However, fifth-order expansions are used in a very small region in the device only. A seventh-order expansion is only assigned right above the band edge inside the intrinsic region. The error plot in Fig. 6.9(c) is very similar to the error plot for the strategy based on the decay of expansion coefficients, hence the same conclusions can be drawn.

Fig. 6.9 further shows that the SHE method converges to a solution as the expansion order increases. As reference, a uniform seventh-order expansion is used. A first-order expansion shows an average error of about $10− 2$ over the device with respect to the seventh-order expansion taken as reference. The average error of a third-order expansion is around $− 3.5 10$ , and approximately $− 5 10$ for a fifth-order expansion. A similar exponential decay of the error with the SHE order has also been observed by Jungemann et al. [53] for the collector current of a bipolar junction transistor, even though at a lower rate.

The number of unknowns obtained with uniform and three adaptive expansions are depicted in Fig. 6.10. Savings of a factor of almost three are obtained with the target quantity based scheme, which yields virtually the same accuracy as uniform expansions. Memory requirements and execution times are consequently reduced by the same factor. Typically, higher savings for the execution times are obtained in practice, because linear solvers with optimal complexity are not available in general. A combination of the three different schemes proposed may give slightly higher savings, which become significant at very high expansion orders only.

(a) Indicator for the first adaptation step (left). Expansion order after the first adaption step (right).

(b) Indicator for the second adaptation step (left). Expansion order after the second adaption step (right).

Figure 6.6: Expansion indicator and expansion orders in the $n+nn+$ diode for three adaption steps using an adaptive scheme based on the decay of expansion coefficients as outlined in Sec. 6.2.1. The scheme starts with a uniform first-order expansion. The error indicator is not computed on the left and right contact and thus plotted with an arbitrary value of $0$ .

(a) Indicator for the first adaptation step (left). Expansion order after the first adaption step (right).

(b) Indicator for the second adaptation step (left). Expansion order after the second adaption step (right).

Figure 6.7: Expansion indicator and expansion orders in the $n+nn+$ diode for three adaption steps using an adaptive scheme based on a target quantity as discussed in Sec. 6.2.2. The scheme starts with a uniform first-order expansion.

(a) Indicator for the first adaptation step (left). Expansion order after the first adaption step (right).

(b) Indicator for the second adaptation step (left). Expansion order after the second adaption step (right).

Figure 6.8: Expansion indicator and expansion orders in the $n+nn+$ diode for three adaption steps using an adaptive scheme based on residuals as described in Sec. 6.2.3. The scheme starts with a uniform first-order expansion.

(a) Rate of decay based scheme.

(b) Target quantity based scheme.

Figure 6.9: Comparison of the relative errors in carrier velocities for different uniform and adaptive SHE orders using different adaptive schemes. A uniform seventh-order expansion is used as reference.

Figure 6.10: Comparison of the number of unknowns for the three different adaptive expansion order control schemes.

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

Chapter 6Adaptive Variable-Order SHE

6.1 Variable-Order SHE

6.2 Adaptive Control of the SHE Order

6.2.1 Rate of Decay of Expansion Coefficients

6.2.2 Target Quantity Driven Adaptive Control

6.2.3 Residual-Based Adaptive Control

6.3 Results

Chapter 6
Adaptive Variable-Order SHE