λ2 = lε2sqCVTm, ǫ2 = 6m∗2
nqVTl2, δ2 = CnI
m, δ2n = δ2Dβ4(φ)
n , δ2p = −δ2βD5(φ)
p , βn = Dln2τ, βp = Dl2pτ, ρ2n= κnqDTmD6(ϕn)
nVTCm , ρ2p = κpqDTmD7(ϕp)
pVTCm .
(3.29)
The boundary conditions (2.21)-(2.30) can be similarly scaled.
4 Energy Relaxation with Non-Parabolic Band Approximation
We have to impose some physical assumptions on our model. The as-sumptions (H1)-(H2) below are imposed in order to get simpler expressions for the variables. The non-degeneracy assumption (H3) is valid for semicon-ductor devices with a doping concentration which is below 1019cm−3. Almost all devices in practical applications satisfy this condition. The assumptions (H4) is the non-parabolic band structure in the sense of Kane [14] which is to abide by the law of conservation of momentum. All these assumptions follow from that of [14].
(H1) The energy-band diagram ε of the semiconductor crystal is spher-ically symmetric and a strictly monotone function of the modulus k = −→k
of the wave vector k. Therefore, the Brillouin zone equals ?3 and ε :? −→?, k −→ ε(k).
(H2) A momentum relaxation time can be defined by
τ (ε) = (φ0(2N0+ 1)εβN(ε))−1, β > −2, φ0 > 0, (4.1) where N(ε) = 4πk2/ |ε′(k)| is the density of states of energy ε = ε(k) [2] and N0 is the phonon occupation number [3].
(H3) The electron density n and the internal energy E are given by non-degenerate Boltzmann statistics.
(H4) Let the energy ε(k) satisfy
ε(1 + αε) = k2
2m∗ (4.2)
The constant m∗is the effective electron mass given by m∗ = m0kBT0/2k02, where m0is the unscaled effective mass, k0is a typical wave vector, and α > 0 is the non-parabolic parameter. Notice that we get a parabolic band diagram if α = 0.
The energy relaxation terms (2.6) and (2.7) can be written as the following general form [14]:
W = T2β−1 µ0τnw
rβ(αT )( g1
pβ(αT, 2) − g2
pβ(αT, 3)). (4.3) Notice that the function rβ is in fact a polynomial:
rβ(αT ) = Γ(β + 2) + 5Γ(β + 3)αT + 8Γ(β + 4)(αT )2+ 4Γ(β + 5)(αT )3. (4.4) The symbol Γ denotes the Gamma function defined by
Γ(s) =
∞
0
us−1e−udu, s > 0. (4.5)
The assumption (H4) implies γ(T u) = 2m∗T u(1 + αT u) [14]. Define the functions
pβ(αT, l) =
∞
0
1 + αT u
(1 + 2αT u)2ul−β−1e−udu, (4.6) q(αT, l) =
∞
0
(1 + αT u)12(1 + 2αT u)u12+le−udu. (4.7) By letting αT u = x, we have
1 + αT u
(1 + 2αT u)2 = 1 + x
(1 + 2x)2 = (1 + x) 1
(1 + 2x)2 (4.8)
= (1 + x)(1 − 4x + 12x2− 32x3+ · · · )
= 1 − 3x + 8x2− 20x3+ · · · and hence
pβ(αT, l) =
∞
0 [1 − 3αT u + 8(αT u)2− 20(αT u)3+ · · · ]ul−β−1e−udu (4.9)
Similarly, we have
(1 + αT u)12(1 + 2αT u) = (1 + x)12(1 + 2x) (4.10)
= [1 + 1/2
1 x +1/2 · (−1/2)
1 · 2 x2+ 1/2 · (−1/2) · (−3/2)
1 · 2 · 3 x3+ · · · ] (1 + 2x)
= (1 + 1 2x − 1
8x2+ 1
16x3+ · · · )(1 + 2x)
= 1 + 1 2x − 1
8x2 + 1
16x3+ 2x + x2 −1
4x3+1
8x4+ · · ·
= 1 + 5 2x + 7
8x2− 3
16x3+ · · · and
q(αT, l) =
∞
0
(1 + 5
2αT u +7
8(αT u)2− 3
16(αT u)3 + · · · )u12+le−udu (4.11) The functions g1 and g2 in (4.3) can be computed in terms of n and T . Under the assumptions (H1)-(H3), we can write
g1(n, T ) = µ(1)β (T )T n, g2(n, T ) = µ(2)β (T )T2n (4.12) with the temperature-dependent mobilities
µ(i)β (T ) = µ0pβ(αT, i + 1)
q(αT, 0) T−1/2−β, i = 1, 2. (4.13) Here, the mobility constant µ0 is given by
µ0 = (3πφ0(2N0+ 1)m∗(2m∗)3/2)−1. (4.14)
Finally, we consider two cases of the parameters α and β, namely, (α = 0, β = 1/2) and (α = 1/2, β = 1/2). Note the momentum relaxation time parameter β = 1/2 corresponding to the so-called Chen model which is shown in [14] to be more effective for handling the spurious velocity overshot peak than the so-called Lyumkis model with β = 0. Note also that, by (4.2), the case of α = 0 corresponds to the parabolic band structure whereas α = 1/2 corresponds to non-parabolic band structure.
Case 1. (α = 0, β = 1/2) Parabolic band.
From (4.3), we have W1s= 1
µ0τnw
r1/2(0)( g1
p1/2(0, 2) − g2
p1/2(0, 3)) (4.15)
s denotes after scaling, namely, the term in (3.27), and
g1 = µ(1)1/2(Ts)Tsns, (4.16(a)) g2 = µ(2)1/2(Ts)Ts2ns,
µ(1)1/2(Ts) = µ0p1/2(0, 2)
q(0, 0) Ts−1, µ(2)1/2(Ts) = µ0p1/2(0, 3)
q(0, 0) Ts−1, (4.17)
p1/2(0, 2) = q(0, 0) =
∞
0
u1/2e−udu =
√π
2 , (4.18(a)) p1/2(0, 3) = r1/2(0) =
∞
0
u3/2e−udu = 3√ π
4 . (4.18(b)) Hence, we have
g1 = µ(1)1/2(Ts)Tsns = µ0Ts−1Tsns = µ0ns, (4.16(b)) where W1 is unscaled.
The energy-transport equation expressed in [14] is given by
−div J2s = −J1s· ▽sVs+ Wns. (4.21) Now, we unscale it to obtain
−ldivJ2 l
where VT = kBTL/q, Dn= unVT and J1is the electron current density defined as that of (2.21) in our model.
From (2.6) where the electron average energy is taken as ωn= 3/2kBTn, our energy relaxation term can be written as
nWn = −nωn− ω0
τnw = −n3/2kBTn− 3/2kBTL
τnw
= −3 2kB
n(Tn− TL) τnw
(4.23) Which is exactly the same as (4.22) used in [14]. This verifies the two different scaling formulations between ours and that of [14] are indeed having the same relaxation term.
Case 2. (α = 1/2, β = 1/2) Non-parabolic band.
Again from (4.3), we have
W2s= 1 µ0τnw
r1/2(Ts
2)( g1
p1/2(T2s, 2) − g2
p1/2(T2s, 3)) (4.24) where
g1 = µ(1)1/2(Ts)Tsns, g2 = µ(2)1/2(Ts)Ts2ns, µ(1)1/2(Ts) = µ0p1/2(T2s, 2)
q(T2s, 0) Ts−1, µ(2)1/2(Ts) = µ0p1/2(T2s, 3)
q(T2s, 0) Ts−1, (4.25)
p1/2(Ts From (4.9) and (4.11), we have
p1/2(Ts
q(Ts
Γ(3 2) =
√π 2 , Γ(5
2) = 3 2
√π 2 , Γ(7
2) = 5 2
3 2
√π 2 , Γ(9
2) = 7 2
5 2 3 2
√π 2 , we have the formula
Γ(k + 1
2) = (2k)!√ π
4kk! . (4.34)
The energy-transport equation expressed in [14] is given by
−div J2s = −J1s· ▽sVs+ Wns. Now, we unscale it to obtain
−ldivJ2 l of (2.21) in our model.
5 Algorithm
Our new model for both DG and ET equations with the seven state variables φ, u, v, ζn, ζp, gn, and gp and their associated boundary conditions (BCs) is
re-organized as follows:
∆φ = F (φ, u, v, ζn, ζp), (5.1) 1
q∇ · Jn = R(φ, u, v, ζn, ζp), (5.2) 1
q∇ · Jp = −R(φ, u, v, ζn, ζp), (5.3)
∆ζn = Zn(φ, u, v, ζn, ζp), (5.4)
∆ζp = Zp(φ, u, v, ζn, ζp), (5.5)
∇ · Gn = Rn(gn), (5.6)
∇ · Gp = Rp(gp), (5.7)
The main ingredients of the algorithm solving the DGET model are adap-tive finite element approximation of the model, node-by-node and monotone iterative solution of the resulting nonlinear algebraic systems, and Gummel’s iteration consecutively on the PDEs as described in [5] for the ET model.
For the sake of clearness, we briefly illustrate the algorithm and refer to [5, 6] for more details on the adaptive finite element formulation, monotone convergence analysis, and practical implementation issues.
Here we use the notation l as Gummel’s (outer) iteration index and m as the monotone (inner) iteration index.
Step 1. Initial Mesh: Create a coarse and structured mesh for which the number of nodes can be chosen as small as possible.
Step 2. Preprocessing: See [5].
Step 3. Gummel and Monotone Iterations on (5.1)-(5.7).
Step 3.0. Set l := 0
Step 3.1. Solve the potential equation in (5.1).
Step 3.1.1. Set m := 0 and set the initial guess φ(m)j =
φj or φj if l = 0, φ(l)j otherwise,
for all (xj, yj) ∈ Ωh,
where φj and φj are constant values that can be easily verified to be an upper and lower solution of φ, respectively, and Ωh denotes the set of mesh points on the closure of the domain.
Step 3.1.2 If l = 0, set u(l) and v(l) by the charge neutrality condition.
Step 3.1.3. Compute φ(m+1)j by solving the discrete potential system of (5.1)
ξjφ(m+1)j + γj(φ) φ(m+1)j =
k∈V (j)ξkφ(m)k
−F (φ(m)j , u(l)j , vj(l), ζ(l)n , ζ(l)p ) + γj(φ) φ(m)j , ∀(xj, yj) ∈ Ωh, φ(m+1)j = VO+ Vb, ∀(xj, yj) ∈ ∂ΩhD,
∂φ(m+1)j
∂n = 0, ∀(xj, yj) ∈ ∂ΩhN ,
(5.8) where
γj(φ) = max
∂F (φj)
∂φ ; ˆφj ≤ φj ≤ ˜φj
, (5.9)
ξk are the matrix elements of the discretization, and Ωh, ∂ΩhD, and
∂ΩhN represent the sets of mesh points in the interior, Dirichlet part, and Neumann part of the domain, respectively.
Step 3.1.4. Set φ(m)j := φ(m+1)j ∀ j and m := m + 1. Go to Step 3.1.3 until the stopping criteria of the inner iteration are satisfied.
Step 3.1.5. Set φ(l+1)j := φ(m+1)j ∀ j.
Step 3.2. Solve the electron continuity equation (5.2).
Step 3.2.1. Set m := 0 and set the initial guess u(m)j =
uj or uj if l = 0, u(l)j otherwise,
for all (xj, yj) ∈ Ωh, where uj andu j are constant values for all (xj, yj) ∈ Ωh that can be easily verified to be an upper and lower solution of u, respectively.
Step 3.2.2. Compute u(m+1)j by solving the discrete electron system of (5.2).
Step 3.2.3. Set u(m)j := u(m+1)j ∀ j and m := m + 1. Go to Step 3.2.2 until the stopping criteria of the inner iteration are satisfied.
Step 3.2.4. Set u(l+1)j := u(m+1)j ∀ j.
Step 3.3. Solve the hole continuity equation (5.3) similarly to that in Step 3.2.
Step 3.4. Solve the DG equation (5.4).
Step 3.4.1. Set m := 0 and set the initial guess [ζn](m)j =
[ζn]j or ![ζn]j if l = 0, [ζn](l)j otherwise,
for all (xj, yj) ∈ Ωh, where [ζn]j ≈ ζn(xj, yj) and [ζn]j and ![ζn]j are constant values for all (xj, yj) ∈
Ωh that can be easily verified to be an upper and lower solution of ζn, respectively.
Step 3.4.2. Compute [ζn](m+1)j by solving the discrete system of (5.4).
Step 3.4.3. Set [ζn](m)j := [ζn](m+1)j ∀ j and m := m + 1. Go to Step 3.4.2 until the stopping criteria of the inner iteration are satisfied.
Step 3.4.4. Set [ζn](l+1)j := [ζn](m+1)j ∀ j.
Step 3.5. Solve the DG (5.5) similarly to that in Step 3.4.
Step 3.6. Update [φqn](l+1)j and [φqp](l+1)j by the equations (2.16)-(2.17).
Step 3.7. Set l := l + 1 and go to Step 3.1 until the stopping criteria of the outer iteration are satisfied.
Step 4. Monotone Iteration on (5.6) and (5.7).
Step 4.1. Solve the energy equation (5.6) for gn similarly to that in Step 3.2.
Step 4.2. Solve the energy equation (5.7) for gp similarly to that in Step 3.2.
Step 5. Error Estimation: See [5].
Step 6. Refinement: See [5]
Step 7. Postprocessing: All computed solutions are then postprocessed for further analysis of physical phenomena.
Note that, in each one of Steps 3.1-3.5, 4.1, and 4.2, a Jacobi (node-by-node) type of solution is performed for the corresponding discrete system (5.8), for example, in which the monotone parameters (5.9) can be easily obtained by means of lower and upper solutions. Two important factors that guarantee a global convergence with this kind of simple solutions as initial guesses are the diagonally dominant property of the matrices due to the self-adjoint formulation and the monotonicity of the parameters by the special nonlinearity of the formulation. The diagonally dominant property for (5.1)-(5.7) can proved in exactly the same manner as that given in [5, 6].
It can also be easily shown that each one of the nonlinear functionals in (5.1)-(5.7) is monotone in its respective state variable. It is thus a straightforward generalization from our previous theoretical analysis that all the nonlinear algebraic systems generated by this algorithm preserve these two factors.