半導體元件模擬之平行計算---次微米元件模擬程式的發展與平行計算

全文

(1)

(2) Development and Parallel Computation of Sub-Micron Device Simulator

(3) 87/08/01 88/07/31 !"#$%&'( E-mail [email protected] ()*+,-./0123 4567809:6780;<-=>?6 @0A6BCDE09FG) HIJKLMN9:678EOP QRSTUVWO,-./78GX#YZ0X #[Z0X#\]OMH678^T_1234 5678Gàbc,!dTe,-Ffg h&iOjk]lmKnopOq,cKfr& i O A,es,tù,GX#Òvwx Ayz{;<-=>?6@|}f~678 y Ecf A6BCDEbctKGt6E O68PQRS U O9FG. This formulation is very suitable to the implementation of parallel computation since each row or column of the grid points in a semiconductor device can be treated independently.

(4) ¡¢STUVWO£¤¥Q¦O§¨ If©ª«¬®Ojk¯f°±²³K´ OVW./µU¶·¸¹GºNjk.8»V W¼½¾¿ÀÁfrÂ./vw ÃÄÅ·Gq Æi°ÇÈ9ÉkÊËIKL³OV W./µUfÌ²ÍÎOvwÏÂÐÍÑÒO ¼ÓG ` a q Ô V W § ¨ Õ ¦ & Ö ¡ × O À Á O 0Ø ÔfÙ Ú K É k Û Ü u ] . / y r 4,-6EÝ Þ c O 6 7 8f¤ ¤ ß Ç à á ] O F v wfq Æ â ß ã ä å i 9 É k \ æ O ç è " é ê I o N 9 c O 678 Þ E G 9F>ëyì×"éOI¿íf q`QÉkÛî7ÊËïðñòófSTUVW ôi]#æ!õOâföÇ²iO"éÊË ÷½õÂÐø\O[ZùúfûiãäN©É kÛ»9üýOrG. (Keywords: numerical simulation, Boltzmann transport equation, balance equation method, Crank-Nicolson scheme, alternating-direction-implicit method, parallel computation.) A numerical simulation program for submicron semiconductor devices has been developed based on the balance equation method. The basic equations for the carrier concentration, the carrier velocity, the carrier energy are derived from the Boltzmann transport equation. To simplify the mathematical derivation and numerical computation, it is convenient to normalize all physical quantities with respect to some proper factors so that we have equations of pure numbers. The transient behavior is descritized according to the Crank-Nicolson scheme. The nonlinear equation system is linearized by Newton iteration and the two dimensional problem is reduced to one dimension by an alternating-direction-implicit method.. The classical carrier transport problem is based on the solution of Boltzmann transport equation v v for the distribution function f (r , p ,t ) in the pov v sition r and momentum p spaces as a function of time t: ∂f  ∂f  v = −v ⋅ ∇f − e∇ ϕ ⋅ ∇ vp f +   (1) ∂t  ∂t C where v is carrier velocity, e is the magnitude of. 1.

