National University of Kaohsiung Repository System:Item 310360000Q/10850

全文

(1)國立高雄大學應用數學研究所碩士論文. 非拋物型能帶對量子校正的能量運輸模型中之能量鬆弛之影響第二部份:模擬 Non-Parabolic Band Effects on Energy Relaxation for the Quantum-Corrected Energy Transport Model Part II : Simulation. 研究生：任家進撰指導教授：劉晉良. 中華民國 99 年 7 月.

(2)

(3) 致謝這篇論文的完成首先需要感謝指導教授劉晉良老師的教導，這兩年的時間承蒙老師不厭其煩的細心教導使本篇論文更加的完善也讓我在這段時間學習到更多知識，並且從老師對追求知識的態度及待人處事上更是讓我獲益良多，僅此致上我最誠摯的敬意及感謝。本篇論文的完成也要感謝高雄師範大學教授陳仁純老師，感謝陳仁純老師這段時間特意撥空指導程式方面的編寫，也由於陳仁純老師這段時間花費大量時間幫忙修正我在程式編寫上的錯誤以及提供寶貴的意見，才使得本篇論文能順利完成。在此學生致上誠摯的謝意。在就學的兩年時間也必須感謝高雄大學應用數學系的每位老師和系辦小姐的教導與協助，感謝每位老師這段時間上在課業的教導，讓我在課業上有更多的認識。另外也要感謝各位同學的陪伴，這兩年的相處點滴都將會是我最珍貴的回憶之一。最後，要特別感謝我的家人，由於你們鼎力的支持、加油、包容，我才能無後顧之憂的完成這兩年的學業。再次感謝所有幫助我關心我的人，謝謝你們。.

(4) Non-Parabolic Band Effects on Energy Relaxation for the Quantum-Corrected Energy Transport Model Part II : Simulation. by Jia-Chin Jen Advisor Jinn-Liang Liu. Department of Applied Mathematics, National University of Kaohsiung Kaohsiung, Taiwan 811, R.O.C. July 2010.

(5) Contents 1. Introduction. ……………………………………………………….1. 2 The QCET Model. 3. ………………………………………………....3. Singularly Perturbed Formulation of QCET Model. ………………10. 4 Energy Relaxation with Non-Parabolic Band Approximation. 5. Algorithm. 6. Numerical Results. 7. Conclusions. ……16. ………………………………………………………..26. ………………………………………………...31. …………………………………………………….…37.

(6) 非拋物型能帶對量子校正的能量運輸模型中之能量鬆弛之影響第二部份:模擬之能量鬆弛之影響第二部份模擬指導教授：劉晉良教授國立新竹教育大學應用數學系. 學生：任家進國立高雄大學應用數學所. 摘要. 量子校正的能量運輸模型是由描述穩定狀態中電子和電洞流的七條非線性對稱的偏微分方程式所組成，它們的能量運輸，古典和量子位能是在一個等比縮放的奈米半導體元件中產生。我們延伸量子校正的能量運輸模型中包含非拋物線型能帶結構之電子與電洞的能量鬆弛項。我們參考 Degond 等人的論文得到我們詳細的量子校正的能量傳輸模型的能量鬆弛項公式。在數值結果說明了在拋物型能帶中電子與電洞的最大溫度大約 3500K，非拋物型能帶大約 1300K。關鍵字：量子校正、半導體、能量鬆弛、非拋物型能帶結構關鍵字.

(7) Non-Parabolic Band Effects on Energy Relaxation for the Quantum-Corrected Energy Transport Model Part II : Simulation Advisor: Professor Jinn-Liang Liu Department of Applied Mathematics National Hsinchu University of Education. Student: Jia-Chin Jen Institute of Department of Applied Mathematics National University of Kaohsiung. ABSTRACT The quantum-corrected energy transport (QCET) model consisting of seven self-adjoint nonlinear PDEs describes the steady state of electron and hole flows, their energy transport, and classical and quantum potentials within a nano-scale semiconductor device. We extend the energy relaxation terms of the QCET model to include the non-parabolic band structure for both electron and hole. We give explicit formulas of the energy relaxation terms proposed by Degond et al. in [14] for our QCET model. Numerical results have shown that the electron and hole temperatures reduce from about 3500 K for the parabolic case to about 1300 K for the non-parabolic case. Keywords: quantum, semiconductor, energy relaxation, non-parabolic band structure.

(8) 1. Introduction. This thesis is a continuation of the theoretical work presented in [10]. Semiconductor devices can be simulated by means of the semiconductor Boltzmann equation. However, this method is too costly and time consuming to model real problems in semiconductor applications. Acceptable accuracy can be reached by solving macroscopic semiconductor models which derived from the Boltzmann equation. The simplest models are drift-diffusion models which consist of the mass continuity equation for the charge carriers and a definition for the particle current density [15]. These models, however, are not accurate enough for submicron device modeling, owing to the rapidly changing fields and temperature effects [14]. The quantum-corrected energy transport (QCET) model consisting of seven self-adjoint nonlinear PDEs describes the steady state of electron and hole flows, their energy transport, and classical and quantum potentials within a nano-scale semiconductor device [7]. We extend the energy relaxation term of QCET model to include the non-parabolic band structure for electron. Since the QCET equations are of parabolic type, the numerical solution needs less effort than the hydrodynamic models. Moreover, the QCET equations can be written in a drift-diffusion formulation, therefore the numerical effort is comparable to the drift-diffusion models [14]. In this thesis, we compute explicitly about the energy relaxation term in terms of the temperature. For non-parabolic bands in the sense of Kane [14],. 1.

