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

全文

(1)國立高雄大學應用數學學系(研究所) 碩士論文. QCET 模型的最大長度之期後誤差估計第一部分: 理論. Maximum Norm A Posteriori Error Estimate for the Quantum-Corrected Energy Transport Model Part 1：Theory. 研究生：林常宇撰指導教授：劉晉良. 中華民國 98 年 7 月.

(2) 論文 etd-0630109-124133 的電子論文授權書. 第 1 頁，共 1 頁. 國立高雄大學博碩士論文電子檔案上網授權書本授權書所授權之論文為授權人在國立高雄大學應用數學系碩士班97學年度第2學期取得碩士學位之論文。論文名稱：指導教授：. QCET 模型的最大長度之期後誤差估計第一部份:理論劉晉良. 茲同意將授權人擁有著作權之上列論文全文﹝含摘要﹞，依下列使用權限之設定，以非專屬、無償授權國立高雄大學圖書館，不限地域、時間與次數，以微縮、光碟或其他各種數位化方式將上列論文重製，並決定是否得將數位化之上列論文及論文電子檔以上載網路方式，提供讀者基於個人非營利性質之線上檢索、閱覽、下載或列印。全文電子檔使用權限：校內外完全公開授權人：親筆簽名或蓋章：. 林常宇民國. 年. 月. 日. 國家圖書館博碩士論文電子檔案上網授權書本授權書所授權之論文為授權人在國立高雄大學應用數學系碩士班97學年度第2學期取得碩士學位之論文。論文名稱：指導教授：. QCET 模型的最大長度之期後誤差估計第一部份:理論劉晉良. 茲同意將授權人擁有著作權之上列論文全文﹝含摘要﹞，以非專屬、無償授權國家圖書館，不限地域、時間與次數，以微縮、光碟或其他各種數位化方式將上列論文重製，並決定是否得將數位化之上列論文及論文電子檔以上載網路方式，提供讀者基於個人非營利性質之線上檢索、閱覽、下載或列印。. 授權人：親筆簽名或蓋章：. 林常宇民國. 年. http://ethesys.nuk.edu.tw/ETD-db/ETD-search-c/authority?URN=etd-0630109-124133. 月. 日. 2009/7/6.

(3)

(4) 致. 謝. 首先我要感謝老師劉晉良教授這兩年來耐心的指導，一開始當我面對這個複雜的模型時，總覺得困惑不解，不知如何下手，經由一次次老師重複的耐心指導，漸漸地越來越能釐清方向。在老師的身上不但學到了專業知識，更能感染到老師開朗的個性和做研究的熱誠，也因此促使我慢慢地體會到做研究的樂趣，面對問題、思考問題、解決問題，享受那份屬於自己的滿足。接著我想要感謝應用數學系這個大家庭，感謝系上所有的老師，感謝系辦兩位親切的大姐姐，讓這個環境充滿了溫馨與互助，感謝這兩年互相打氣的同學，大家一起熬夜，一起狂歡，一起經歷研究生活的點點滴滴，感謝大家最後我要感謝我最愛的父母，這兩年來早出晚歸，好像把家當成了旅館，也鮮少有時間陪父母互動，然而愛我的爸媽卻沒有絲毫的抱怨，默默地支持我，爸媽我永遠愛您們。. 林常宇撰.

(5) Maximum Norm A Posteriori Error Estimate For the Quantum-Corected Energy Transport Model Part I: Theory. By Chang-Yu Lin Advisor Jinn-Liang Liu. Department of Applied Mathematics, National University of Kaohsiung Kaohsiung, Taiwan 811, R.O.C. July 2009.

(6) Contents 1 Introduction. 1. 2 The QCET Model. 3. 3 A Posteriori Error Estimation for a General Semilinear PDE 12 4 Error Estimators for QCET. 27. 5 Conclusion. 36. i.

(7) QCET 模型的最大長度之期後誤差估計第一部份: 模型的最大長度之期後誤差估計第一部份:理論. 指導教授：劉晉良教授國立高雄大學應用數學系. 學生：林常宇國立高雄大學應用數學系. 摘要. 這個量子校正的能量運輸模型是由七條自身伴隨非線性的偏微分方程式所組成，它所探討的是奈米半導體裝置中，電子流和電洞流的穩定性，能量轉換的穩定性，再加上古典電位能和量子電位能的穩定性。我們引用一套期後誤差估計的方法在我們的模型上，這個誤差估計的方法能幫助我們在未來網格點適應性加切時，能當作為誤差指標。關鍵字：量子校正的能量運輸模型、期後誤差估計、適應性加切、誤差指標關鍵字. ii.

(8) Maximum Norm A Posteriori Error Estimate for the Quantum-Corrected Energy Transport Model Part I: Theory. Adviser: Professor Jinn-Liang Liu Institute of Department of Applied Mathematics, National University of Kaohsiung. Student: Chang-Yu Lin Institute 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 develop a second-order maximum norm a posteriori error estimate proposed by Kopteva [11] for the QCET which after scaling involves the scaled Debye length, intrinsic carrier density, Planck constant, and thermal conductivity as the singular perturbation parameters. This estimate can be used as an error indicator for the refinement process in an adaptive algorithm. Keywords: QCET model, posteriori error estimate, error indicator. iii.