(5) electronic charge, and ϕ is the electrostatic poten-. ∂wd v v  ∂w  = − vd ⋅ ∇ wd + e ∇ϕ ⋅ vd +  d  ∂t  ∂t c. tial arising from the space charges in semiconductor according to the Poisson equation. e ∇ 2ϕ = − (N D − n) (2). (7)  v 2 κ  m * v2d  − ∇ ⋅  nvd − ∇  wd −  3n 2  kB   The collision terms in the balance equations can be represented by the ensemble relaxation rates. (∂n / ∂t )C = −ν n( n − n0 ) ,. εs. where εs is the permittivity, N D is the dopant concentration, and n is the electron concentration. The balance equations for the carrier density n, v the momentum density np d , and the energy den-. (∂vvd. (∂wd / ∂t )C. sity nwd can be obtained by integrating both. νn = 0. ν p = ν imp + (ν p3wd3 + ν p2wd2 + ν p1wd + ν p0 )(1 − ewd ) ν w = (ν w1wd + ν w0 ) ⋅ (1 − ewd ). 2  1  − ∇ nwd − m * nvd2  (4) 3  2  ∂(nwd ) v v  ∂( nwd )  = −∇ ⋅ (nvd wd ) + en∇ ϕ ⋅ vd +   ∂t  ∂t C. (. whereν imp = 5 × 1013 N D / 1025. 0.4. ,. ν p2 = 7.590 × 1013 , ν p3 = 1.174 × 1013 , ν w0 = 2.500 × 1012 , ν w1 = 2.338 × 1011 , and wc = 0.03.. m* is the effective mass, kB is the Boltzmann constant, and κ ≅ 5k B2nT / 2m * ν p ( w) is the thermal. To simplify the mathematical derivation and numerical computation, it is convenient to normalize all physical quantities with respect to some factors so that we have equations of pure numbers. In general, the dopant and carrier concentrations are normalized to the intrinsic carrier concentration ni ,. conductivity. v Since the quantities npd and nwd on the left-hand-side of (4) and (5) are the products of two unknown quantities, these equations are slightly inconvenient. We must derive the equav tions in terms of the average carrier velocity vd and. 2 1   ∇ nwd − m * nvd2  3nm *  2 . ).  − q( w − w c )  ),1 , ew = min  exp( kBT   12 ν p0 = 6.782 × 10 ,ν p1 = 6.063 × 1013 ,. 2  v κ  m * v2d  − ∇ ⋅  nvd − ∇  wd −  (5) kB  3  2  v v where vd = pd / m * is the average drift velocity,. −. = −ν w( wd − w0) .. The Monte Carlo method is generally used to compute the ensemble relaxation rates as a function of energy. The results for silicon can be expressed as:. sides of the Boltzmann transport equation (1) for 1, v p and w over the momentum space. ∂n v  ∂n  = −∇ ⋅ (vd n ) +   (3) ∂t  ∂t C v v ∂(npd )  ∂ (npd )  v v = −∇ ⋅ (nvd pd ) + en ∇ϕ +   ∂t  ∂t C. average carrier energy wd v ∂vd e  ∂v  v v = − vd ⋅ ∇ vd + ∇ϕ +  d  ∂t m*  ∂t c. v / ∂t )C = −ν p vd ,. the potential to the thermal potential VkT = k B T / e , the energy to the thermal energy w0 = kB T , the velocity to v0 = 2kbT / m * , the distance to the Debye length x 0 = ε s kB T / e 2 ni ,. (6). the time to t 0 = x 0 / v0 . Therefore, (2), (3), (6), and (7) become ∇ 2ϕ + ( N d − n) = 0 ∂n v = −∇ ⋅ (vd n) − ν n (n − n 0 ) ∂t. 2. (8) (9).