(9) the coefficients can be computed analytically. We use the Gamma function to compute the energy relaxation term. We assume that the energy-band diagram of the semiconductor crystal is → − spherically symmetric and monotone in the modulus of the wave vector k , that non-degenerate Boltzmann statistics can be used and that a momentum relaxation time τ can be defined by τ (ε) ∼ ε−β N (ε)−1 , where ε is the energy, N(ε) denotes the density of states, and β > −2 is a parameter [14]. Then, using the general formulas for the coefficients and densities from [2], we get more explicit expressions than those of [2], involving the energy-band function → − ε( k ) and depending on the temperature. Furthermore, we get analytical expressions under the additional assumption of non-parabolic bands in the sense of Kane [14]:. ε(1 + αε) =. →2 2 − k , 2m0. where is the reduced Planck constant, m0 the effective electron mass, and α > 0 the non-parabolicity parameter. The numerical experiments are performed by employing the QCET model which is defined by different non-parabolic parameter α. The numerical results in [14] show that the spurious velocity overshoot spike at the anode junction becomes smaller in the non-parabolic band case (α = 1/2), compared to the parabolic case (α = 0). We give explicit formulas of the energy relaxation terms proposed by Degond et al. in [13] for our QCET model. Numerical results have shown. 2.

(10) that the electron and hole temperatures reduce from about 3500 K for the parabolic case to about 1300 K for the non-parabolic case.. 2. The QCET Model. Based on the density gradient (DG) theory [1], the QCET model proposed in [7] is −∆φ = f1 (φ) :=. . q  εs. φ+φ −nI exp( VTqn )u −φ−φ +nI exp( VT qp )v. +C.  . q (−n + p + C) , (2.1) εs −ζ n −∆ζ n = f2 (ζ n ) := VT ln(ζ 2n ) − VT ln(nI u) − φ , (2.2) 2bn ζp −∆ζ p = f3 (ζ p ) := −VT ln(ζ 2p ) + VT ln(nI v) − φ ,(2.3) 2bp −∇ · D4 (φ)∇u = f4 (u) . φ −φ q neq peq − n2I exp( qnVT qp )uv : =   −φ−φ φ+φ nI exp( VT qp )v + nI exp( VTqn )u  τ √ εt −εi +2 neq peq exp( kB TL ) =. =. q (neq peq − np). , √ i τ n + p + 2 neq peq exp( kεtB−ε ) TL. −∇ · D5 (φ)∇v = f5 (v) := −f4 (u) ,. −∇ · D6 (ϕn )∇gn = f6 (gn ) := Jn ·E + nWn , −∇ · D7 (ϕp )∇gp = f7 (gp ) := Jp ·E + pWp ,. 3. (2.4) (2.5). ωn − ω0 , (2.6) τ nω ωp − ω0 Wp := − , (2.7) τ pω Wn := −.

