半導體奈米結構之數值模擬(2/2)
計畫類別： 個別型計畫
計畫編號： NSC92-2115-M-009-008-
執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日
執行單位： 國立交通大學應用數學系(所)
計畫主持人： 劉晉良
報告類型： 完整報告
處理方式： 本計畫可公開查詢
中 華 民 國 93 年 12 月 6 日
UNCORR
ECTE
D P
ROO
F
2
A quantum corrected energy-transport model for
3
nanoscale semiconductor devices
4
Ren-Chuen Chen, Jinn-Liang Liu
*5 Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan Received 29 June 2004; received in revised form 4 October 2004; accepted 4 October 2004
8 Abstract
9 An energy transport model coupled with the density gradient method as quantum mechanical corrections is
pro-10 posed and numerically investigated. This new model is comprehensive in both physical and mathematical aspects. It
11 is capable of describing hot electron transport as well as significant quantum mechanical effects for advanced devices
12 with dimensions comparable to the de Broglie wave-length. The model is completely self-adjoint for all state variables
13 and hence provides many appealing mathematical features such as global convergence, fast iterative solution, and
14 highly parallelizable. Numerical simulations on diode and MOSFET with the gate length down to 34 nm using this
15 model have been performed and compared with that using the classical transport model. It is shown that the I–V
char-16 acteristics of this short-channel device is significantly corrected by the density-gradient equations with current drive
17 reduced by up to 60% comparing with that of the classical model along. Moreover, a 2D quantum layer, which is only
18 a fraction of the length scale of inversion layer, is also effectively captured by this new model with very fine mesh near
19 the interface produced by an adaptive finite element method.
20 2004 Published by Elsevier Inc.
21
22 1. Introduction
23 Numerical simulation of charge transport in device structures is widely used for analysis of physical 24 processes in the semiconductor devices and estimation of their electrical parameters. The major part of 25 the activities in this field is based on drift-diffusion (DD) equations. However, there is a growing realization 26 that technologist cannot ignore quantum effects much longer. The combination of thin gate oxides and hea-27 vy doping in the conventional MOSFETs, and the thin silicon body of the double-gate structures, will result 28 in substantial quantum mechanical (QM) threshold voltage shift and transconductance degradation[21].
0021-9991/$ - see front matter 2004 Published by Elsevier Inc. doi:10.1016/j.jcp.2004.10.006
*
Corresponding author.
E-mail address:[email protected](J.-L. Liu).
Journal of Computational Physics xxx (2004) xxx–xxx
UNCORR
ECTED
PROO
F
29 Computationally efficient methods to include QM effects are required for the purpose of practical Compu-30 ter Aided Design of this generation of devices.
31 Some numerical methods employing full quantum models such as non-equilibrium GreenÕs function[23]
32 and WignerÕs function[9]suffer from unsolved robustness problems and are still much too costly for device 33 or circuit simulations. Another approach for including QM effects is to add quantum corrections to clas-34 sical models[4,3,15–20,22,28,30,32,33]. In particular, the density gradient (DG) model developed by Anc-35 ona et al. is a more rigorous macroscopic transport model which avoids ad hoc assumption to the material 36 parameters or imposing an artificial shape function[34]. It is demonstrated in[1,7,8,29]that this model is 37 feasible and efficient to accurately and generally simulate multi-dimensional devices with gate lengths rang-38 ing from 30 nm down to 6 nm when combined with the DD model.
39 In this paper, we further extend the DG model to combine with the energy transport (ET) model pro-40 posed in our previous work[11]and show that this new combination (DGET) is capable of describing hot 41 electron transport as well as significant QM effects for advanced devices. Our model is able to explain that 42 electron temperature essentially differs from the lattice temperature. It is clear that this effect cannot be de-43 scribed by the DG model along. Quantum hydrodynamic (QHD) models give accurate simulation results, 44 but the numerical methods to solve this system are too costly and time consuming to model real problems in 45 semiconductor production mode where simulation results are needed in hours or minutes. The DGET mod-46 el is of parabolic type so that its numerical solution needs less effort than QHD models which contain 47 hyperbolic modes.
48 Moreover, our model is completely self-adjoint for all state variables and hence provides many appealing 49 mathematical features such as global convergence, fast iterative solution, and highly parallelizable as dem-50 onstrated in our previous papers[11,12,24]. The global convergence is a consequence of monotone iterative 51 methods used in solving the discrete systems of nonlinear algebraic equations resulting from adaptive finite 52 element approximation of the model. It is shown here that the convergence analysis of these methods given 53 in[11,12]can be straightforwardly carried over to the present model. Our numerical experiments on various 54 device structures with high drain bias have shown that the monotone iteration do not suffer from the con-55 vergence difficulties as frequently encountered by the commonly used NewtonÕs iteration since the Jacobian 56 is either close to singular or poorly conditioned[29]. This is a fundamental issue constantly faced by the 57 practitioner in device and circuit modeling. Numerical simulations on diode and MOSFET with the gate 58 length down to 34 nm using the DGET model have been performed and compared with that using the 59 ET model. It is shown that the I–V characteristics of this short-channel device is significantly corrected 60 by the density-gradient equations with current drive reduced by up to 60% compared with that of the clas-61 sical model along. Moreover, a 2D quantum layer, which is only a fraction of the length scale of inversion 62 layer, is also effectively captured by this new model with very fine mesh near the interface produced by the 63 adaptive finite element method.
64 The paper is divided into the following sections: Section 2 briefly recalls the ET model considered in[11]
65 and the DG model. A full self-adjoint formulation of both models is then given in Section 3. For the sake of 66 clearness, we also extend our previous adaptive finite-element algorithm[11]to the present model in Section 67 4. In Section 5, numerical results of simulation on various diodes to compare with the results in the liter-68 ature and on MOSFET device structures to demonstrate the effectiveness of the proposed model. A short 69 concluding remark is given in Section 6.
70 2. The energy transport and density gradient models 71 As in[11], we consider the following ET model
D/¼q es ðn p þ N A N þ DÞ; ð1Þ
UNCORR
ECTE
D P
ROO
F
1 qr Jn¼ R; ð2Þ 1 qr Jp¼ R; ð3Þ r Sn¼ Jn E n xn x0 snx ; ð4Þ r Sp¼ Jp E p xp x0 spx ; ð5Þ87 where / is the electrostatic potential, n and p are the electron and hole concentrations, q is the elemen-88 tary charge, es is the permittivity constant of semiconductor, NA and NþD are the densities of ionized
89 impurities, Jn and Jp are the current densities, R is the function describing the balance of generation
90 and recombination of electrons and holes, Sn and Sp are the energy fluxes for carriers, E is the electric
91 field, snxand spx are the carrier energy relaxation times, x0 is the thermal energy, and xn and xp are
92 the carrier average energies. These physical variables are tightly coupled together with the following aux-93 iliary relationships E¼ r/; ð6Þ Jn¼ qlnnr/ þ qDnrn ¼ qnvn; ð7Þ Jp¼ qlppr/ qDprp ¼ qpvp; ð8Þ Sn¼ Jn qxnþ Jn qkBTnþ Qn; ð9Þ Sp¼ Jp þqxpþ Jp þqkBTpþ Qp; ð10Þ x0¼ 3 2kBTL; ð11Þ xn ¼ 3 2kBTnþ 1 2m njvnj 2 ; ð12Þ xp ¼ 3 2kBTpþ 1 2m pjvpj 2 ; ð13Þ Qn¼ jnrTn; ð14Þ Qp¼ jprTp; ð15Þ jn¼ 2 kB q 2 nqlnTL; ð16Þ jp¼ 2 kB q 2 pqlpTL; ð17Þ
UNCORR
ECTED
PROO
F
R¼ np n 2 i s0 nðp þ pTÞ þ s0pðn þ nTÞ ; ð18Þ126 where Qnand Qpare the heat fluxes for carries, kBis BoltzmannÕs constant, Tn, Tp, and TLare the electron,
127 hole and lattice temperatures, lnand lpare the field-dependent electron and hole mobilities, Dnand Dpare
128 the electron and hole diffusion coefficients expressed by the Einstein relation with the mobilities, m nand mp
129 are the electron and hole effective masses, vnand vpare the electron and hole velocities, jnand jpare the
130 electron and hole heat conductivities, and(18)is the Shockley–Read–Hall (SHR) generation-recombination 131 model with nibeing the intrinsic carrier concentration, s0nand s
0
p the electron and hole lifetimes, and pTand
132 nTthe electron and hole densities associated with energy levels of the traps. In the above equations, vectors
133 are denoted by bold letters.
134 The DG theory was developed by observing that the electron gas is energetically sensitive not only to its 135 density but also to the gradient of the density. It captures the nonlocality of quantum mechanics to the low-136 est-order of h2 where h is the reduced Planck constant and can be rigorously derived from quantum 137 mechanics[4,3]. Specifically, a third order derivative term of quantum correction is added to the carrier cur-138 rent density as Jn ¼ qlnnr/ þ qDnrn 2qlnbnnr Dpffiffiffin ffiffiffi n p ; ð19Þ Jp ¼ qlppr/ qDprp þ 2qlpbppr D ffiffiffipp ffiffiffi p p ; ð20Þ 0 100 200 300 400 500 600 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 X–Axis (nm) Electrostatic Potential (V) 0 100 200 300 400 500 600 0.2 0.3 0.4 0.5 1.0 2.0 3.0 4.0 5.0 X–Axis (nm) Electron Concentration (10 17 cm –3 ) ET DGET 0 100 200 300 400 500 600 –4E–4 –2E–4 0 2E–4 4E–4 6E–4 6E–4 10E–4 12E–4 14E–4 X–Axis (nm)
Electron Quantum Potential (V)
0 100 200 300 400 500 600 0 500 1000 1500 2000 2500 3000 3500 X–Axis (nm) Electron Temperature (K) (a) (b) (d) (c)
UNCORR
ECTE
D P
ROO
F
145 where the coefficients bn¼ h
2 12m nqand bp¼ h2 12m
pq are the material parameters measuring the strength of the 146 gradient effects in the gas. To alleviate the difficulty in discretization caused by this higher order term, addi-147 tional variables called the quantum potentials
/qn 2bn Dpffiffiffin ffiffiffi n p ; ð21Þ /qp 2bp D ffiffiffipp ffiffiffi p p ð22Þ 154 have been introduced in [29]and thus can be lumped with the classical drift term to obtain
Jn¼ qlnnrð/ þ /qnÞ þ qDnrn; ð23Þ
Jp¼ qlpprð/ þ /qpÞ qDprp: ð24Þ
160 We thus have a complete set of seven PDEs(1)–(5)and(21) and (22)describing both ET and DG models 161 with the seven state variables /, n, p, /qn, /qp, Tn, and Tp.
162 Note that the coefficients in(21)and(22)result in a boundary layer near the silicon/silicon-oxide interface 163 for short-channel devices. The layer is only a fraction of the length scale of the inversion layer, in which the 164 electron density typically drops from its peak value of order 1018at about 0.5–1.5 nm away from the interface 165 to zero at the interface[1,7]. Numerical treatments for this boundary layer problem are evidently subtle and 166 challenging. A more detailed description of our approach to this problem will be given in Section 5.
0 20 40 60 80 100 120 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 X–Axis (nm) Electrostatic Potential (V) 0 20 40 60 80 100 120 0.3 0.4 0.5 1.0 2.0 3.0 4.0 5.0 X–Axis (nm) Electron Concentration (10 18 cm –3 ) ET DGET 0 20 40 60 80 100 120 –2E–3 0 2E–3 4E– 3 6E–3 8E–3 X–Axis (nm)
Electron Quantum Potential (V)
0 20 40 60 80 100 120 0 500 1000 1500 2000 2500 3000 3500 4000 4500 X–Axis (nm) Electron Temperature (K) (a) (b) (d) (c)
UNCORR
ECTED
PROO
F
167 Remark 2.1. Taking first three moments of the Boltzmann transport equation (BTE) with conservation of 168 particles, momentum, and energy, the classical hydrodynamic (CHD) model can be expressed as (for 169 simplicity, we list the equations of electrons only)[6,17]:
on otþ r ðnvnÞ ¼ on ot c ; opn ot þ vnr pnþ pn rvn¼ qnE rðnkBTÞ þ opn ot c ; oxn ot þ r ðvnxnÞ ¼ qnvn E r ðvnnkBTÞ r Qnþ oxn ot c ;
172 where pn ¼ mnnvnis the momentum density. Considering the steady state and employing the collision terms
opn ot c ¼ pn spn ; oxn ot c ¼ xn x0 snx ; 175 we have[11] Jn ¼ qln kBTn q rn þ nr kBTn q /
178 and Eq.(4). Similarly the three conservation equations of the QHD model are[16,18]
0 5 10 15 20 25 30 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 X–Axis (nm) Electrostatic Potential (V) 0 5 10 15 20 25 30 0.2 0.3 0.4 0.5 1.0 2.0 3.0 4.0 5.0 X–Axis (nm) Electron Concentration (10 19 cm –3 ) ET DGET 0 5 10 15 20 25 30 –0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 X–Axis (nm)
Electron Quantum Potential (V)
0 5 10 15 20 25 30 0 500 1000 1500 2000 2500 3000 3500 X– Axis (nm) Electron Temperature (K) (a) (b) (d) (c)
UNCORR
ECTE
D P
ROO
F
on otþ r ðnvnÞ ¼ on ot c ; opn ot þ vnr pnþ pn rvnþ n 3rQ ¼ qnE rðnkBTÞ þ opn ot c ; oxn ot þ r ðvnxnÞ ¼ qnvn E r ðvnnkBTÞ r Qnþ oxn ot c :181 The quantum correction to the momentum equation is related to the quantum potential of Bohm[27]
Q¼ h 2 2m n Dpffiffiffin ffiffiffi n p ;
184 and to the energy density given by xn ¼ 3 2kBTnþ 1 2m njvnj2 h2n 24m n DlogðnÞ:
187 Following the previous deductive procedure the quantum correction current density equation is Jn¼ qln kBTn q rn þ nr kBTn q / h 2 6m nq r D ffiffiffi n p ffiffiffi n p ¼ qlnnrð/ þ /qnÞ þ qDnrn þ lnkBnrTn: –0.01 –0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 X–Axis
Electron Quantum Potential
φqn (V) 30nm 120nm 600nm Junction Junction
UNCORR
ECTED
PROO
F
190 Compared to the DG model, there is a mechanism that can cause an increase in diffusion, i.e., particle 191 diffusion is enhanced when Tnis significantly greater than TL. Once we obtain the information of Tnfrom
192 the DGET model, we will use this formulation to estimate the drain current and to sketch the I–V curves. 193 On the other hand, since we do not add the quantum correction to the energy density, the difference of the 194 temperature distribution between the ET and DGET models is not very significant.
195 Remark 2.2. A quantum energy balance model appears to be first proposed by Grubin and Kreskovsky in
196 [18]for 1D mesoscopic structures. In their model, quantum correction terms are explicitly included in the
197 carrier average energies(11)and(12). As a result, third order derivative terms of correction are associated 198 not only with the carrier densities (see(19)and(20)) but also with the carrier energies. Putting these cor-199 rection terms into our model, i.e. into(9)and(10), we will obtain a product of the correction terms in(9)
200 and(10)which obviously makes computations more formidable for 2D simulation. Instead, the correction 201 terms in our model are only explicitly added to the carrier density. The energy balance equations are implic-202 itly and thus less corrected by the quantum effects via the carrier current densities.
203 3. A self-adjoint formulation of the DGET model
204 PDEs in self-adjoint form are analytically as well as numerically appealing. In[11,12], we give a rather 205 thorough study of the self-adjoint DD and ET models in terms of mathematical analysis and numerical jus-206 tification. We now consider the self-adjoint formulation of the above DGET model and, for this purpose, 207 introduce the following new variables
UNCORR
ECTE
D P
ROO
F
u¼ exp un VT ; ð25Þ v¼ exp up VT ; ð26Þ fn ¼ ffiffiffi n p ; ð27Þ fp ¼ ffiffiffip p ; ð28Þ gn ¼ Tn=exp 5un 4VT ; ð29Þ gp ¼ Tp=exp 5up 4VT ; ð30Þ221 where VT= (kBTL)/q is the thermal voltage and unand upare the generalized quasi-Fermi potentials that
222 include the QM effects as shown below. Assuming a Maxwell–Boltzmann energy distribution of carriers, we 223 have the quantum correction expressions of the carriers
0 20 40 60 80 100 0 20 40 60 80 100 Transverse Distance (nm) Depth (nm)
UNCORR
ECTED
PROO
F
n¼ niexp / unþ /qn VT ¼ niexp /þ /qn VT u¼ f2n; ð31Þ p¼ niexp up / /qp VT ¼ niexp / /qp VT v¼ f2p; ð32Þ229 and rewrite the quantum potentials as /qn¼ VTln f2n uni /; ð33Þ /qp¼ VTln f2p vni ! /: ð34Þ
234 For Eq.(1) we have
D/¼ F ð/; u; v; fn;fpÞ; ð35Þ 237 where Fð/; u; v; fn;fpÞ ¼ qni es uexp /þ /qn VT v exp / /qp VT þqðN A NþDÞ es : ð36Þ
240 Substituting(31)into the electron current equation(23), we obtain
UNCORR
ECTE
D P
ROO
F
Jn¼ qlnnrð/ þ /qnÞ þ qDnr niexp /þ /qn VT u ð37Þ ¼ qlnnrð/ þ /qnÞ þ q Dn VT niexp /þ /qn VT u rð/ þ /qnÞ þ qDn niexp /þ /qn VT ru ¼ qDnniexp /þ /qn VT ru; ð38Þ243 which defines the generalized quasi-Fermi potential as in
Jn¼ qlnnrun; ð39Þ
247 with the quantum correction in electron concentration. Boundary conditions for this potential can be easily 248 specified. Similar expressions also exist for hole.
249 For the energy fluxes, we rewrite(9)more precisely as Sn¼ 5Jn 2qkBTn jnrTnþ Jn q 1 2m njvnj 2 : ð40Þ
252 Substituting(16),(29), and(39)into this equation, we have Sn¼ 5Jn 2qkBgnexp 5un 4VT jn exp 5un 4VT rgnþ 5 4VT gnexp 5un 4VT run þJn q 1 2m njvnj 2 ¼ jnexp 5un 4VT rgnþ Jn q 1 2m njvnj2 : ð41Þ
UNCORR
ECTED
PROO
F
255 Hence, we obtain the following self-adjoint form r jnexp 5un 4VT rgn ¼ RnðgnÞ; ð42Þ 258 where RnðgnÞ ¼ n xn x0 snx Jn E 1 qr 1 2m n jJnj 2 q2n2Jn ! : ð43Þ
261 We also have a similar equation for hole.
262 Our new model for both DG and ET equations with the seven state variables /, u, v, fn, fp, gn, and gp
263 and their associated boundary conditions (BCs) is re-organized as follows:
D/¼ F ð/; u; v; fn;fpÞ; ð44Þ 1 qr Jn ¼ Rð/; u; v; fn;fpÞ; ð45Þ 1 qr Jp ¼ Rð/; u; v; fn;fpÞ; ð46Þ Dfn¼ Znð/; u; v; fn;fpÞ; ð47Þ 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 Transverse Distance (nm) Depth (nm)
Hole Concentration (log., cm
–3
)
UNCORR
ECTE
D P
ROO
F
Dfp¼ Zpð/; u; v; fn;fpÞ; ð48Þ r Gn¼ RnðgnÞ; ð49Þ r Gp¼ RpðgpÞ; ð50Þ 285 where Fð/; u; v; fn;fpÞ ¼ qni es uexp /þ /qn VT v exp / /qp VT þqðN A N þ DÞ es ; ð51Þ Jn¼ þqDnniexp /þ /qn VT ru; ð52Þ Jp¼ qDpniexp / /qp VT rv; ð53Þ Rð/; u; v; fn;fpÞ ¼ n2 i uvexp /qn/qp VT 1 h i s0 n nivexp //qp VT þ pT h i þ s0 p niuexp /þ/qn VT þ nT h i ; ð54Þ Znð/; u; v; fn;fpÞ ¼ fn 2bn VTlnðf2nÞ VTlnðuniÞ / ; ð55Þ 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 3500 Transverse Distance (nm) Depth (nm) Electron Temperature (K)UNCORR
ECTED
PROO
F
Zpð/; u; v; fn;fpÞ ¼ fp 2bp VTlnðf2pÞ þ VTlnðvniÞ / h i ; ð56Þ /qn¼ VTlnðf2nÞ VTlnðuniÞ /; ð57Þ /qp¼ VTlnðf2pÞ þ VTlnðvniÞ /; ð58Þ Gn¼ jnexp 5un 4VT rgn; ð59Þ Gp¼ jpexp 5up 4VT rgp; ð60Þ RnðgnÞ ¼ n xn x0 snx Jn E 1 qr 1 2m n jJnj 2 q2n2Jn ! ; ð61Þ RpðgpÞ ¼ p xp x0 spx Jp E þ 1 qr 1 2m p jJpj 2 q2p2Jp ! : ð62Þ315 The boundary conditions are changed accordingly to
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 Transverse Distance (nm) Depth (nm) Hole Temperature (K)
UNCORR
ECTE
D P
ROO
F
/¼ VOþ Vb; u¼ exp VO VT ; v¼ exp VO VT ; f2n¼1 2 ðN þ D NAÞ þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðNþ D NAÞ 2 þ 4n2 i q ; fp¼ ni=fn; gn¼ 300 exp 5VO 4VT ; gp¼ 300 exp 5VO 4VT on oXD 318 and o/ on¼ ou on¼ ov on¼ ofn on ¼ ofp on ¼ ogn on ¼ ogp on ¼ 0 onoXN;321 where VOdenotes the applied voltage and Vbrepresents the built-in potential. Here, X R2 denotes the
322 bounded domain of the silicon. The boundaryoX = oXD[ oXNis piecewise smooth consisting of Dirichlet
323 oXDand Neumann oXNparts. The Dirichlet part corresponds to the ohmic contacts on the device. Note
UNCORR
ECTED
PROO
F
324 that the above Neumann BCs for fnand fpdo not hold on the entireoXNexcluding the oxide interface at
325 which a zero Dirichlet BC is imposed. As mentioned in[8], the quantum potentials would have to be infinite 326 at the interface to force the carrier densities to exactly zero there. Thus, a suitable constraint on the values 327 of the quantum potentials at the interface is also not available. A small but non-zero value of the carrier 328 densities is instead used in that paper. Our implementation of such non-exact zero Dirichlet BC at the inter-329 face will be specified in Section 5.
330 It should be noted that effective approximation of the gradient of current densities in formulas(61) and
331 (62)is in general very difficult to acquire. Simplified models for these formulas based on physical
consid-332 eration are possible. For example, by assuming that the drift energy is only a small part of the total kinetic 333 energy[10],(61) and (62)can be reduced to
RnðgnÞ ¼ n xn x0 snx Jn E; RpðgpÞ ¼ p xp x0 spx Jp E;
336 which will be used in our numerical simulations.
337 Remark 3.1. As observed in[5], the SRH generation-recombination model(18)should be modified for the 338 DG model since this standard expression will produce spurious generation and recombination near the 339 oxide barrier. We thus consider here the modified SRH (MSRH) proposed in[5]and extend it into the self-340 adjoint context as follows:
UNCORR
ECTE
D P
ROO
F
np neqpeq s0 n pþpffiffiffiffiffiffiffiffiffiffiffiffineqpeqexp etei kBT þ s0 p nþpffiffiffiffiffiffiffiffiffiffiffiffineqpeqexp etei kBT ¼ n2 iuvexp /qn/qp VT neqpeq s0 n nivexp //qp VT þpffiffiffiffiffiffiffiffiffiffiffiffineqpeqexp etei kBTL h i þ s0 p niuexp /þ/qn VT þpffiffiffiffiffiffiffiffiffiffiffiffineqpeqexp etei kBTL h i :343 Note that the term n2
i in(18)is replaced by neqpeqin this MSRH model where etand eiare the trapped and
344 intrinsic energies. The quantities neqand peqare the spatially dependent equilibrium densities obtained from
345 a separate numerical solution of the same DG problem, but with all voltages and R set to zero. Following 346 that paper, we choose s0
n¼ s 0 p¼ 10
8 s with e
t= ei in our simulation. A comparison of numerical results
347 based on both SRH and MSRH models will be given in Section 5. 348 Remark 3.2. For simplicity, we use fixed mobilities of ln= 1500 cm
2
/V s and lp= 500 cm 2
/V s which are 349 roughly equal to the intrinsic values at room temperature for silicon as considered in[8]. In our numerical 350 experiences in[11,12], the field-dependent mobility model of the Caughey–Thomas expression still can be 351 used in the DGET simulation.
352 Remark 3.3. The above self-adjoint formulation is based on Maxwell–Boltzmann statistics. However, it is 353 unclear to us whether the self-adjointness can also be derived for the case of Fermi–Dirac statistics which is 354 more exact but more complicated to implement. Evidently, this issue deserves further investigation in the 355 future. 90 91 92 93 94 95 96 97 98 99 100 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Y–Axis (nm) Electron Concentration (10 18 cm –3 ) ET DG DGET
UNCORR
ECTED
PROO
F
356 4. An adaptive finite element algorithm for the DGET model
357 The main ingredients of the algorithm solving the DGET model are adaptive finite element approxima-358 tion of the model, node-by-node and monotone iterative solution of the resulting nonlinear algebraic sys-359 tems, and GummelÕs iteration consecutively on the PDEs as described in[11]for the ET model. For the sake 360 of clearness, we briefly illustrate the algorithm and refer to[11,12]for more details on the adaptive finite 361 element formulation, monotone convergence analysis, and practical implementation issues.
362 Here, we use the notation l as GummelÕs (outer) iteration index and m as the monotone (inner) iteration 363 index.
364 Step 1. Initial mesh: create a coarse and structured mesh for which the number of nodes can be chosen as 365 small as possible.
366 Step 2. Preprocessing: see[11].
367 Step 3. Gummel and Monotone iterations on(44)–(48). 368 Step 3.0. Set l: = 0
369 Step 3.1. Solve the potential equation in(44). 370 Step 3.1.1. Set m: = 0 and set the initial guess
/ðmÞj ¼ f/j or c/j if l¼ 0; /ðlÞj otherwise; ( for allðxj; yjÞ 2 X h ;
373 where f/j and c/j are constant values that can be easily verified to be an upper and lower solution of /, 374 respectively, and Xhdenotes the set of mesh points on the closure of the domain.
375 Step 3.1.2 If l = 0, set u(l)and v(l)by the charge neutrality condition.
UNCORR
ECTE
D P
ROO
F
376 Step 3.1.3. Compute /ðmþ1Þj by solving the discrete potential system of(44)
nj/ðmþ1Þj þ cjð/Þ/ ðmþ1Þ j ¼ P k2V ðjÞnk/ ðmÞ k F ð/ m ð Þ j ; u ðlÞ j ; v ðlÞ j ;f ðlÞ n ;f ðlÞ p Þ þ cjð/Þ/ ðmÞ j ; 8ðxj; yjÞ 2 X h; /ðmþ1Þj ¼ VOþ Vb; 8ðxj; yjÞ 2 oX h D; o/ðmþ1Þj on ¼ 0; 8ðxj; yjÞ 2 oX h N; 8 > > > < > > > : ð63Þ 380 where cjð/Þ ¼ max oFð/jÞ o/ ; ^/j6/j6 ~/j ; ð64Þ
384 nkare the matrix elements of the discretization, and Xh,oXhD, andoX h
Nrepresent the sets of mesh points in
385 the interior, Dirichlet part, and Neumann part of the domain, respectively.
Step 3.1.4. Set /ðmÞj :¼ /ðmþ1Þj 8j 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 8j. 389
Step 3.2. Solve the electron continuity equation(45). Step 3.2.1. Set m: = 0 and set the initial guess
uðmÞj ¼ euj orbuj if l¼ 0; uðlÞj otherwise; ( for allðxj; yjÞ 2 X h ;
UNCORR
ECTED
PROO
F
394 whereeujandbujare constant values for allðxj; yjÞ 2 X h
that can be easily verified to be an upper and lower 395 solution of u, respectively.
Step 3.2.2. Compute uðmþ1Þj by solving the discrete electron system of(45).
Step 3.2.3. Set uðmÞj :¼ u ðmþ1Þ
j 8j 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 8j.
400
401 Step 3.3. Solve the hole continuity equation(46)similarly to that in Step 3.2. 402 Step 3.4. Solve the DG equation(47).
Step 3.4.1. Set m: = 0 and set the initial guess ½fn ðmÞ j ¼ g ½fnj or d½fnj if l¼ 0; ½fn ðlÞ j otherwise; 8 < : for allðxj; yjÞ 2 X h ;
406 where [fn]j fn(xj,yj) and g½fnjand d½fnjare constant values for allðxj; yjÞ 2 X h
that can be easily verified to 407 be an upper and lower solution of fn, respectively.
Step 3.4.2. Compute½fnðmþ1Þj by solving the discrete system of(47).
Step 3.4.3. Set½fnðmÞj :¼ ½fnðmþ1Þj 8j 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½fnðlþ1Þj :¼ ½fnðmþ1Þj 8j.
412
413 Step 3.5. Solve the DG(48)similarly to that in Step 3.4. 414 Step 3.6. Update ½/qn
ðlþ1Þ
j and½/qp ðlþ1Þ
j by the Eqs. (57) and (58).
415 Step 3.7. Set l: = l + 1 and go to Step 3.1 until the stopping criteria of the outer iteration are satisfied. 416
417 Step 4. Monotone Iteration on(49)and(50).
UNCORR
ECTE
D P
ROO
F
418 Step 4.1. Solve the energy equation(49)for gn similarly to that in Step 3.2.
419 Step 4.2. Solve the energy equation(50)for gp similarly to that in Step 3.2.
420
421 Step 5. Error estimation: See[11]. 422 Step 6. Refinement: See[11]
423 Step 7. Postprocessing: All computed solutions are then postprocessed for further analysis of physical
424 phenomena.
425
426 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 per-427 formed for the corresponding discrete system (63), for example, in which the monotone parameters(64)
428 can be easily obtained by means of lower and upper solutions. Two important factors that guarantee a glo-429 bal convergence with this kind of simple solutions as initial guesses are the diagonally dominant property of 430 the matrices due to the self-adjoint formulation and the monotonicity of the parameters by the special non-431 linearity of the formulation. The diagonally dominant property for(44)–(50)can proved in exactly the same 432 manner as that given in [11,12]. It can also be easily shown that each one of the nonlinear functionals in
433 (44)–(50)is monotone in its respective state variable. It is thus a straightforward generalization from our
434 previous theoretical analysis that all the nonlinear algebraic systems generated by this algorithm preserve 435 these two factors. We thus summarize these important results in the following theorem.
436 Theorem. For each one of the PDEs(44)–(50)with associated boundary conditions, the matrices resulting to 437 the adaptive finite element approximation are diagonally dominant. Moreover, starting with suitable lower and 438 upper solutions of the corresponding PDE, the Jacobi iteration in each of Steps 3.1–3.5, 4.1, and 4.2 generates 439 a pair of lower and upper sequences which converge monotonically to the exact solution of the nonlinear 440 algebraic system of equations of that PDE.
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 Transverse Distance (nm) Depth (nm) Hole Temperature (K)
UNCORR
ECTED
PROO
F
441 5. Numerical examples
442 To demonstrate the effectiveness and accuracy of the DGET model, several numerical studies have been 443 made for sample diode and MOSFET device structures. A benchmark device, namely, an abrupt n+–n–n+ 444 silicon diode is first used to verify our methods and formulation with the results reported in literature. 445 Numerical experiments are performed first on a 600 nm silicon diode at 300 K with n+= 5.0· 1017
cm3 446 and n = 2.0· 1015
cm3. The length of the n-region is approximately 400 nm. The steady state results 447 for this problem are illustrated by the dotted and solid curves with respective to the DGET and ET models 448 inFigs. 1(a)–(d) where the applied voltage VOis taken as 2.0 V. The dotted curve coincides with the solid
449 curve. This represents that the new model can be applied to devices with larger size, i.e., where the QM 450 effects are negligible. These results agree also very well with that previously reported in the literature
451 [2,14,17,31].
452 To verify QM effects with our model, we then reduce the scale down to two cases. Case (1) is a 120 nm 453 silicon diode with n+= 5.0· 1018cm3and n = 2.0· 1015cm3. The length of the n-region is approximately 454 80 nm. The applied voltage VOis taken as 1.2 V. Case (2) is a 30 nm silicon diode with n+= 5.0· 1019cm3
455 and n = 2.0· 1015cm3. The length of the n-region is approximately 20 nm. The applied voltage VOis
ta-456 ken as 0.8 V.Figs. 2 and 3show the significant change of the electron density predicted by the new model 457 but for the electron temperature the change is not very significant. The maximal temperatures of ET and 458 DGET models are T = 3423 K and T = 3309 K, respectively. The corresponding thermal energies are 459 Eth¼32kBT ¼ 0:442 eV and Eth= 0.428 eV. Therefore, the temperature reduced by the QM corrections
460 of the DGET model is very similar to that by the nonparabolicity effects presented in[14]. Fig. 4 shows 461 a visible tendency of the quantum potential /qntoward a large variation when the channel length is
de-462 creased. Here, we scale the figures into the same size for comparison.
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 Transverse Distance (nm) Depth (nm) Hole Temperature (K)
UNCORR
ECTE
D P
ROO
F
463 The second example of our simulation test on the model is a MOSFET device structure which has an 464 elliptical 1019cm3Gaussian doping profiles in the source and drain regions and 1016cm3in the p-sub-465 strate region as shown inFig. 5. The junction depth is 20 nm, the lateral diffusion under gate is 8 nm, the 466 channel length is 34 nm, and the gate oxide thickness is 2 nm. With VBS= 0 V, VDS= 1.0 V and VGS= 0.8
467 V,Figs. 6–13present the final adaptive mesh, electrostatic potential, electron concentration, hole concen-468 tration, electron temperature, hole temperature distribution, electron quantum potential, and hole quantum 469 potential, respectively. Across the junction,Figs. 12 and 13clearly show similar quantum potential profiles 470 as that inFig. 4for the 1D diode device. Furthermore, in the direction perpendicular to the interface, a very 471 thin boundary layer of about 6 nm appear in the inversion layer as shown inFigs. 12–17. The boundary 472 layer as shown inFig. 6is effectively captured by the a posteriori error estimation with 1-irregular refine-473 ment strategy developed in[11,12,25,26].
474 As mentioned earlier, a suitable constraint on the values of the quantum potentials at the oxide interface 475 is not available. One solution to this lack of quantum potential BCs is to solve the DGET model in the 476 oxide as well as in the adjoining silicon and poly gate. This will allow us to simulate the tunneling effects 477 across the oxide[13]. This issue is not addressed here and will be reported elsewhere in our future works. 478 We do not impose zero Dirichlet BCs for the variables fnand fpexactly at the interface but instead at the
479 grid points in silicon that are very close (about 0.13 nm) to the interface. In effect, these BCs are very similar 480 to that in[8]where a small but non-zero value of the carrier densities is set at the interface. We found that if 481 the BCs are prescribed exactly at the interface, the result of temperature will be very poor although the 482 algorithm is still convergent.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 V DS (V) IDS (mA/ µ m) ET DG DGET
UNCORR
ECTED
PROO
F
483 As noted in Remark 2.2, the quantum corrections are implicitly defined in the energy fluxes(9)and(10)
484 via the carrier densities(19)and(20), the temperature distribution of carriers will not be very accurate as 485 shown inFigs. 10 and 11. More specifically, the temperature peak for electron appears to be near the drain 486 but the peak for hole is in the middle of the channel. We found that the generation-recombination model 487 will influence the hole temperature distribution significantly. If the standard SHR model(54)is used instead 488 of the MSHR model of Remark 3.1, the hole temperature is even much worse as shown inFigs. 18 and 19. 489 To our knowledge, there are no numerical results of quantum corrected carrier temperatures available in 490 the literature to be compared with our results. Evidently, efficient and effective numerical methods for han-491 dling energy fluxes with explicit quantum corrections are needed for future investigations.
492 The electron density profile shown inFig. 14is a cross section of the 2D profile at the middle point of the 493 interface. The peak of the density is about 1.5 nm away from the interface, which agrees well with that in
494 [7], see also[1].Fig. 15is the electron current density computed by the ET model, which clearly shows that 495 the classical density is sharply peaked comparing with the smoothly peaked inFigs. 16 and 17obtained by 496 the DG and DGET models. The substantial QM effect of transconductance degradation is also evidently 497 displayed in these figures. From these figures, we observe that the carrier temperature provides a mecha-498 nism to increase carrier diffusion as noted in Remark 2.1. This is the main justification to consider DGET 499 instead of DG along.
500 Finally,Fig. 20shows the simulated I–V curves in which the drain current obtained by the DGET model 501 is about 20–60% below that predicted by the ET model for the gate biases of 0.7, 0.8, and 0.9 V. This result 502 is also in good agreement with that of[8]where a MOSFET with 30 nm gate length and 2 nm gate oxide 503 thickness is considered. Admittedly, this represents a serious decrease in the current drive capability of the 504 device. Note also that the difference of the maximal temperatures between the ET (T = 3677 K) and DGET 505 (T = 3649 K) models is not very significant. The corresponding thermal energies are 0.475 and 0.471 eV. 506 The figure also shows that ET over estimates the current whereas DG under estimates.
507 6. Conclusion
508 A self-adjoint model combining both ET and DG models is proposed here for nanoscale semiconductor 509 devices. This model is capable of describing hot carrier and quantum correction effects.
510 Moreover, due to the self-adjointness and monotonic nonlinearity, the present model enjoys many fav-511 orable mathematical properties such as global convergence with simple initial guesses, highly parallelizable, 512 and fast iterative solution. Numerical convergence is a fundamental issue constantly faced by the practi-513 tioner in device and circuit simulation. This model and monotone iterative methods may offer an alternative 514 in handling the convergence difficulties frequently associated with NewtonÕs methods.
515 Our numerical simulations on diode and MOSFET with the gate length down to 34 nm using the DGET 516 model have been performed and compared with that using the ET model. And the results are shown to be in 517 good agreement with those reported in the literature. It is shown that the I–V characteristics of this short-518 channel device is significantly corrected by the density-gradient equations with current drive reduced by up 519 to 60% comparing with that of the classical model along. Furthermore, a 2D quantum layer, which is only a 520 fraction of the length scale of inversion layer, is also effectively captured by this model with very fine mesh 521 near the interface generated by an adaptive finite element method.
522 Nevertheless, many improvements on our preliminary model can be further studied in future works. For 523 example, the self-adjoint formulation of the present paper is based on Maxwell–Boltzmann statistics. It is 524 however unclear to us whether the self-adjointness can be also derived for the case of Fermi–Dirac statistics 525 which is more exact but more complicated to implement. Moreover, efficient and effective numerical meth-526 ods for handling energy fluxes with explicit quantum corrections are also deserved for future investigations.
UNCORR
ECTE
D P
ROO
F
527 References
528 [1] H. Abebe, E. Cumberbatch, Quantum mechanical effects correction models for inversion charge and I–V characteristics of the
529 MOSFET device, in: Proceedings of the 2003 Nanotechnology Conference, vol. 2, 2003, pp. 218–221.
530 [2] N.R. Aluru, A. Raefsky, P.M. Pinsky, K.H. Law, R.J.G. Goossens, R.W. Dutton, A finite element formulation for the
531 hydrodynamic semiconductor device equations, Comput. Methods Appl. Mech. Engrg. 107 (1993) 269.
532 [3] M.G. Ancona, G.J. Iafrate, Quantum correction to the equation of state of an electron gas in a semiconductor, Phys. Rev. B 39
533 (1989) 9536.
534 [4] M.G. Ancona, H.F. Tiersten, Macroscopic physics of the silicon inversion layer, Phys. Rev. B 35 (1987) 7959.
535 [5] M.G. Ancona, Z. Yu, R.W. Dutton, P.J.V. Voorde, M. Cao, D. Vook, Density-gradient analysis of MOS tunneling, IEEE Trans.
536 Electron. Dev. 47 (2000) 2310.
537 [6] Y. Apanovich, E. Lyumkis, B. Polsky, A. Shur, P. Blakey, Steady-state and transient analysis of submicron devices using energy
538 balance and simplified hydrodynamic models, IEEE Trans. CAD 13 (1994) 702.
539 [7] A. Asenov, A.R. Brown, J.R. Watling, The use of quantum potentials for confinement in semiconductor devices, Model. Simul.
540 Microstruct. (MSMÕ2002) (2002) 490.
541 [8] B.A. Biegel, M.G. Ancona, C.S. Rafferty, Z. Yu, Efficient multi-dimensional simulation of quantum confinement effects in
542 advanced MOS devices, NAS Tech. Report NAS-04-008, 2004.
543 [9] B.A. Biegel, J.D. Plummer, Comparison of self-consistency iteration options for the Wigner function method of quantum device
544 simulation, Phys. Rev. B 54 (1996) 8070.
545 [10] D. Chang, J.G. Fossum, Simplified energy-balance model for pragmatic multi-dimensional device simulation, Solid-State
546 Electron. 41 (1997) 1795.
547 [11] R.-C. Chen, J.-L. Liu, An iterative method for adaptive finite element solutions of an energy transport model of semiconductor
548 devices, J. Comput. Phys. 189 (2003) 579.
549 [12] R.-C. Chen, J.-L. Liu, Monotone iterative methods for the adaptive finite element solution of semiconductor equations, J.
550 Comput. Appl. Math. 159 (2003) 341.
551 [13] C.-H. Choi, Z. Yu, R.W. Dutton, Modelling of MOS scaling with emphasis on gate tunnelling and source/drain resistance,
552 Superlattices Microstruct. 27 (2000) 191.
553 [14] P. Degond, A. Ju¨ngel, P. Pietra, Numerical discretization of energy-transport models for semiconductors with non-parabolic
554 band structure, SIAM Sci. Comput. 22 (2000) 986.
555 [15] D.K. Ferry, J.-R. Zhou, Form of the quantum potential for use in hydrodynamic equations for semiconductor device modeling,
556 Phys. Rev. B 48 (1993) 7944.
557 [16] C.L. Gardner, The quantum hydrodynamic model for semiconductor devices, SIAM J. Appl. Math. 54 (1994) 409.
558 [17] C.L. Gardner, J.W. Jerome, D.J. Rose, Numerical methods for the hydrodynamic device model: subsonic flow, IEEE Trans.
559 CAD 8 (1989) 501.
560 [18] H.L. Grubin, J.P. Kreskovsky, Quantum moment balance equations and resonant tunnelling structures, Solid-State Electron. 32
561 (1989) 1701.
562 [19] W. Ha¨nsch, T. Vogelsang, R. Kircher, M. Orlowski, Carrier transport near the Si/SiO2interface of a MOSFET, Solid-State 563 Electron. 32 (1989) 839.
564 [20] T. Ho¨hr, A. Schenk, A. Wettstein, W. Fichtner, On density-gradient modeling of tunneling through insulators, IEICE Trans.
565 Electron. E86-C (2003) 379.
566 [21] S. Jallepalli, J. Bude, W.K. Shih, M.R. Pinto, C.M. Maziar, A.F. Tasch Jr., Electron and hole quantization and their impact on
567 deep submicron p- and n-MOSFET characteristics, IEEE Trans. Electron Devices 44 (1997) 297.
568 [22] J.P. Kreskovsky, H.L. Grubin, VLSI Design 3 (2) (1995) 179.
569 [23] R. Lake, G. Klimeck, R.C. Bowen, D. Jovanovic, J. Appl. Phys. 81 (1997) 7845.
570 [24] Y. Li, J.-L. Liu, S.M. Sze, T.-S. Chao, A new parallel adaptive finite volume method for the numerical simulation of
571 semiconductor devices, Comput. Phys. Commun. 142 (2001) 285.
572 [25] J.-L. Liu, On weak residual error estimation, SIAM J. Sci. Comput. 17 (1996) 1249.
573 [26] J.-L. Liu, I.-J. Lin, M.-Z. Shih, R.-C. Chen, M.-C. Hsieh, Object oriented programming of adaptive finite element and finite
574 volume methods, Appl. Numer. Math. 21 (1996) 439.
575 [27] C. Philippidis, D. Bohm, R.D. Kaye, The Aharonov–Bohm effect and the quantum potential, Il Nuovo Cimento 71B (1982)
576 75.
577 [28] R. Pinnau, A. Unterreiter, The stationary current–voltage characteristics of the quantum drift-diffusion model, SIAM J. Numer.
578 Anal. 37 (1999) 211.
579 [29] C.S. Rafferty, B. Biegel, Z. Yu, M.G. Ancona, J. Bude, R.W. Dutton, Multi-dimensional quantum effect simulation using a
580 density-gradient model and script-level programming techniques, Proc. SISPAD (1998) 137–140.
UNCORR
ECTED
PROO
F
582 [31] M. Rudan, F. Odeh, J. White, Numerical solution of the hydrodynamic model for a one-dimensional device, COMPEL 6 (1987)
583 151.
584 [32] L. Shifren, R. Akis, D.K. Ferry, Correspondence between quantum and classical motion: comparing Bohmian mechanics with a
585 smoothed effective potential approach, Phys. Lett. A 274 (2000) 75.
586 [33] M.J. van Dort, P.H. Woerlee, A.J. Walker, A.H. Juffermans, H. Lifka, A simple model for quantization effects in heavily-doped
587 silicon MOSFETÕs at inversion conditions, IEEE Trans. Electron. Dev. 39 (1992) 932.
588 [34] Z. Yu, R.W. Dutton, R.A. Kiehl, Circuit/device modeling at the quantum level, IEEE Trans. Electron. Dev. 47 (2000) 1819.