(9) 1. Introduction. The model we discuss in this thesis was derived by the physical phenomenon of nano-scale semiconductor device simulation. This thesis discusses a scaling reformulation of the quantum-corrected energy transport model (QCET) presented in [12] into a singularly perturbed system which allows us to investigate the layer solutions of carrier densities, currents, temperatures, and potentials. We have to consider the impact of quantum effects because of the advance of nano-scale devices. The scaled Planck constant is moderately small revealing the importance of quantum effects in nano-devices. In this thesis, we use finite element method to discretize the model. The contribution of this thesis is to apply a new error estimate for the QCET model. This error estimate proposed by Kopteve [11] has only one nonlinear equation with one parameter whereas we discuss seven nonlinear equations with four parameters here. In addition, our boundary condition values are not exactly as the same as [11]. Our model has two boundary condition values which are Dirichlet boundary condition values and Neumann boundary condition values. However, in [11] , the model only has only Dirichlet boundary condition values. Our error estimate looks like below.. 1.

(10) B U − u. . ∞. .        2 2  ≤ C0  max hi M1,ij + max ki M2,ij      i = 0, ..., N  i = 1, ..., N    j = 0, ..., M j = 1, ..., M. (1.1). where, roughly speaking,. M1,ij ≈ Dx2 Uij ln (2 + ε/κ) + 1. M2,ij ≈ Dy2 Uij ln (2 + ε/κ) + 1. (1.2). (1.3). with κ := min{mini {hi }, min{kj }}, ε is the a singularly perturbed parameter. The following is the arrangement of this thesis. In Section 2, we mention the QCET model presented in [12]. Also, the scaling formulation and singular parameters for the model are also given in this section. In Section 3, we divide the seven equations into two types. We can discretize the seven equations by using finite element method and then derive the general form of the matrix algebra system. By applying the theorem of [11] into our seven equations, we can derive a new theorem. In Section 4, we discuss the maximum norm a posteriori error estimate respectively. A short concluding remark is given in Section 5.. 2.

(11) 2. The QCET Model. Based on the density gradient (DG) theory [1], the QCET model proposed in [4] is. −∆φ = f1 (φ) := =. . . φ+φ −nI exp( VTqn )u. q   εs +n exp( −φ−φqp )v + C I VT. (2.1). q (n − p + C) , εs. −∆ζ n = f2 (ζ n ) :=.

(12) −ζ n VT ln(ζ 2n ) − VT ln(nI u) − φ , 2bn. −∆ζ p = f3 (ζ p ) :=.

(13) ζp −VT ln(ζ 2p ) + VT ln(nI v) − φ , 2bp. −∇ · D4 (φ)∇u = f4 (u) :. =. =. . q neq peq − n2I exp(. φqn −φqp )uv VT. −φ−φ. nI exp( VT qp )v + nI exp(  τ √ i +2 neq peq exp( kεtB−ε ) TL. . φ+φqn )u VT. 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 , 3. Wn := −. (2.2). (2.3). (2.4)  . (2.5). ωn − ω0 , τ nω. (2.6).