(6) v ∂vd v v 1 = −vd ⋅ ∇ vd − ∇ nwd − nvd2 ∂t 3n 1 v + ∇ ϕ − ν pvd 2 v v ∂wd 2 = −vd ⋅ ∇wd − ∇ nvd wd − v2d ∂t 3n v 5 2 + ∇ wd − v2d + ∇ϕ ⋅ vd −ν w (wd − w0) 6ν p. (. [ (. (. Φv y (ηi , j −1,ηi −1, j ,ηi, j ,ηi +1, j ,ηi , j +1 ) =. ). (10). )]. ). (11). 2(ηi, j −η i0, j ) − Gi, j + Gi0, j (21) = ∆t where  ∂Gn,i , j ∂Gn,i, j ∂Gn, i, j ∂Gn,i , j ∂Gn,i, j   ∂vxk ,l ∂vyk ,l ∂wk ,l ∂ϕ k ,l   ∂nk ,l  ∂Gvx,i , j ∂Gvx ,i , j ∂Gvx ,i, j ∂Gvx, i, j ∂Gvx ,i , j   ∂vxk ,l ∂vyk ,l ∂wk ,l ∂ϕ k ,l   ∂nk ,l ∂Gi, j ∂Gvy, i, j ∂Gvy ,i , j ∂Gvy ,i, j ∂Gvy, i, j ∂Gvy ,i , j  ≡ ∂vxk ,l ∂vyk ,l ∂wk ,l ∂ϕ k ,l  ∂ηk ,l  ∂nk ,l  ∂Gw,i, j ∂Gw,i , j ∂Gw,i , j ∂Gw,i, j ∂Gw,i , j   ∂n ∂vxk ,l ∂vyk ,l ∂wi, j ∂ϕ k ,l  k ,l    ∂Gϕ , i, j ∂Gϕ ,i , j ∂Gϕ ,i, j ∂Gϕ ,i, j ∂Gϕ ,i , j   ∂n ∂vxk ,l ∂vyk ,l ∂wk ,l ∂ϕ k ,l   k ,l A simple alternating-direction-implicit (ADI) method is employed and the two dimensional problem can be reduced to one dimension under this condition.. [. )]. ). In the two dimensional problem, the velocity v vd should resolve along x and y directions as vx v and vy , while the function Gv as Gvx and Gvy . The equation system can be written in the following forms: Φ n (ηi , j −1, ηi −1, j ,η i , j ,η i +1, j , ηi , j +1 ) =. Φv x. ∆t. (η. i, j −1. − (Gn,i , j + Gn0,i, j ) = 0. ∆t. (16). − (Gvx ,i, j + Gvx0 ,i, j ) = 0. ]. (a) Along y-direction (j inner loop, i outer loop). ,ηi −1, j ,η i, j ,ηi +1, j ,ηi, j +1) =. 2(vx,i , j − vx0,i , j ). (20). The nonlinear equation system is solved by the Newton iteration method.  ∂G ∂Gi , j ∂Gi, j 2 δη i , j −1 + δη i−1, j +  i , j − δη i, j ∂ηi −1, j ∂η i, j −1  ∂ηi, j ∆t  ∂Gi, j ∂Gi, j δη i +1, j + δη i, j +1 + ∂ηi +1, j ∂ηi , j +1. Gϕ = ∇ 2ϕ + ( N d − n). 2(ni , j − ni0, j ). (19). vx ,i, j , vy,i , j , wi, j and ϕ i , j at the grid point (i,j).. v v 2 G w = − vd ⋅ ∇wd − ∇ nvd w d − vd2 3n 5 2 v + ∇ wd − vd2 + ∇ϕ ⋅ vd − ν w (wd − w0) 6ν p. (. 2(ni, j − ni0, j ) − (Gw,i , j + G0w, i, j ) = 0 t ∆ Φϕ (ηi , j −1,ηi−1, j ,ηi , j ,ηi +1, j ,ηi, j +1) =. where ηi , j represents the unknown variables ni, j ,. ). [ (. (18). Gϕ ,i, j + Gϕ0,i , j = 0. The Crank-Nicolson method will be employed to solve the transient problem. ni , j − ni0, j Gn,i, j + Gn0,i , j (12) = ∆t 2 v v v v vi, j − vi0, j Gv ,i, j + Gv0,i, j (13) = ∆t 2 wi , j − w0i, j Gw,i, j + Gw0,i, j (14) = ∆t 2 Gϕ ,i, j + Gϕ0,i, j =0 (15) 2 where v G n = −∇ ⋅ (vd n) − ν n ( n − n 0 ) v 2 v v Gv = −vd ⋅ ∇vd − ∇ nwd − nvd2 3n 1 v + ∇ϕ −ν p (wd )vd 2. (. 2(vy ,i, j − v0y ,i, j ) − (Gvy ,i, j + Gvy0 ,i, j ) = 0 ∆t ( Φ w η i, j −1,ηi −1, j ,ηi, j ,η i+1, j ,ηi , j +1 ) =. (17). 3.

(7) ∂Gi, j ∂η i, j −1 = −.  ∂Gi, j 2  ∂Gi , j δη i , j+1 − δη i , j +  ∂η  η t ∆ ∂ i, j +1  i,j . of Numerical Algorithms in Semiconductor Device Simulation”, Solid State Electron, Vol. 30, No. 8, pp.813-820, 1987. [3] J. M. Ortega, “Introduction to Parallel and Vector Solution of Linear System”, Second Edition, 1988. [4] T. W. Tang, “History and Recent Develop-. δη i , j −1 + . 2(ηi, j −η i0, j ). ∆t. [. − Gi, j + Gi0, j. ]. ∂Gi , j ∂Gi, j δη i0−1, j − δη i0+1,j ∂ ηi +1, j ∂η i−1, j. (22). (b) Along x-direction (i inner loop, j outer loop)  ∂G ∂Gi, j ∂Gi , j 2  δη i −1, j +  i, j − δη i , j + δη i +1, j  ∂η  ∂η i−1, j ∂ηi +1, j  i , j ∆t . = −. 2(ηi, j −η i0, j ). ∆t. [. − Gi, j + Gi0, j. ]. ∂Gi , j ∂Gi, j δη i0, j −1 − δη i,0j +1 ∂ηi, j +1 ∂η i, j −1. (23). This formulation is suitable to the implementation of parallel computation since each row or column of the grid points in a semiconductor device can be treated independently.

(8) - !"#$%&'()*+,-./ 012 3456789: );<= >[email protected])CDEFGH & 'IJKLMN OP)QR ST U V4 GWXYZ[\ &'I]^ )_àbn + nn + U = cd 6 UV4e < f e g h i j k l L f m L n " 0.5v o p h q r s t u)J e v w j k n " 0.3v o p x h q r s t u 4e y j k x z {op|}~tu ~ 4. ment of Hydrodynamic Transport Models for Semiconductor Device Simulations”, ISSDT '95, South Africa, 1995. e< l L h q r s t u.

(9) [1] K. Tomizawa, “Numerical Simulation of Submicron Semiconductor Devices”, Artech House, 1993. [2] A. Yoshi, M. Tomizawa and K. Yokoyama. “Investigation. eg m L h q r s t u. 4.

(10) ev x h q r s t u ey x z { o p r s t u. 5.

(11)