(11) with the seven unknown functions φ, ζn =. √ n,. (2.8). ζp =. √ p,. (2.9). φ + φqn −ϕn ), n = nI exp( )u, VT VT −φ − φqp ϕp )v, v = exp( ), p = nI exp( VT VT 5ϕ gn = Tn exp(− n ), 4VT 5ϕp ), gp = Tp exp( 4VT u = exp(. (2.10) (2.11) (2.12) (2.13). and the auxiliary relations E = −∇φ. (2.14). φqn = VT ln(ζ 2n ) − VT ln(unI ) − φ,. (2.15). φqp = −VT ln(ζ 2p ) + VT ln(vnI ) − φ, φ + φqn ), D4 (φ) = qDn nI exp( VT −φ − φqp ), D5 (φ) = −qDp nI exp( VT 5ϕ D6 (ϕn ) = κn exp( n ), 4VT 5ϕp D7 (ϕp ) = κp exp(− ), 4VT Jn = −qµn n∇(φ + φqn ) + qDn ∇n = D4 (φ)∇u,. (2.16). Jp = −qµp p∇(φ + φpn ) − qDn ∇p = D5 (φ)∇v,. (2.22). (2.17) (2.18) (2.19) (2.20) (2.21). Gn = D6 (ϕn )∇gn ,. (2.23). Gp = D7 (ϕp )∇gp ,. (2.24) 4.

(12) where φ is the electrostatic potential, n and p the electron and hole densities, ϕn and ϕp the generalized quasi-Fermi potentials, Tn and Tp the electron and hole temperatures, φqn and φqp the quantum potentials, E the electric field, Jn and Jp the electron and hole current densities, κn and κp the heat conductivities, εs the permittivity constant of the semiconductor, C the doping profile (impurity concentration), bn =. h2 12m∗n q. and bp =. h2 12m∗p q. the material. parameters measuring the strength of the gradient effects in the electron gas [1], τ nω and τ pω the carrier energy relaxation times, ω n and ω p the carrier average energies, εt trap energy, εi intrinsic energy, and other symbols with their values given in Table 1. The system (2.1)-(2.7) models the steady state of electron and hole flows through the device by augmenting the macroscopic energy transport model (2.1), (2.4)-(2.6) [5] with the DG equations (2.2) and (2.3) [1]. The square roots of carrier densities in (2.8) and (2.9) were introduced in [19] as extra unknown functions to define the quantum (Bohm) potentials (2.15) and (2.16) by means of the generalized quasi-Fermi potentials ϕn and ϕp in (2.10) and (2.11) [7, 19]. These quantum potentials represent the first order quantum corrections of the drift-diffusion fluxes as shown in (2.17) and (2.18). We observe from (2.2), (2.3), (2.15), and (2.16) that the QCET reduces to the classical ET model of [5] in the semiclassical limit h → 0. This model is also more general than the quantum drift diffusion models considered, for example, in [4, 11, 12, 13, 16, 17, 19] by including the energy transport equations (2.6) and (2.7) to deal with hotspot problems in nano-scale device design [18]. 5.

(13) Note particularly that the right-hand side nonlinear functionals fi , i = 1, · · · 7, in (2.1)-(2.7) are all expressed in terms of their respective unknown variables φ, u, v, ζ n , ζ p , gn , and gp to illustrate that each functional can be straightforwardly differentiated with respect to its variable and that each PDE is semilinear. All functionals are nonlinear due to the Slotboom-type transformations (2.10)-(2.13). Furthermore, these transformations also result in that all divergence operators in the left-hand sides of the system (2.1)-(2.7) are self-adjoint. This self-adjoint and semilinear formulation provides many appealing approximation features such as single finite element subspace for all seven variables, global and optimal convergence, and fast iterative solution [5, 8] with suitable conditioning scalings of the discrete system [6]. The Slotboom-type variables gn and gp in (2.12) and (2.13) were first introduced in [5]. The relation between n and u in (2.10) implicitly defines the Slotboom variable u that generalizes the classical Slotboom variable to include the quantum potential. A similar generalization was also introduced in [17]. Let Ω ⊂ ℜ2 denote the bounded domain of the silicon in Fig. 1. The boundary ∂Ω = ∂ΩO ∪∂ΩI ∪∂ΩN is piecewise smooth with ∂ΩO = BC∪DE∪ AF denoting the Ohmic contacts, ∂ΩI = CD the silicon/oxide interface, and ∂ΩN = AB ∪ EF the Neumann boundary parts. The boundary conditions for the unknown functions are. 6.

(14) 7.

(15) φ = VO + Vb on ∂ΩO ,

(16) . 1 2 2 2 ζn = C + C + 4nI on ∂ΩO , 2 ζ p = nI /ζ n on ∂ΩO , ζ p = ζ n = 0 on ∂ΩI , −VO ) on ∂ΩO , u = exp( VT VO v = exp( ) on ∂ΩO , VT 300 on ∂ΩO , gn = O exp( 5V ) 4VT 300 on ∂ΩO , gp = O exp(− 5V ) 4VT ∂ζ p ∂φ ∂u ∂v ∂ζ ∂gn ∂gp = = = n = = = = 0 on ∂ΩN , ∂n ∂n ∂n ∂n ∂n ∂n ∂n ∂u ∂v ∂gn ∂gp = = = = 0 on ∂ΩI , ∂n ∂n ∂n ∂n. (2.25) (2.26) (2.27) (2.28) (2.29) (2.30) (2.31) (2.32) (2.33) (2.34). where VO is the applied voltage, Vb is the built-in potential, and n is an outward normal unit vector to ∂ΩN . The condition (2.28) on the interface is so chosen that we do not consider tunneling effects [4] [7] across the interface.. 8.

(17) Table 1. Physical constants Symbol Physical meaning. Value. Unit. kB. Boltzmann constant. 1.38 × 10−23. J/K. q. elementary charge. 1.602 × 10−19. C. m0. electron rest mass. 9.11 × 10−31. Kg. m∗n. electron effective mass. 0.98 × m0. Kg. m∗p. hole effective mass. 0.16 × m0. Kg. ε0. permittivity of vacuum. 8.85 × 10−14. F/cm. εs. silicon dielectric constant. 11.7 × ε0. F/cm. . reduced Planck constant. 1.054 × 10−34. J-s. VT. thermal voltage. 0.0259. V. nI. intrinsic carrier concentration 1.5 × 1010. cm−3. µn. electron mobility. 1350. cm2 /V-s. µp. hole mobility. 480. cm2 /V-s. Dn. electron diffusion coefficient. 34.659. cm2 /s. Dp. hole diffusion coefficient. 12.432. cm2 /s. TL. lattice temperature. 300. K. τ. electron (hole) lifetime. 10−8. s [2]. ω0. thermal energy. 6.21 · 10−21. J. 9.

(18) 3. Singularly Perturbed Formulation of QCET Model. Let l be the diameter (characteristic length) of the device. The standard singular perturbation parameter is the scaled Debye length λ defined as λ2 =. εs VT l2 qCm. (3.1). where. Cm = max |C (x, y)| .. (3.2). (x,y)∈Ω. The magnitude of the doping remains roughly the same since the invention of transistor in 1947 [20] whereas the gate length of CMOS technology has been shrunk from 10000 nm (lC = 2.5 · 10−5 m) in 1971 to 45 nm (lQ = 2.5 · 45 · 10−9 m) in 2008. This means that the scaled Debye length which is a crucial parameter in singular perturbation analysis has been enlarged from that of the classical device (with Cm = 1019 cm−3 = 1025 m−3 ). λ2C =. εs VT 11.7 · 8.85 · 10−12 · 0.0259 = = 0.268 · 10−9 2 −5 −19 25 lC2 qCm (2.5 · 10 ) · 1.602 · 10 · 10. (3.3). to that of the modern quantum device (with Cm = 1025 m−3 ). λ2Q. εs VT l2 = 2 = C2 ·λ2C = lQ qCm lQ. . 2.5 · 10−5 2.5 · 45 · 10−9 10. 2. ·0.268·10−8 = 0.132·10−3 . (3.4).

(19) Apparently choosing the scaled Debye length as a singular perturbation parameter is no longer valid for nano-devices. However, singular effects are still present for such devices as shown by the numerical results below. The question is thus what kind of parameter is more feasible for nano-devices. Consider the following scaled intrinsic carrier concentration as a parameter. δ2 =. nI . Cm. (3.5). For both classical and quantum devices, we have. δ 2Q = δ 2C =. 1.5 · 1010 · 106 = 0.15 · 10−8 1025. (3.6). which is obviously more suitable for being a singular perturbation parameter. From the quantum potential equations (2.2) and (2.3), we have another parameter, namely, the scaled Planck constant ǫ defined by. 2. ǫ =.   . 2bn VT l 2. =. 2bp VT l 2. =. 2 6m∗n qVT l2 2 6m∗p qVT l2. for electron,. (3.7). for hole.. Since the effective masses of the electron and hole are of the same order, the scaled Planck constant is denoted by this single symbol. It’s values for the classical and quantum devices are. 1.0542 · 10−68 6 · 0.98 · 9.11 · 10−31 · 1.602 · 10−19 · 0.0259 · 2.52 · 10−10 = 0.8 · 10−9 (3.8) 2 l2 2.5 · 10−5 = C2 · ǫ2C = · ǫ2C = 0.395 · 10−4 (3.9) lQ 2.5 · 45 · 10−9. ǫ2C =. ǫ2Q. 11.

(20) The scaled Planck constant is moderately small indicating that for the 45nm device the quantum effect cannot be ignored any more. We introduce the following scalings where the new dimensionless quantities are marked with the subscript s x y xs = , ys = , ∇s = l l   φ =    s                      . φqns =.    . ϕms. ns =. n , Cm. nIs =. nI Cm. us =. φ VT. , ϕns =. φqn , VT. ϕn , VT. ∂ ∂ , ∂xs ∂ys. . (3.10). ϕp , VT. ϕps =. φqp , VT. Es = −∇s φs , = ϕVm , ϕm = max(x, y)∈Ω |ϕn | , ϕp , T. ps =. φqps =. . p , Cm. ζ ns =. = δ 2 , Cs =. C , Cm. ζ √n , Cm. ζ ps =. ζ √p , Cm. neq , Cm. peqs =. neqs =. (3.11). peq , Cm. u exp(−ϕms ). = exp (−ϕns + ϕms ) , v vs = exp(ϕ = exp ϕ − ϕ , ps ms ms ) n = Cm δ 2 exp φs + φqns − ϕms us , p = Cm δ 2 exp −φs − φqns + ϕms vs.                φqns = ln(ζ 2ns ) − ln(δ 2 us ) + ϕms − φs ,       φqps = − ln(ζ 2ps ) + ln(δ 2 vs ) + ϕms − φs ,   4 (φs ) = exp φs + φqns − ϕms ,  D       D 5 (φs ) = − exp −φs − φqps + ϕms ,  4 (φs )∇s us ,  Jns = qDlJnnnI = D       Jps = lJp = D 5 (φs )∇s vs , qDp nI 12. (3.12). (3.13).

(21)                          . Tns =. Tn , Tm. Tps =. Tp , Tm. TL , Tm. Tm = max(x, y)∈Ω {Tn , Tp } , gns = T expgn−5ϕm = Tns exp −5 (ϕns − ϕms ) , 4 m ( 4 ) gp gps = T exp 5ϕm = Tps exp 54 ϕps − ϕms , m ( 4 ) 6 (ϕns ) = exp 5 (ϕns − ϕms ) , D 4 7 (ϕps ) = exp −5 ϕps − ϕms , D  4     lG 6 (ϕns )∇s gns ,  Gns = κn Tnm = D      p  7 (ϕps )∇s gps , Gps = κlG =D   p Tm    l2 Wp  l2 W   Wns = qDn VnT , Wps = qDp VT .      τ nωs = τ nω , τ pωs = τ pω KB TL. TLs =. (3.14). KB TL. Then the potential equation (2.1) can be written as .  VT qCm  −δ exp φs + φqns − ϕms us +  − 2 ∆s φs = 2 l εs δ exp −φs − φqns + ϕms vs + Cs 2. which can be rearranged in a simplified form as −λ2 ∆s φs = f1 (φs )   2 −δ exp φs + φqns − ϕms us + . :=  2 δ exp φs + φqns − ϕms vs + Cs. The DG equation (2.2) is transformed to. 1 Cm ∆s ζ ns 2 l√ − Cm ζ ns = VT ln Cm ζ 2ns − VT ln Cm δ 2 us − VT φs , 2b √ n − Cm ζ ns VT 2 = ln ζ ns − ln δ 2 us − φs , 2bn −. 13. (3.15).

(22) which can be written as −ǫ2 ∆s ζ ns = f2 (ζ ns ) := −ζ ns ln ζ 2ns − ln δ 2 us − φs .. (3.16). Similarly, Equation. (2.3) is scaled to −ǫ2 ∆s ζ ps = f3 (ζ ps ) := ζ ps − ln ζ 2ps + ln δ 2 vs − φs .. (3.17). For the electron continuity equation (2.4), we have qDn Cm δ 2 4 (φs )∇s us ∇s · D l2 2 2 4 q Cm neqs peqs − Cm δ exp φqns − φqps us vs =  , 2 δ exp −φs − φqps + ϕms vs  τ Cm  √ +δ 2 exp φs + φqns − ϕms us + 2 neqs peqs exp εktB−εTi −. or. 4 (φs )∇s us = β n f4 (us ) −δ 2 ∇s · D β n neqs peqs − δ 4 exp φqns − φqps us vs : =  (3.18) 2 δ exp −φs − φqps + ϕms vs +   √ 2 εt −εi δ exp φs + φqns − ϕms us + 2 neqs peqs exp kB T where β n =. l2 . Dn τ. Similarly, Equation. (2.5) is scaled to. where β p =. 5 (φs )∇s vs = β p f5 (vs ) := −β p f4 (us ) −δ 2 ∇s · D. (3.19). l2 . Dp τ. For the ET equation (2.6), we obtain qDn Cm δ 2 VT κn Tm qDn VT Cm − ∇s · D (ϕ )∇ g = Jns · Es + ns Wns 6 s ns ns 2 2 l l l2 14.

(23) 6 (ϕns ) κn Tm D −∇s · ρ2n ∇s gns = δ 2 Jns · Es + ns Wns , ρ2n = . qDnVT Cm. (3.20). Similarly, we have for (2.7). −∇s ·. 7 (ϕps ) 2 κp Tm D ρp ∇s gps = δ 2 Jps · Es + ps Wps , ρ2p = . qDp VT Cm. (3.21). Dropping the subscript s in (3.15)-(3.21), the scaled QCET model is simplified to .  2 −δ exp φ + φ − ϕ u+ qn m , −λ2 ∆φ = f1 (φ) =  2 δ exp φ + φqn − ϕm v + C −ǫ2 ∆ζ n = f2 (ζ n ) = −ζ n ln ζ 2n − ln δ 2 u − φ , −ǫ2 ∆ζ p = f3 (ζ p ) = ζ p − ln ζ 2p + ln δ 2 v − φ ,. −∇ · (δ 2n ∇u) = f4 (u). . 4. . . . (3.22) (3.23) (3.24) (3.25). neq peq − δ exp φqn − φqp uv , 2 δ exp −φ − φqp + ϕm v+   √ 2 εt −εi δ exp φ + φqn − ϕm u + 2 neq peq exp kB T. = . . −∇ · (δ 2p ∇v) = −f5 (v) = f4 (u) −∇ · ρ2n ∇gn = f6 (gn ) = δ 2 Jn · E + nWn , −∇ · ρ2p ∇gp = f7 (gp ) = δ 2 Jp · E + pWp ,. with. 15. (3.26) (3.27) (3.28).

(24)  2  VT  λ2 = lε2sqC , ǫ2 = 6m∗qVT l2 , δ 2 = CnmI ,   m n  2 2 2 δ 2n = δ Dβ4 (φ) , δ 2p = −δ βD5 (φ) , β n = Dln τ , β p = n p    7 (ϕ ) 6 (ϕ ) κp Tm D  2 κn Tm D p 2 n  ρn = qDn VT Cm , ρp = qD . p VT Cm. l2 , Dp τ. (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 assumptions (H1)-(H2) below are imposed in order to get simpler expressions for the variables. The non-degeneracy assumption (H3) is valid for semiconductor devices with a doping concentration which is below 1019 cm−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 16.

(25) 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πk 2 / |ε′ (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 nondegenerate Boltzmann statistics. (H4) Let the energy ε(k) satisfy. ε(1 + αε) =. k2 2m∗. (4.2). The constant m∗ is the effective electron mass given by m∗ = m0 kB T0 /2 k02 , where m0 is the unscaled effective mass, k0 is 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 =. T 2β−1 g1 g2 rβ (αT )( − ). µ0 τ nw pβ (αT, 2) pβ (αT, 3). Notice that the function rβ is in fact a polynomial: 17. (4.3).

(26) 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) =. . ∞. us−1 e−u du, s > 0.. (4.5). 0. The assumption (H4) implies γ(T u) = 2m∗ T u(1 + αT u) [14]. Define the functions. pβ (αT, l) = q(αT, l) =. . ∞. 0 ∞. 1 + αT u l−β−1 −u u e du, (1 + 2αT u)2 1. 1. (1 + αT u) 2 (1 + 2αT u)u 2 +l e−u du.. (4.6) (4.7). 0. By letting αT u = x, we have. 1+x 1 1 + αT u = = (1 + x) 2 2 (1 + 2αT u) (1 + 2x) (1 + 2x)2 = (1 + x)(1 − 4x + 12x2 − 32x3 + · · · ). (4.8). = 1 − 3x + 8x2 − 20x3 + · · · and hence. pβ (αT, l) =. . 0. ∞. [1 − 3αT u + 8(αT u)2 − 20(αT u)3 + · · · ]ul−β−1 e−u du (4.9) 18.

(27) Similarly, we have. 1. 1. (4.10) (1 + αT u) 2 (1 + 2αT u) = (1 + x) 2 (1 + 2x) 1/2 · (−1/2) 2 1/2 · (−1/2) · (−3/2) 3 1/2 x+ x + x + ···] = [1 + 1 1·2 1·2·3 (1 + 2x) 1 1 1 = (1 + x − x2 + x3 + · · · )(1 + 2x) 2 8 16 1 2 1 1 3 1 1 = 1 + x − x + x + 2x + x2 − x3 + x4 + · · · 2 8 16 4 8 5 7 2 3 3 = 1 + x + x − x + ··· 2 8 16 and. q(αT, l) =. . ∞ 0. 1 5 7 3 (1 + αT u + (αT u)2 − (αT u)3 + · · · )u 2 +l e−u du (4.11) 2 8 16. 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 (1). (2). g1 (n, T ) = µβ (T )T n, g2 (n, T ) = µβ (T )T 2 n. (4.12). with the temperature-dependent mobilities. (i). µβ (T ) = µ0. pβ (αT, i + 1) −1/2−β T , i = 1, 2. q(αT, 0). (4.13). Here, the mobility constant µ0 is given by. µ0 = (3πφ0 (2N0 + 1)m∗ (2m∗ )3/2 )−1 . 19. (4.14).

(28) 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 g1 g2 r1/2 (0)( − ) µ0 τ nw p1/2 (0, 2) p1/2 (0, 3). (4.15). s denotes after scaling, namely, the term in (3.27), and (1). g1 = µ1/2 (Ts )Ts ns ,. (4.16(a)). (2). g2 = µ1/2 (Ts )Ts2 ns , (1). µ1/2 (Ts ) = µ0. p1/2 (0, 2) −1 (2) p1/2 (0, 3) −1 Ts , µ1/2 (Ts ) = µ0 T , q(0, 0) q(0, 0) s. √ π u e du = p1/2 (0, 2) = q(0, 0) = , 2 0 √ ∞ 3 π 3/2 −u u e du = p1/2 (0, 3) = r1/2 (0) = . 4 0 . Hence, we have. 20. ∞. 1/2 −u. (4.17). (4.18(a)) (4.18(b)).

(29) (1). g1 = µ1/2 (Ts )Ts ns = µ0 Ts−1 Ts ns = µ0 ns ,. (4.16(b)). 3 3 (2) g2 = µ1/2 (Ts )Ts2 ns = µ0 Ts−1 Ts2 ns = µ0 Ts ns , 2 2 and. W1s. W1. √ 1 3 π µ0 n s 3/2µ0 Ts ns √ = (√ − ) µ0 τ nws 4 π/2 3 π/4 3 ns = (TLs − Ts ), 2 τ nws 3 n/Cm (TL /TL − T /TL ) = kB TL , 2 τ nw. (4.19) (4.20). where W1 is unscaled.. The energy-transport equation expressed in [14] is given by. −div J2s = −J1s · ▽s Vs + Wns .. (4.21). Now, we unscale it to obtain. −ldivJ2. l l ▽V 3 n/Cm (TL /Tm − T /Tm ) l2 = −J · l + k T 1 B L qµ0 VT2 Cm qµ0 VT Cm VT 2 τ nw qDn VT 3 n(TL − T ) −divJ2 = −J1 · ▽V + kB 2 τ nw 3 n(T − TL ) −divJ2 = −J1 · ▽V − kB , (4.22) 2 τ nw 21.

(30) where VT = kB TL /q, Dn = un VT and J1 is 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/2kB Tn , our energy relaxation term can be written as. ωn − ω0 3/2kB Tn − 3/2kB TL = −n τ nw τ nw 3 n(Tn − TL ) = − kB 2 τ nw. nWn = −n. (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 Ts g2 g1 r1/2 ( )( − ) Ts µ0 τ nw 2 p1/2 ( 2 , 2) p1/2 ( T2s , 3). (4.24). where (1). g1 = µ1/2 (Ts )Ts ns , (2). g2 = µ1/2 (Ts )Ts2 ns , (1) µ1/2 (Ts ). p1/2 ( T2s , 2) −1 (2) p1/2 ( T2s , 3) −1 = µ0 T , µ1/2 (Ts ) = µ0 T , q( T2s , 0) s q( T2s , 0) s. 22. (4.25).

(31) Ts p1/2 ( , 2) = 2 Ts p1/2 ( , 3) = 2 Ts q( , 0) = 2. . ∞. (1 +. 0. 1 + T2s u 1/2 −u u e du, 2 0 (1 + Ts u) ∞ 1 + T2s u 3/2 −u u e du, 2 0 (1 + Ts u) ∞. Ts 1/2 u) (1 + Ts u)u1/2 e−u du. 2. (4.26(a)) (4.26(b)). (4.27). From (4.9) and (4.11), we have. Ts p1/2 ( , 2) = 2. . Ts p1/2 ( , 3) = 2. . ∞. 3 5 (1 − Ts u + 2Ts2 u2 − Ts3 u3 + · · · )u1/2 e−u du 2 2 ∞ ∞ 0∞ 1 3 5 3 −u −u 2 = u 2 e du − Ts u 2 e du + 2Ts u 2 e−u du 2 0 0 0 5 3 ∞ 7 −u − Ts u 2 e du 2 0 3 3 5 7 5 9 (4.28) = Γ( ) − Ts Γ( ) + 2Ts2 Γ( ) − Ts3 Γ( ), 2 2 2 2 2 2. ∞. 3 5 (1 − Ts u + 2Ts2 u2 − Ts3 u3 + · · · )u3/2 e−u du 2 2 0∞ ∞ ∞ 3 5 7 3 −u −u 2 = u 2 e du − Ts u 2 e du + 2Ts u 2 e−u du 2 0 0 0 ∞ 9 5 3 − Ts u 2 e−u du 2 0 5 3 7 9 5 11 = Γ( ) − Ts Γ( ) + 2Ts2 Γ( ) − Ts3 Γ( ), (4.29) 2 2 2 2 2 2. 23.

(32) . Ts q( , 0) = 2. ∞. 5 7 3 3 3 (1 + Ts u + Ts2 u2 − Ts u + · · · )u1/2 e−u du 4 32 128 0∞ ∞ 1 3 5 7 2 ∞ 5 −u −u −u = u 2 e du + Ts u 2 e du + Ts u 2 e du 4 32 0 0 0 3 3 ∞ 7 −u T u 2 e du − 128 s 0 5 7 3 3 9 3 5 7 T Γ( ), (4.30) = Γ( ) + Ts Γ( ) + Ts2 Γ( ) − 2 4 2 32 2 128 s 2. and hence. p1/2 ( T2s ,2). µ0 Ts ns. p1/2 ( T2s ,3). 1 Ts q( T2s ,0) r1/2 ( )( − ) µ0 τ nws 2 p1/2 ( T2s , 2) p1/2 ( T2s , 3) 1 Ts µ ns µ Ts ns = r1/2 ( )( T0s − 0Ts ) µ0 τ nws 2 q( 2 , 0) q( 2 , 0) ns Ts TLs Ts = r1/2 ( )( Ts − Ts ), τ nws 2 q( 2 , 0) q( 2 , 0) T /TL TL /TL n/Cm T /TL − T /TL ), )( T /TL = r1/2 ( τ nw /kB TL 2 q( , 0) q( , 0). W2s =. W2. µ0 n s. q( T2s ,0). 2. (4.31) (4.32). 2. where, by (4.4),. r1/2 (. Ts 5 5 7 9 1 11 ) = Γ( ) + Ts Γ( ) + 2Ts2 Γ( ) + Ts3 Γ( ). 2 2 2 2 2 2 2. Since. Γ(x + 1) = xΓ(x), for example,. 24. (4.33).

(33) 3 Γ( ) 2 5 Γ( ) 2 7 Γ( ) 2 9 Γ( ) 2. = = = =. √ π , 2√ 3 π , 2 2√ 53 π , 2 2 2√ 753 π , 222 2. we have the formula √ (2k)! π 1 . Γ(k + ) = 2 4k k! The energy-transport equation expressed in [14] is given by. −div J2s = −J1s · ▽s Vs + Wns . Now, we unscale it to obtain. 25. (4.34).

(34) l l ▽V l2 = −J · l + W 1 n qu0 VT2 Cm qµn VT Cm VT qDn VT l l ▽V −ldivJ2 = −J1 ·l 2 qu0 VT Cm qµn VT Cm VT TL /TL n/Cm T /TL T /TL l2 + r1/2 ( )( T /TL − T /TL ) τ nw /kB TL 2 q( 2 , 0) q( 2 , 0) qDn VT l ▽V = −J1 ·l qµn VT Cm VT TL T l2 n T /TL )( T /T − T/T ) + kB r1/2 ( τ nw 2 q( L , 0) q( L , 0) qDn VT Cm −ldivJ2. 2. 2. TL T /TL T )( T /T − T /TL ) L τ nw 2 q( 2 , 0) q( 2 , 0) T /T kB r1/2 ( 2 L ) (TL − T ) T /T q( 2 L ,0) = −J1 · ∇V + n τ nw T /TL kB r1/2 ( 2 ) (T − TL ) T /T q( 2 L ,0) = −J1 · ∇V − n (4.23) τ nw. −divJ2 = −J1 · ∇V +. n. kB r1/2 (. where Dn = un VT and J1 is the electron current density defined as that 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. 26.

(35) re-organized as follows: ∆φ = F (φ, u, v, ζ n , ζ p ), 1 ∇ · Jn = R(φ, u, v, ζ n , ζ p ), q 1 ∇ · Jp = −R(φ, u, v, ζ n , ζ p ), q ∆ζ n = Zn (φ, u, v, ζ n , ζ p ), ∆ζ p = Zp (φ, u, v, ζ n , ζ p ),. (5.1) (5.2) (5.3) (5.4) (5.5). ∇ · Gn = Rn (gn ),. (5.6). ∇ · Gp = Rp (gp ),. (5.7). The main ingredients of the algorithm solving the DGET model are adaptive 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 27.

(36) Step 3.1. Solve the potential equation in (5.1). Step 3.1.1. Set m := 0 and set the initial guess   φ or φ if l = 0, h j j (m) for all (xj , yj ) ∈ Ω , φj = (l)  φ otherwise, j. where φj and φj are constant values that can be easily verified to be h. an upper and lower solution of φ, respectively, and Ω 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. (m+1). Step 3.1.3. Compute φj. by solving the discrete potential system. of (5.1)  (m+1) (m+1) (m)   ξ j φj + γ j (φ) φj = k∈V (j) ξ k φk       −F (φ(m) , u(l) , v (l) , ζ (l) , ζ (l) ) + γ (φ) φ(m) , ∀(xj , yj ) ∈ Ωh , j n p j j j j (m+1)   = VO + Vb , ∀(xj , yj ) ∈ ∂ΩhD , φj    (m+1)    ∂φj = 0, ∀(xj , yj ) ∈ ∂ΩhN , ∂n (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.. 28.

(37) (m). Step 3.1.4. Set φ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. (l+1). Step 3.1.5. Set φ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 j or u j if l = 0, h (m) for all (xj , yj ) ∈ Ω , where u j and u j are uj = (l)  u otherwise, j. h. constant values for all (xj , yj ) ∈ Ω that can be easily verified to. be an upper and lower solution of u, respectively. (m+1). Step 3.2.2. Compute uj. by solving the discrete electron system of. (5.2). (m). Step 3.2.3. Set uj. (m+1). := uj. ∀ j and m := m + 1. Go to Step 3.2.2. until the stopping criteria of the inner iteration are satisfied. (l+1). Step 3.2.4. Set uj. (m+1). := uj. ∀ 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   [ζ ]j or [ζ!]j if l = 0, h n n (m) [ζ n ]j = for all (xj , yj ) ∈ Ω , where [ζ n ]j ≈  [ζ ](l) otherwise, n j. ζ n (xj , yj ) and [ζ n ]j and [ζ! n ]j are constant values for all (xj , yj ) ∈ 29.

(38) h. Ω that can be easily verified to be an upper and lower solution of ζ n , respectively. (m+1). Step 3.4.2. Compute [ζ n ]j (m). Step 3.4.3. Set [ζ n ]j. by solving the discrete system of (5.4).. (m+1). := [ζ n ]j. ∀ j and m := m + 1. Go to Step. 3.4.2 until the stopping criteria of the inner iteration are satisfied. (l+1). Step 3.4.4. Set [ζ n ]j. (m+1). := [ζ n ]j. ∀ j.. Step 3.5. Solve the DG (5.5) similarly to that in Step 3.4. (l+1). Step 3.6. Update [φqn ]j. (l+1). and [φqp ]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.. 30.

(39) Note that, in each one of Steps 3.1-3.5, 4.1, and 4.2, a Jacobi (node-bynode) 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.. 6. Numerical Results. The MOSFET device is shown in Fig. 1. In Fig. 2 and 3 we shows temperature of parabolic band and temperature of non-parabolic band. It can be seen that the maximal temperature value is Tn = 3423 K for α = 0 (parabolic) and Tn = 1260 K for α =. 1 2. (non-parabolic). It shows taht the temperature. of electron is overvalued in parabolic band than that in non-parabolic band. The maximal temperature value is Tm = 1260 K for non-parabolic case and the energy relaxation time τ nw = 0.4 · 10−12 s [7].. 31.

(40) 32.

(41) 33.

(42) 34.

(43) 35.

(44) 36.

(45) 7. Conclusions. In this thesis we have presented the QCET model and its scaled form. The energy relaxation term can be written analytically in terms of the electron density and the temperature for non-parabolic band structure in the sense of Kane. There are two parameters, namely, the non-parabolicity parameter α and the momentum relaxation time parameter β. The variation of the energy relaxation term is smaller in parabolic bands than that in non-parabolic bands. This shows that the QCET model describes the charge flow of electrons in a ballistic diode with reasonable accuracy. The spurious velocity overshoot spike at the anode junction becomes smaller in the non-parabolic band case, compared to the parabolic case [14]. We give explicit formulas of the energy relaxation terms proposed by Degond et al. in [14] for our QCET model. Numerical results have shown that the electron and hole temperatures reduce from about 3500 K for the parabolic case to about 1300 K for. 37.

(46) the non-parabolic case. Therefore, the results we have observed are similar to that of [14].. References [1] M. G. Ancona, G. J. Iafrate, Quantum correction to the equation of state of an electron gas in a semiconductor, Phys. Rev. B 39 (1989) 9536-9540. [2] N. Ben Abdallah and P. Degond, On a hierarchy of macroscopic models for semiconductors, J. Math. Phys. 37 (1996) 3306. [3] N. Ben Abdallah, P. Degond, P. Markowich and C. Schmeiser, High field approximations of the spherical harmonics expansion model for semiconductors, Z. Angew. Math. Phys. 52 (2001), no 2, 201-230. [4] B. A. Biegel, M. G. Ancona, C. S. Rafferty, Z. Yu, Efficient multidimensional simulation of quantum confinement effects in advanced MOS devices, NAS Tech. Report NAS-04-008, 2004. [5] R.-C. Chen, J.-L. Liu, An iterative method for adaptive finite element solutions of an energy transport model of semiconductor devices, J. Comput. Phys. 189 (2003) 579-606. [6] R.-C. Chen, J.-L. Liu, Monotone iterative methods for the adaptive finite element solution of semiconductor equations, J. Comput. Applied Math. 159 (2003) 341-364. 38.

(47) [7] R.-C. Chen, J.-L. Liu, A quantum corrected energy transport model for nanoscale semiconductor devices, J. Comput. Phys. 204 (2005) 131-156. [8] R.-C. Chen, J.-L. Liu, An accelerated monotone iterative method for the quantum-corrected energy transport model, J. Comp. Phys. 227 (2008) 6266-6240. [9] R.-C. Chen, J.-L. Liu, C.T. Lee, A singularly perturbed formulation of the quantum-corrected energy transport model, preprint, (2009). [10] K. C. Chen, Non-Parabolic Band Effects on Energy Relaxation for the Quantum-Corrected Energy Transport Model Part I: Theory, Master Thesis, National University of Kaohsiung 2009. [11] D. Connelly, Z. Yu, D. Yergeau, Macroscopic simulation of quantum mechanical effects in 2-D MOS devices via the density gradient method, IEEE Trans. Electron Devices 49 (2002) 619-626. [12] C. de Falco, E. Gatti, A. L. Lacaita, R. Sacco, Quantum-corrected driftdiffusion models for transport in semiconductor devices, J. Comput. Phys. 204 (2005) 533-561. [13] P. Degond, S. Gallego, F. Méhats, An entropic quantum drift-diffusion model for electron transport in resonant tunneling diodes, J. Comput. Phys. 221 (2007) 226-249.. 39.

(48) [14] P. Degond, A. Jügel, P. Pietra, Numerical discretization of energytransport models for semiconductors with non-parabolic band structure, SIAM on Scientific Computing 22 (2000) 986-1007. [15] P. A. Markowich. The Stationary Semiconductor Device Equations. Springer, Wien, 1986. [16] S. Odanaka, Multidimensional discretization of the stationary quantum drift-diffusion model for ultrasmall MOSFET structures, IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 23 (2004) 837—842. [17] R. Pinnau, Uniform convergence of an exponentially fitted scheme for the quantum drift diffusion model, SIAM J. Numer. Anal. 42 (2004) 1648-1668. [18] E. Pop, S. Sinha, K. E. Goodson, Heat generation and transport in nanometer-scale transistors, Proc. IEEE 94 (2006) 1587-1601. [19] C. S. Rafferty, B. Biegel, Z. Yu, M. G. Ancona, J. Bude, R. W. Dutton, Multi-dimensional quantum effect simulation using a density-gradient model and script-level programming techniques, Proc. SISPAD (1998) 137-140. [20] W. Shockley, Transistor Technology Evokes New Physics, Nobel lecture, December 11, (1956) 344-345.. 40.

(49)