(14) −∇ · D7 (ϕp )∇gp = f7 (gp ) := Jp ·E + pWp ,. Wp := −. ωp − ω0 , τ pω. (2.7). with the seven unknown functions φ,. u = exp(. −ϕn ), VT. v = exp(. ϕp ), VT. ζn =. √ n,. (2.8). ζp =. √ p,. (2.9). φ + φqn )u, VT. (2.10). −φ − φqp )v, VT. (2.11). n = nI exp(. p = nI exp(. gn = Tn exp(−. gp = Tp exp(. 5ϕn ), 4VT. 5ϕp ), 4VT. (2.12). (2.13). and the auxiliary relations. E = −∇φ. (2.14). φqn = VT ln(ζ 2n ) − VT ln(unI ) − φ,. (2.15). 4.

(15) φqp = −VT ln(ζ 2p ) + VT ln(vnI ) − φ,. D4 (φ) = qDn nI exp(. D5 (φ) = −qDp nI exp(. (2.16). φ + φqn ), VT. (2.17). −φ − φqp ), VT. (2.18). 5ϕn ), 4VT. (2.19). D6 (ϕn ) = κn exp(. D7 (ϕp ) = κp exp(−. 5ϕp ), 4VT. (2.20). Jn = −qµn n∇(φ + φqn ) + qDn ∇n = D4 (φ)∇u,. (2.21). Jp = −qµp p∇(φ + φpn ) − qDn ∇p = D5 (φ)∇v,. (2.22). Gn = D6 (ϕn )∇gn ,. (2.23). Gp = D7 (ϕp )∇gp ,. (2.24). 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. 5.

(16) 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 [16] 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) [4, 16]. 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 [3, 8, 9, 10, 13, 14, 16] by including the energy transport equations (2.6) and (2.7) to deal with hotspot problems in nano-scale device design [15]. Note particularly that the righthand 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. 6.

(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. 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 approxima7.

(18) tion features such as single finite element subspace for all seven variables, global and optimal convergence, and fast iterative solution [5, 7] 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 [14]. Let Ω ⊂

(19) 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. φ = VO + Vb on ∂ΩO ,. ζ 2n. 1 2 2 = C + C + 4nI on ∂ΩO , 2. (2.25). (2.26). ζ p = nI /ζ n on ∂ΩO ,. (2.27). ζ p = ζ n = 0 on ∂ΩI ,. (2.28). u = exp(. −VO ) on ∂ΩO , VT. 8. (2.29).

(20) v = exp(. VO ) on ∂ΩO , VT. (2.30). gn =. 300 on ∂ΩO , O exp( 5V ) 4VT. (2.31). gp =. 300 on ∂ΩO , O exp(− 5V ) 4VT. (2.32). ∂ζ 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.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 [3] [4] across the interface.. 9.

(21) The QCET model is scaled in [12] as .  2 −δ exp φ + φ − ϕ u+ qn m , −λ2 ∆φ = f1 (φ) =  2 δ exp −φ − φqn + ϕm v + C

(22) −2 ∆ζ n = f2 (ζ n ) = −ζ n ln ζ 2n − ln δ 2 u − φ ,.

(23) −2 ∆ζ p = f3 (ζ p ) = ζ p − ln ζ 2p + ln δ 2 v − φ ,. 10. (2.35). (2.36). (2.37).

(24) −∇ · (δ 2n ∇u) = f4 (u). (2.38) neq peq − δ exp φqn − φqp uv =   δ 2 exp −φ − φqp + ϕm v+ .  √ 2 εt −εi δ exp φ + φqn − ϕm u + 2 neq peq exp kB T . 4. −∇ · (δ 2p ∇v) = f5 (v) = −f4 (u). −∇ ·. with. with. . (2.39). ρ2n ∇gn = f6 (gn ) = δ 2 Jn · E + nWn ,. −∇ · ρ2p ∇gp = f7 (gp ) = δ 2 Jp · E + pWp ,.  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      ρ2n = κn Tm D 6 (ϕn ) , ρ2p = κp Tm D7 (ϕp ) . qDn VT Cm qDp VT Cm.  2  VT  λ2 = l2εsqC , 2 = 6m∗qVT l2 , δ 2 = CnmI ,   m n  2 2 2 δ 2n = δ Dβ4 (φ) , δ 2p = −δ βD5 (φ) , β n = Dln τ , β p = n p      ρ2n = κn Tm D 6 (ϕn ) , ρ2p = κp Tm D7 (ϕp ) . qDn VT Cm qDp VT Cm. (2.40). (2.41). l2 , Dp τ. (2.42). l2 , Dp τ. (2.43). where λ is called the scaled Debye length, the scaled Planck constant, δ the standard intrinsic carrier density, δ n, p the modified intrinsic electron and 11.

(25) hole densities, and ρn, p the electron and hole thermal conductivities. The order of these parameters is summarized in Table 2 for which the applied voltages of the device are VBS = 0 V, VDS = 1 V, and VGS = 0.8 V [12]. Table 2. Order of parameters λ2 = O(10−3 ) 2 = O(10−3 ) δ 2 = O(10−8 ) O(10−27 ) ≤ δ 2n ≤ O(105 ) O(10−14 ) ≤ δ 2p ≤ O(1018 ) O(10−32 ) ≤ ρ2n ≤ O(102 ) O(10−19 ) ≤ ρ2p ≤ O(1018 ) The boundary conditions (2.21)-(2.30) can be similarly scaled.. 3. A Posteriori Error Estimation for a General Semilinear PDE. The signularly perturbed semilinear reaction-diffusion problem considered in [11] read as T u := − (ε u) + b(x, y, u) = 0, (x, y) ∈ Ω = (0, 1) × (0, 1),. u(x, y) = gD , (x, y) ∈ ∂ΩD ,. 12. (3.1). (3.2).

(26) ∂u(x, y) = gN , (x, y) ∈ ∂ΩN , ∂n. (3.3). where the boundary of a bounded open set domain Ω is denoted as ∂Ω, ε is a small positive parameter, and b(x, y, u) is nonlinear and has an unique solution in u. Part of the boundary of Ω will satisfy equation (3.2) while the rest of them will satisfy equation (3.3). The Dirichlet function gD is the well-known assumption of the boundary condition number for Poisson’s equation. n in equation ( 3.3) is an outward normal vector. The parameter ε in equation (3.1) will be considered in two cases. For the first case, ε is a constant while ε is a function of (x,y) on Ω for the second case. Note that in [11] Kopteva does not consider the Neumann boundary condition and the parameter ε is a constant. For our QCET model, this is too restrict. Our presentation of (3.1) is more general. We discrete the equations (3.1) - (3.3) by finite element method. The discussion falls into two parts. We will show the weak problem and finite element problem of equations (3.1) - (3.3).. 3.1 3.1.1. Weak Problem Case 1: ε = constant. Since ε is constant, we can rewrite (3.1) as the following. −εu + b(x, y, u) = 0, on Ω,. 13. (3.5).

(27) u(x, y) = gD. on ∂ΩD ,. (3.6). on ∂ΩN. (3.7). ∂u(x, y) = gN ∂n. where = ∂/∂x2 + ∂/∂y 2 is the Laplace operator. From the strong form of (3.5), we know. −ε. . Ω. (u)vdΩ = −. . bvdΩ. (3.8). Ω. is true for arbitrary test function v(x, y). Integration by parts gives ε − ∇ (v∇u) dΩ + ∇u∇vdΩ = − bvdΩ Ω. Ω. (3.9). Ω. Green’s theorem applied to the first item of (3.9) gives. −ε. . ∂Ω. (v∇u) · nds + ε. . Ω. ∇u∇vdΩ = −. . bvdΩ. (3.10). Ω. Hence we can obtain the following weak formulation of (3.5):. ε. . Ω. ∇u∇vdΩ = ε. . ∂ΩD. (v∇u) · nds +. . ∂ΩN. (v∇u) · nds − bvdΩ (3.11) Ω. Therefore, we can obtain the weak problem by organizing the strong problem of (3.5)-(3.7) as follow. Weak Problem: Assume nonlinear function b in u, u is unique solution. Then given the constants gD and gN , the weak formulation form is. 14.

(28) B(u, v) = F (u, v) ∀v ∈ H01 (Ω). (3.12). H 0 (Ω) = L2 (Ω). (3.13). 2

(29) 2 H (Ω) : = v(x) : v + (∇v) dΩ < ∞. (3.14). where. 1. Ω. H01 (Ω) : = v(x) ∈ H 1 (Ω) : v = 0 on ∂ΩD 1 HD (Ω) : = v(x) ∈ H 1 (Ω) : v = gD on ∂ΩD B(u, v) = ε. F (u, v) = ε. . Ω. ∂ΩN. . ∂ΩD. 3.1.2. u • vdΩ. gN vds −. . bvdΩ. (3.15). (3.16). (3.17). (3.18). Ω. (v∇u) · nds = 0. (1). Case 2: ε = ε(x, y). In this case, the first item of (3.1) is not a constant. We derive the weak form (3.1) in the same way as that for case 1. However, we only need to focus on the integration part of the Poisson equation. The conditions of boundary. 15.

(30) condition number of (3.2) and (3.3) are the same as Case1. Firstly, we multiply the equation (3.1) by the test function v.. −. . (ε(x, y) u) vdΩ +. Ω. . bvdΩ = 0. (3.19). Ω. Integration by parts gives − ∇ [(ε(x, y)∇u) v] dΩ + (ε(x, y)∇u) · ∇vdΩ = − bvdΩ Ω. Ω. (3.20). Ω. Green’s theorem applied to the first item of (3.9) gives. −. . ∂Ω. ((ε(x, y)∇u) v) · nds +. . Ω. ε(x, y)∇u · ∇vdΩ = −. . bvdΩ. (3.21). Ω. Hence we can obtain the following weak formulation of (3.21) . Ω. =. . ε(x, y)∇u · ∇vdΩ ∂ΩD. (3.22). ((ε(x, y)∇u) v) · nds +. . ∂ΩN. ((ε(x, y)∇u) v) · nds − bvdΩ Ω. Therefore, we can obtain the weak problem of case 2 as follows. Weak Problem: Assume nonlinear function b in u, u is an unique solution. Then given the constants gD and gN , the weak formulation form is. B(u, v) = F (u, v) ∀v ∈ H01 (Ω). B(u, v) =. . Ω. ε(x, y)∇u · ∇vdΩ 16. (3.23). (3.24).

(31) F (u, v) =. . ∂ΩN. 3.2 3.2.1. ε(x, y)gN vds −. . bvdΩ. (3.25). Ω. Finite Element Problem Case 1: ε = constant. We denote the linear finite element subspace. 1 S h := span{φi }N i=1 ⊂ H (Ω). (3.26). where φi are basis functions. Then we can approximate u, v by uh , v h respectively. We rewrite the equation (3.23) as the following finite element problem. Finite Element problem: Assume nonlinear function b in u, u is an h unique solution. Then given the constants gD and gN , find uh ∈ SD such. that. B(uh, v h ) = F (uh , vh ) ∀v h ∈ S0h ⊂ H01 (Ω). (3.27). where. h. h. B(u , v ) =. . ε uh • v h dΩ. (3.28). Ω. h 1 SD := {v h(x) ∈ S h : v = gD on ∂ΩD } ⊂ HD (Ω). uh =. 17. Uj φj ,. (3.29). (3.30).

(32) where Ui are unknown scalars and φi (x, y) is the basis function on the nodal point i. The FE approximation of problem (3.5)-(3.7) can be expressed as. u(x, y) ≈ uh(x, y) =. N . Ui φi (x, y),. (3.31). i=1. The test function v(x, y) can also be expressed as the linear combination of the shape functions, i.e.,. v(x, y) ≈ v h (x, y) =. M . Vj φj (x, y), Vj : arbitrary scalars. (3.32). j=1. Substituting (3.31) and (3.32) into (3.12), we have . ε uh • . Ω. =. . ε. ∂ΩN. N j=1. j=1. =. Vj. j=1. . ∂ΩN. . =. Ω. Vj. . Vj φj dΩ. N. Vj φj (x, y)bdΩ. (3.34). Ω j=1. ε uh • φj dΩ. (3.35). Ω. εφj (x, y)gN ds −. N j=1. Vj. . φj (x, y)bdΩ. εφj (x, y)gN ds − 18. (3.36). Ω. ε uh • φj dΩ. ∂ΩN. (3.33). j=1. Vj φj (x, y)gN ds −. N N . N . (3.37) . Ω. φj (x, y)bdΩ ∀j.

(33) B(uh, φj ) = F (uh, φj ), ∀j. (3.38). Substituting (3.31) into (3.38), we obtain B(. N . N . Ui φi , φj ), ∀j. (3.39). N . Ui φi , φj ), ∀j. (3.40). Ui φi , φj ) = F (. i=1. N . i=1. Ui B(φi , φj ) = F (. i=1. i=1. .     B(φ1 , φj ) B(φ2 , φj ) · · · B(φn , φj )    . . U1   N U2   Ui φi , φj ), ∀j = F (  ..  .  i=1  Un (3.41). Since. B(φi , φj ) =. . ε φi • φj dΩ = B(φj , φi ), ∀j,. (3.42). Ω. the matrix formulation of (3.41) can be written as. AU = F (U ). (3.43). . (3.44). where A is symmetric. aij = B(φi , φj ) =. Ω. 19. ε φi • φj dΩ.

(34) N . [F (U )]j = F (. Ui φi , φj ) =. i=1. . ∂ΩN. εφj (x, y)gN ds −. . φj (x, y)bdΩ (3.45). Ω. Note that F (U ) is a nonlinear vector since b is nonlinear in U . 3.2.2. Case 2: ε = ε(x, y). For this case, all the equations in the above equations (3.27)-(3.45) remain the same except that the parameter ε is now a function instead of a constant.. 3.3. Maximum norm a posteriori error estimate. Theorem 3.1.. Let u(x, y) be a solution of the problem (3.1)-(3.3) there _. _. exist constants β and β such that 0<β <bu < β, U be an uniquely solution of the problem (3.43) on an arbitrary mesh {(xi , yj )}, and U B (x, y) be the _. piecewise bilinear interpolation of U. i.e., U B = U B (x, y) is continuous in Ω, bilinear on each [xi−1 , xi ] × [yj−1 , yj ] , and equal to Uij at the mesh nodes. Then. B U − u. . ∞. .        2 2  max hi M1,ij + max ki M2,ij  ≤C     i = 0, ..., N  i = 1, ..., N    j = 0, ..., M j = 1, ..., M (3.46). 20.

(35) where. . 2 . M1,ij := min Dx2 Ui−1,j , Dx2 Uij ln(2 + ε/κ) + ε Dx− Dx2 Uij + Dx− Uij + 1, M2,ij. (3.47). 2 . − 2. − 2 2 := min Dy Ui,j−1 , Dy Uij ln(2 + ε/κ) + ε Dy Dy Uij + Dy Uij + 1, (3.48). with κ := min {mini {hi } , minj {kj }} , while the constant C is independent of ε and the mesh. Proof: We imitate Theorem 2.1 in [11] whereas we make boundary conditions different. First of all, we substitute U B into (3.1) and define a discrete function Vij by. Dx− Vij =. Dx2 Vij. Dx− Vi+1,j − Dx− Vi,j = , (hi + hi+1 )/2. Dy− Vij =. Dy2 Vij. Vij − Vi−1,j , hi. Vij − Vi,j−1 , kj. Dy− Vi,j+1 − Dy− Vi,j = . (kj + kj+1 )/2. (3.49). (3.50). (3.51). (3.52). We thus have −ε2 Dx2 UijB − ε2 Dy2 UijB + b xi , yj , UijB = 0,. (3.53). for i = 1, ...N − 1, j = 1, ...M − 1, where U B = gD on ∂ΩD , ∂U B /∂n = gN on ∂ΩN . We now extend Dx2 Uij to the mesh node i = 0, N as follows. Firstly, we 21.

(36) formally enlarge the discrete equation (3.53) to i = 0 and i = N by using the Neumann boundary condition. Since ∂U B /∂n = gN which is ∂U B /∂x = gN , then ∂ 2 U B /∂x = 0. This yields. Dy2 U0,j := ε−2 b(0, yj , U0,j ). (3.54). Dy2 UN,j := ε−2 b(1, yj , UN,j ). (3.55). j = 0, ..., M . Similarly, we use the Dirichlet boundary condition on the mesh node j = 0, M . We consider that Ui,0 =Ui,M =gD , and Dx2 Uij = 0. Similarly, it also yields Dy2 Ui,0 := ε−2 b(xi , 0, Ui,0 ). (3.56). Dy2 Ui,M := ε−2 b(xi , 1, Uj,M ). (3.57). Now that D2x Uij and D2y Uij are extended to all i, j. Our discrete equation (3.53) holds true for all i = 0, ..., N and j = 0, ..., M. By (3.1), we have

(37). T U B − T u = −ε2 ∂ 2 /∂x2 + ∂ 2 /∂y 2 U B + b(x, y, U B ). (3.58). where ∂U B /∂x2 and ∂U B /∂y 2 can be understood in the sense of distributions. We then define an auxiliary function. q(x, y) := b(x, y, U B (x, y)). (3.56). and let q B denote its piecewise bilinear interpolant on {(xi, yj )} . Hence. 22.

(38)

(39).

(40) T U B − T u = −ε2 ∂ 2 /∂x2 + ∂ 2 /∂y 2 U B + q B + q − q B .. (3.57). By (3.53). qij = b (xi , yj , Uij ) = ε2 Dx2 Uij + ε2 Dy2 Uij. (3.58). i = 0, ..., N ; j = 0, ..., M Then q B (x, y) = bB (xi , yj , Uij ) = ε2 Dx2 UijB + ε2 Dy2 UijB _.

(41) J

(42) I = ε2 Dx2 UijI + ε2 Dy2 UijJ , (x, y) ∈ Ω. (3.59). Therefore B. TU − Tu =. . −ε. 2∂. 2. J. U I (x, yj ) + q1I (x, yj ) ∂x2. +. (3.60). I.

(43) U J (xi , y) I + q2 (xi , y) + q − qB −ε 2 ∂y.

(44) ∂ ∂ = F1 (x, y) + F2 (x, y) + q − q B ∂x ∂y 2∂. 2. where F1 and F2 are functions of the current mesh and computed solution. This will enable us to estimate the error U B − u in the maximum norm by linearizing the operator T and invoking its stability properties. (See [11] for more details) We have. 23.

(45) . B U − u ∞. "  |q1,ij |    i = 1, ..., N         j = 0, ..., M  ε  ! (3.61) ≤  " ln 2 +   kj2 κ  +  max |q2,ij |  ε2    i = 0, ..., N     j = 1, ..., M  ! 2 "  hi − max |Dx q1,ij | ε2    i = 1, ..., N        j = 0, ..., M    ! k2. "  +C . j. −.   +  max 2 Dy q2,ij ε     i = 0, ..., N     j = 1, ..., M +β −1 q − qB !. max. ∞. 24. h2i ε2.

(46) by Theorem 4.1 in [11] yields . !. max. h2i ε2. |q1,ij |. " .    i = 1, ..., N        B  j = 0, ..., M  U − u ≤ C   ln(2 + ε ) ! " 2   kj κ  + |q2,ij |  max ε2     i = 0, ..., N     j = 1, ..., M  ! 2 "  hi − max |Dx q1,ij | ε     i = 1, ..., N       j = 0, ..., M   "  ! k2. +C .  j. −.    + D q max 2,ij y ε     i = 0, ..., N     j = 1, ..., M +β −1 q − qB . (3.62). ∞. and. . q − q B ∞.       ≤C      . ! "  h2i 1 + |Dx− Uij |2   i = 1, ..., N    j = 0, ..., M  ! " . . 2 2 −.  + max kj 1 + Dy Uij   i = 0, ..., N   j = 1, ..., M max. (3.63). So we obtain a version of the desired a posteriori error estimate (3.46) in. . which the quantities min{|Dx2 Ui−1,j | , |Dx2 Uij |} and min Dy2 Ui,j−1 , Dy2 Uij. are replaced by |Dx2 Uij | and Dy2 Uij , respectively. The quantities |q1,ij | and 25.

(47) |q2,ij | in (3.61) can be replaced by min{|q1,i−1,j | , |q1,ij |} and min{|q2,i,j−1 | , |q2,ij |}, respectively. Combining this sharper version of (3.61) with (3.58) and (3.62) yields the desired estimate (3.18).. . Theorem 3.2. Let u(x, y) be a solution of the problem (3.1)-(3.3), and bu ≤ 0 on some interval, and Uij a solution of problem (3.43) on an arbitrary mesh {(xi , yj )}, and U B (x, y) be the piecewise bilinear interpolant of U . i.e., _. U B = U B (x, y) is continuous in Ω, bilinear on each [xi−1 , xi ] × [yj−1 , yj ] , and equal to Uij at the mesh nodes. Then. B U − u. . ∞. .        2 2   max hi M1,ij + max ki M2,ij  ≤C    i = 0, ..., N   i = 1, ..., N   j = 0, ..., M j = 1, ..., M. (3.63). Proof: We take bu =.  . bu. if bu > 0.  max b if u. and then we prove as Theorem 3.1. .  . bu ≤ 0 . ,. Remark: The right hand sides of (3.46) and (3.63) give us precise definitions of a posteriori error indicators for all elements and the error estimator of the computed solution. Since the terms M1,ij , and M2,ij involve the computed solution Uij at each mesh point (xi , yj ), and the singularly perturbed parameters ε which is know. 26.

(48) 4. Error Estimators for QCET. We now apply the previous two theorem to define the error estimator define by η for each one of the seven PDEs (2.35)-(2.41). Note particularity that the notations U B , u, b, and ε represent differently for each PDE of the QCET model. However the definition of κ in (3.46) is the same through out this thesis.. 4.1. Potential Equation (2.35)  2 −δ exp φ + φ − ϕ u+ qn m  −λ2 ∆φ = f1 (φ) =  2 δ exp −φ − φqn − ϕm v + C . −λ2 ∆φ − f1 (φ) = 0. (4.1). (4.2). −f1 (φ) −f1 (φ) = (4.3) ∂φ ∂φ.

(49) = δ 2 exp φ + φqn − ϕm u + exp −φ − φqn + ϕm v > 0 Corresponding to (3.5), we ε = λ, b = −f1 (φ), u = φ. The discrete system of nonlinear algebraic equations (3.43) is expressed as. A1 U1 − F (U1 ) = 0. (4.4). So we know (4.2), (4.3) satisfy Theorem 3.1. The maximum norm a posteriori error estimate of Theorem 3.1. implies that. 27.

(50) ˜ e := U1B − φ ≤ Cη,. (4.5). η := max {η 0 , η 1 , η2 , η 3 }. ηl := max.          . max.    i = 1, ..., N       j = 0, ..., M. !. " (l) h2i M1,ij ;. max i = 0, ..., N. (4.6). ! (l) ki2 M2,ij. j = 1, ..., M. ˜. (0) (0) where l = 0, 1, 2, 3, C = C ln(2 + λ/κ), while M1,ij = M2,ij = 1,. 2 (1) M1,ij : = Dx− Uij ,. 2 (1) M2,ij : = Dy− Uij ,. . (2) M1,ij := min Dx2 Ui−1,j , Dx2 Uij ,. . (2) M2,ij := min Dy2 Ui,j−1 , Dy2 Uij ,. (3) M1,ij := ε Dx− Dx2 Uij ,. (3) M2,ij := ε Dy− Dy2 Uij , Dx− Uij =. Uij − Ui−1,j , hj 28.          "         . ,. (4.7). (4.8) (4.9) (4.10) (4.11). (4.12). (4.13). (4.14).

(51) Dy− Uij =. Uij − Ui,j−1 . kj. (4.15). Here η l and M (l) , involve discrete analogues of lth-order derivatives and. Dx2 Uij =. 2 [hi+1 Ui−1,j − (hi + hi+1 )Uij + hi Ui+1,j ] , (hi + hi+1 ) hi+1 hi. (4.16). Dy2 Uij =. 2 [kj+1 Ui,j−1 − (kj + kj+1 )Uij + kj Ui,j+1 ] , (kj + kj+1 ) kj+1 kj. (4.17). Dx2 Ui−1,j =. 2 [hi Ui−2,j − (hi−1 + hi )Ui−1,j + hi−1 Ui,j ] , (hi−1 + hi ) hi−1 hi. (4.18). Dy2 Ui,j−1 =. 2 [kj Ui,j−2 − (kj−1 + kj )Uij−1 + kj−1 Ui,j ] (kj−1 + kj ) kj−1 kj. (4.19). 2 Ui−2,j hi−1 hi (hi + hi+1 ) & ' 1 1 + Ui−1,j +2 hi−1 hi (hi + hi+1 ) hi+1 & ' 1 1 1 + Ui,j + Ui+1,j − hi+1 hi−1 + hi hi + hi+1. Dx− Dx2 Uij = −. 4.2. (4.20). Density Gradient Equation (2.36). Equation (2.36) can be written as. −2 ∆ζ n − f2 (ζ n ) = 0,

(52) f2 (ζ n ) = −ζ n ln ζ 2n − ln δ 2 u − φ 29. (4.21). (4.22).

(53) A2 U2 = F2 (U2 ).. (4.23). We know that the partial derivative of f2 with respect to ζ n is not all positive on the domain (see numerical results present in [7] ), so we know (4.21) satisfy Theorem 3.2. The maximum norm a posteriori error estimate of Theorem 3.2 implies that ˜ e := U2B − φ ≤ Cη,. where. (4.24). η := max {η 0 , η 1 , η2 , η 3 }. η l := max.          . max.    i = 1, ..., N       j = 0, ..., M. ! " (l) h2i M1,ij ;. max i = 0, ..., N. (4.25). ! (l) ki2 M2,ij. j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + /κ).. 4.3.          "         . , (4.26). Density Gradient Equation (2.37). Again, we have. −2 ∆ζ p − f3 (ζ p ) = 0,

(54) f3 (ζ p ) = −ζ p ln ζ 2p − ln δ 2 v + φ 30. (4.27). (4.28).

(55) A3 U3 = F3 (U3 ).. (4.29). We know that the partial derivative of f3 with respect to ζ p is not all positive on the domain, so we know (4.37) satisfy Theorem3.2. The maximum norm a posteriori error estimate of Theorem 3.2 implies that ˜ e := U3B − ζ p ≤ Cη,. where. (4.30). η := max {η 0 , η 1 , η2 , η 3 }. η l := max.          . max.    i = 1, ..., N       j = 0, ..., M. ! " (l) h2i M1,ij ;. max i = 0, ..., N. (4.31). ! (l) ki2 M2,ij. j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + /κ).. 4.4.          "         . , (4.33). Electron Continuity Equation (2.38). Equation (2.38) is rewritten as follow. −∇ · (δ 2n ∇u) − f4 (u) = 0,. 31. (4.34).

(56) f4 (u) = β n S δ 2 exp φ + φqn − ϕm u, δ 2 exp −φ − φqp + ϕm v, neq peq. (4.35). The finite element method and the Scharfetter-Gummel exponential fitting scheme of [5] are used to approximate the boundary value problem (2.37) and yield a system of nonlinear algebraic equations in matrix form as A4 U4 = F4 (U4 ).. (4.36). We know that the partial derivative of −f4 with respect to ζ p is all positive on the domain, so we know (4.48) satisfy Theorem3.1. The maximum norm a posteriori error estimate of Theorem 3.1 implies that ˜ e := U4B − u ≤ Cη,. where. (4.37). η := max {η 0 , η 1 , η2 , η 3 }. η l := max.          . max.    i = 1, ..., N       j = 0, ..., M. ! " (l) h2i M1,ij ;. max i = 0, ..., N j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + δ n /κ). 32. (4.38). ! (l) ki2 M2,ij.          "         . , (4.39).

(57) 4.5. Hole Continuity Equation (2.39). Similarly, we have. −∇ · (δ 2p ∇v) + f5 (v) = 0,. (4.40). f5 (v) = −β p S δ 2 exp φ + φqn − ϕm u, δ 2 exp −φ − φqp + ϕm v, neq peq. (4.41). The finite element method and the Scharfetter-Gummel exponential fitting scheme of [5] are used to approximate the boundary value problem (2.39) and yield a system of nonlinear algebraic equations in matrix form as A5 U5 = F5 (U5 ).. (4.42). We know that the partial derivative of f5 with respect to v is all positive on the domain, so we know (4.54) satisfy Theorem3.1. The maximum norm a posteriori error estimate of Theorem 3.1 implies that. where. ˜ e := U4B − u ≤ Cη,. η := max {η 0 , η 1 , η2 , η 3 }. 33. (4.43). (4.44).

(58) η l := max.          . ! " (l) h2i M1,ij ;. max.    i = 1, ..., N       j = 0, ..., M. max i = 0, ..., N. ! (l) ki2 M2,ij. j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + δ p /κ).. 4.6.          "         . , (4.45). Energy Transport Equation (2.40). Similarly, we have −∇ ·. . ρ2n ∇gn − f6 (gn ) = 0,. f6 (gn ) = δ 2 Jn · E + nWn. (4.46). (4.47). The finite element method and the Scharfetter-Gummel exponential fitting scheme of [5] are used to approximate the boundary value problem (2.40) and yield a system of nonlinear algebraic equations in matrix form as A6 U6 = F6 (U6 ).. (4.48). We know that the partial derivative of f6 with respect to gn is all positive on the domain, so we know (4.60) satisfy Theorem 3.1. The maximum norm a posteriori error estimate of Theorem 3.1 implies that ˜ e := U4B − u ≤ Cη,. 34. (4.49).

(59) where. η := max {η 0 , η 1 , η2 , η 3 }. η l := max.          . max.    i = 1, ..., N       j = 0, ..., M. ! " (l) h2i M1,ij ;. max i = 0, ..., N j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + ρn /κ).. 4.7. (4.50). ! (l) ki2 M2,ij.          "         . , (4.51). Energy Transport Equation (2.41). Finally, we have −∇ · ρ2p ∇gp − f7 (gp ) = 0, f7 (gp ) = δ 2 Jp · E + pWp ,. (4.66). (4.67). The finite element method and the Scharfetter-Gummel exponential fitting scheme of [5] are used to approximate the boundary value problem (2.41) and yield a system of nonlinear algebraic equations in matrix form as A7 U7 = F7 (U7 ).. 35. (4.68).

(60) We know that the partial derivative of f7 with respect to gp is all positive on the domain, so we know (4.66) satisfy Theorem3.1. The maximum norm a posteriori error estimate of Theorem 3.1. implies that ˜ e := U4B − u ≤ Cη,. where. (4.69). η := max {η 0 , η 1 , η2 , η 3 }. η l := max.          . max.    i = 1, ..., N       j = 0, ..., M. ! " (l) h2i M1,ij ;. i = 0, ..., N j = 1, ..., M. ˜. l = 0, 1, 2, 3, C = C ln(2 + ρp /κ).. 5. max. (4.70). ! (l) ki2 M2,ij.          "         . , (4.71). Conclusion. We have extended the a posteriori error theory developed by Kopteva [11] to our singularly perturbed quantum-corrected energy transport model which consists of seven semilinear PDEs with the scaled Debye length, intrinsic carrier density, Planck constant, and Thermal conductivity as the singular perturbation parameters. For these seven equations, we present explicit formulas for computing the error indicators and estimates, which are essential 36.

(61) elements are indispensable in the state-of-the-art mathematical software for engineering and scientific computing including the semiconductor device simulation for advanced nano-devices in industry.. 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] M.G. Ancona, Z. Yu, R.W. Dutton, P.J.V. Voorde, M. Cao, D. Vook, Density-gradient analysis of MOS tunneling, IEEE Trans. Electron. Dev. 47 (2000) 2310. [3] 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. [4] R.-C. Chen, J.-L. Liu, A quantum corrected energy transport model for nano-scale semiconductor devices, J. Comput. Phys. 204 (2005) 1347 (2000) 2310. [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.. 37.

(62) [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. [7] 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. [8] 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. [9] 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. [10] 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. [11] N. Kopteva, Maximum norm a posteriori error estimate for a 2d singularly perturbed reaction-diffusion problem, SIAM J. Numer. Anal., 46 (2008), 1602-1618. [12] J.-L. Liu R.-C. Chen C.-T. Lee. A singular perturbed formulation of the quantum-corrected energy transport model, preprint, 2009.. 38.

(63) [13] 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. [14] R. Pinnau, Uniform convergence of an exponentially fitted scheme for the quantum drift diffusion model, SIAM J. Numer. Anal. 42 (2004) 1648-1668. [15] E. Pop, S. Sinha, K. E. Goodson, Heat generation and transport in nanometer-scale transistors, Proc. IEEE 94 (2006) 1587-1601. [16] 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.. 39.

(64)