Three-dimensional quantum simulation of multigate nanowire field effect transistors

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 79 (2008) 1060–1070

Three-dimensional quantum simulation of multigate nanowire field effect transistors Mincheol Shin School of Engineering, Information and Communications University, 119 Munji-ro, Yuseong, Daejeon 305-732, Republic of Korea Received 15 October 2007; accepted 15 October 2007 Available online 23 October 2007

Abstract Detailed numerical methods for the three-dimensional quantum simulation of the multigate nanowire field effect transistors in the ballistic transport regime are presented in this work. The device has been modeled based on the effective mass theory and the non-equilibrium Green’s function formalism, and its simulation consists of solutions of the three-dimensional Poisson’s equation, two-dimensional Schr¨odinger equations on the cross-sectional planes, and one-dimensional transport equation. Details on numerical techniques for each of the simulation steps are described, with a special attention to the solution of the most CPU demanding two-dimensional Schr¨odinger equation. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Quantum simulation; Silicon nanowire; Field effect transistors; Non-equilibrium Green’s function

1. Introduction The multigate nanowire field effect transistors (NWFETs), which have multiple gates around the silicon nanowire channel, are promising candidates for the next generation transistors and have drawn considerable attention recently [7,4]. In addition to the effective suppression of short channel effects due to the improved gate strength, the multigate NWFETs show excellent current drive and have the merit that they are compatible with conventional CMOS processes [8]. To simulate the devices, accurate modeling and calculations based on quantum mechanics are necessary to assess their performance limits, since cross-sections of the multigate NWFETs are expected to be a few nanometers wide in their ultimate scaling. We have developed an in-house quantum mechanical simulator for the multigate NWFETs based on the effective mass theory and non-equilibrium Green’s function (NEGF) method in the ballistic regime. In this paper, details on the numerical methods employed in our simulator are presented. In Section 2, the approach that has been taken to simulate the device, which is similar to pervious works [9,3], is given. The self-consistent scheme and how the threedimensional (3D) problem is split into two-dimensional (2D) Schr¨odinger and one-dimensional (1D) transport problems are described in the Section. In Section 3, we present how we have implemented numerical solutions involved. The conventional k-space approach together with the product-space approach that we have developed are described with the mass discontinuity across the silicon/oxide boundaries taken into account. The solution of the sparse matrices in a E-mail address: [email protected]. 0378-4754/$32.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.10.007

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

1061

Fig. 1. Simulated device: multigate nanowire field effect transistors.

parallel computing environment and the numerical scheme to efficiently evaluate the charge densities and current in the NEGF formalism are also important topics in the Section. In Section 4, we briefly demonstrate results of simulations on the multigate NWFETs. Finally, we conclude in Section 5. 2. Approach The simulated device is a three-dimensional nanowire FET with source/drain, channel, and multiple gates as shown in Fig. 1. The source/drain are modeled as semi-infinite leads which are heavily n-doped silicon, with typical doping concentration of 1020 cm−3 . The silicon channel may be either n-type or p-type, and gates are assumed to be metallic with the mid-gap work-function. The effective mass Hamiltonian of the device is written as H 3D Ψ (x, y, z) = EΨ (x, y, z), where H

3D

h ¯ 2 ∂2 h ¯2 ∂ =− ∗ 2 − 2mx ∂x 2 ∂y

(1)



1 ∂ m∗y ∂y

 −

h ¯2 ∂ 2 ∂z



1 ∂ m∗z ∂z

 + V (x, y, z),

(2)

where m∗x , m∗y , and m∗z are effective masses in the x, y, and z directions, respectively, and V (x, y, z) is the conduction band-edge profile, given by V (x, y, z) = Ec0 (y, z) − φ(x, y, z),

(3)

where Ec0 (y, z) is half of the band gap of silicon or oxide, depending on whether the point (y, z) belongs to the silicon region or oxide region, respectively, and φ(x, y, z) is the vacuum level potential. Since the electron transport occurs in the x direction, m∗x is simply the silicon effective mass in the transport direction, but m∗y = m∗y (y, z) and m∗z = m∗z (y, z), which can be either silicon or oxide effective masses, depending on whether the point (y, z) belongs to the silicon region or oxide region, respectively. The wavefunction penetration into the oxide region is taken into account as a consequence. In the cross-sectional plane located at x, the wavefunction Ψ (x, y, z) of the 3D Hamiltonian can be expanded as  Ψ (x, y, z) = ϕm (x)ψm (y, z; x), (4) m

where ψm (y, z; x) is the mth mode eigenfunction of the two-dimensional Schr¨odinger equation: H 2D ψm (y, z; x) = Em (x)ψm (y, z; x), where H

2D

h ¯2 ∂ =− 2 ∂y



1 ∂ ∗ my (y, z) ∂y



h ¯2 ∂ − 2 ∂z

(5) 

1 ∂ ∗ mz (y, z) ∂z

 + V (y, z; x),

(6)

subject to the boundary condition that the wave functions vanish at the boundaries of the 2D cross-sectional plane. In the so called uncoupled mode space approach [9], the coupling between different modes (or subbands) is ignored, and

1062

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

ϕm (x) satisfies:   h ¯ 2 d2 − ∗ 2 + Em (x) ϕm (x) = Eϕm (x). 2mx dx

(7)

Therefore the original 3D problem is split into problems of solving the 2D Schr¨odinger equations in the cross-sectional planes and the 1D Schr¨odinger equation in the transport direction. Note that the uncoupled mode-space approach is well suited for the device simulation considered here, because the silicon-on-insulator structure ensures the confinement of the cross-sectional wave functions in the silicon channel and consequently the shape of the 2D wave functions do not change much in the lateral direction. The Schr¨odinger equation in Eq. (7) with open boundary conditions, which is a 1D transport problem, has been solved using the NEGF method [5] as follows. We first write the 1D Green’s function Gm for subband m:

−1 Gm = E − Hm1D − ΣS,m − ΣD,m , (8) where Hm1D = −

h ¯ 2 d2 + Em (x), 2m∗x dx2

(9)

and ΣS,m and ΣD,m are source (S) and drain (D) self-energies of subband m, respectively, which will be defined later. The 1D charge density n1D m (x) in the mth subband is then obtained via  1 (10) dE(fS Gm ΓS,m G†m + fD Gm ΓD,m G†m ), n1D m (x) = 2π x where x is the 1D lattice spacing, ΓS,m and ΓD,m are defined by

† ΓS,m = i ΣS,m − ΣS,m ,

† ΓD,m = i ΣD,m − ΣD,m ,

(11) (12)

and fS(D) are the Fermi distribution functions at the source and drain given by fS(D) (E) =

1 1+e

S(D)

(E−EF

)/kB T

,

(13)

S(D)

where EF are fermi energies at the source and drain, respectively. Once the 1D charge density for each subband is obtained, the 3D quantum charge density is calculated via  2 n3D (x, y, z) = n1D m (x)|ψm (y, z; x)| ,

(14)

m

which is used in the Poisson’s equation: q ∇ 2 φ(x, y, z) = − (ND (x, y, z) − n3D (x, y, z)) 

(15)

to be solved for the potential φ(x, y, z), where ND (x, y, z) is the doping profile. Equations are solved iteratively until the self-consistent potential and charge distribution is obtained. See Fig. (2) for the simulation flow chart. If the self-consistency is reached, the current is calculated by using the Landauer–B¨uttiker formula:  2q  Id = dETm (E)(fS (E) − fD (E)), (16) h m where the transmission probability Tm (E) for subband m is given by Tm (E) = Tr(ΓS,m Gm ΓD,m G†m ).

(17)

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

1063

Fig. 2. Simulation flowchart.

3. Numerical methods In this Section, we describe how we have implemented numerical solutions of the 3D Poisson’s equation, the 2D Schr¨odinger equation, and the 1D NEGF equations in detail. We first describe the widely used k-space method for solution of the 2D Schr¨odinger equation in the cross-sectional planes, including the mass discontinuity across the silicon–oxide interface. As it turns out, this part of the numerical procedure is most demanding in terms of the computation time. We therefore developed an efficient method named as the ‘product space’ method to solve the 2D Schr¨odinger equations. Details on the product space solution will be given in this Section, followed by description on numerical techniques that have been employed to solve the 3D Poisson equation and 1D NEGF equations. 3.1. Solution of 2D Schr¨odinger equations in the k-space Let us first write ψ(y, z) = ψm (y, z; x) in Eq. (5) as  ψ(y, z) = AK |K,

(18)

K

where {|K} is a basis set and AK ’s are expansion coefficients. Inserting Eq. (18) in Eq. (5) and multiplying L| to the both sides of the equation, we obtain an eigenvalue problem:  2D HLK AK = Em AL , (19) K 2D = L|H 2D |K. In the widely used k-space solution [1], we use: where HLK   2 2 sin(ki y) sin(kj z), |K = Ly Lz

(20)

1064

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

Fig. 3. Cross sectional plane of a NWFET.

where Ly and Lz are side lengths of the cross-section in the y and z directions, respectively, and ki =

πi Ly

(i = 1,...,Ny )

(21)

kj =

πj Lz

(j = 1,...,Nz ),

(22)

where Ny and Nz are number of meshes in the y and z directions, respectively. Note that the index K in Eq. (20) is mapped to indices i and j through K = Ny (i − 1) + j. To be specific, let us consider a rectangular cross-section where silicon/oxide interfaces are located at yα and yβ and at zα and zβ in the y and z directions, respectively (see Fig. 3). For the Hamiltonian of Eq. (6) with discontinuous 2D is given by effective masses across the silicon/oxide interfaces, HLK (y)

(z)

2D 0 = HLK + hLK + hLK , HLK

where 0 HLK

4 = Ly Lz

 0

Ly

 0

(23) 

Lz

dy dz sin(kp y) sin(kq z)



h ¯ 2 kj2 h ¯ 2 ki2 + + V (y, z) 2m∗y (y, z) 2m∗z (y, z)

× sin(ki y) sin(kj z), (24)

and (y) hLK

2¯h2 ki =− (sin(kp yα ) cos(ki yα ) − sin(kp yβ ) cos(ki yβ )) Ly Lz

(z) hLK

2¯h2 kj =− (sin(kq zα ) cos(kj zα ) − sin(kq zβ ) cos(kj zβ )) Ly Lz

 

1 1 − ∗ ∗ msi,y mox 1 1 − ∗ ∗ msi,z mox



 ×



 ×

zβ

zα

dz sin(kj z) sin(kq z), (25)

yβ

yα

dy sin(ki y) sin(kp y),

where the index L is mapped to indices p and q in a similar fashion to the index K, m∗si,y and m∗si,z are the effective masses of silicon in the y and z directions, respectively, and m∗ox is the effective mass of oxide. Note that the mass (y) (z) 0 in Eq. (24) can be efficiently discontinuity is accounted for in the hLK and hLK terms of Eq. (25). Numerically, HLK evaluated using the FFT routines [1]. The eigenvalue problem of Eq. (19) can be solved using the standard eigenvalue routines such as dsyevx() of LAPACK, and the size of the problem is Ny Nz × Ny Nz . For the nanowire cross-section of area 5 nm × 5 nm, the size of the k-space meshes should be at least 32 × 32 ∼ 1000, and only the first 10 eigenvalues or so are necessary (because only they contribute to the transport). Therefore, the numerical task becomes that of finding the first few eigenvalues of a matrix of size ∼ 1000 × 1000, and since the 2D Schr¨odinger equations should be solved in every cross-sectional plane, this part of the numerical procedure is quite demanding computationally.

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

1065

3.2. Solution of 2D Schr¨odinger equations in the product space

by

First let us consider the effective masses in the rectangular cross-section of Fig. 3. The effective mass m∗y is given  m∗y (y, z)

=

where



m∗y (y) =

m∗y (y)

if zα ≤ z ≤ zβ

m∗ox

if z < zα or z > zβ

m∗si,y

if yα ≤ y ≤ yβ

m∗ox

if y < yα or y > yβ

,

.

(26)

(27)

However, in a good approximation, we may write: m∗y (y, z) = m∗y (y)

for 0 ≤ z ≤ Lz ,

(28)

for, when we consider the solution of the Schr¨odinger equation in the y direction only, wave functions are expected to penetrate into the oxide region to the left and right of the silicon region (i.e. the oxide regions of zα < z < zβ ) but the amplitude of the wave functions in the upper and lower oxide regions (z < zα and z > zβ , respectively) is expected to be small, since the oxide band gap is much higher than that of silicon. Therefore, the effective masses in the upper and lower oxide region can be regarded as not playing significant role so we may set them as above. The same approximation can be applied to the effective mass in the z direction, m∗z : m∗z (y, z) = m∗z (z) where m∗z (z)

 =

for 0 ≤ y ≤ Ly ,

m∗si,z

if zα ≤ z ≤ zβ

m∗ox

if z < zα or z > zβ

(29)

.

(30)

Now the effective masses in the left and right oxide regions are immaterial, by the same reasoning as given above. Therefore we may write, in a good approximation:     h ¯2 ∂ 1 ∂ 1 ∂ h ¯2 ∂ 2D H =− − + V (y, z). (31) 2 ∂y m∗y (y) ∂y 2 ∂z m∗z (z) ∂z In the product-space solution that we have developed, we use the basis set: |K = χi (y)ζj (z),

(32)

where χi (y) is the ith eigenfunction of the following 1D Schr¨odinger equation in the y direction:     h ¯2 d 1 d − + V¯ (y) χi (y) = i χi (y), 2 dy m∗y dy where V¯ (y) =

1 zβ − z α



zβ

dzV (y, z),

V¯ (z) =

1 yβ − y α

(34)

zα

and ζj (z) is the jth eigenfunction of the following 1D Schr¨odinger equation in the z direction     h ¯2 d 1 d − + V¯ (z) ζj (z) = j ζj (z), 2 dz m∗z dz where

(33)



(35)

yβ

dy V (y, z). yα

(36)

1066

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

If we put Eq. (32) in Eq. (31) and make use of Eqs. (33)–(36), we obtain: H 2D |K

= =

   h ¯2 d 1 dχi (y) 1 dζj (z) (y) − χ + V (y, z)χi (y)ζj (z) i m∗y (y) dy 2 dz m∗z (z) dz (i + j + V (y, z) − V¯ (y) − V¯ (z))|K.

h ¯2 d −ζj (z) 2 dy



(37)

By multiplying L| to Eq. (37), we obtain: 2D HLK = L δLK + L|(V (y, z) − V¯ (y) − V¯ (z))|K,

(38)

where L = i + j .

(39)

In our product-space solution, one therefore needs to solve the two 1D Schr¨odinger equations in the y and z directions, and then solve the above eigenvalue problem. Numerically, the 1D Schr¨odinger equations in Eqs. (33) and (35) can be solved quite easily, and the one-dimensional version of the k-space solution method that has been described in the previous subsection was employed in our simulation. Note that the potentials in the 1D equations are averaged ones in the silicon part of the plane. See Eqs. (34) and (36). For example, in the y-direction solution, V¯ (y) is obtained by averaging over the silicon region of zα < z < zβ . Since the conduction band profile in the oxide region is about 4 eV higher, in the case of silicon oxide, than that of the silicon region, the averaged potential will be obscured by the high conduction band profile of oxide if the oxide region is included in the averaging. In solving the 1D Schr¨odinger equation in one direction, however, the entire silicon and oxide regions are included, together with the mass discontinuity at the silicon/oxide interfaces. An advantage of the product-space approach is that only the first few eigenvalues, which correspond to the subband modes that contribute to the transport, are sufficient to be included in the eigenvalue problem. In other words, if M is the number of subbands that participate in the transport, the size of the matrix [H 2D ] in Eq. (38) is reduced to M by M. As mentioned previously, M is about 10 for the nanowire transistor of cross-sectional area of 5 nm × 5 nm. Therefore, compared to the k-space solution, where the first 10 eigenvalues of a 1000 by 1000 matrix were sought, one now needs to find the first 10 eigenvalues of a 10 by 10 matrix. As the area of the cross-sectional plane increase, M increases in proportion to the area. We have also found that the matrix [H 2D ] is effectively a banded matrix of band size NB which ranges from 1 to 10 for the cross sectional areas of 5 nm × 5 nm to 20 nm × 20 nm, which also contributes to the efficiency of the numerical calculation. In Fig. 4, we have shown typical cross-sectional wave functions obtained by using the product-space approach.

Fig. 4. Typical wave functions obtained by using the product-space approach: (a) the first mode and (b) the sixth mode wave functions of a top-gated NWFET.

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

1067

Fig. 5. A typical convergence behavior of the Poisson-NEGF self-consistent iterations.

3.3. Solution of 3D Poisson’s equation In our simulation, we solved the following Poisson’s equation: q k k−1 ∇ 2 φk = − (ND − nk3D e(φ −φ )/kB T ), 

(40)

instead of the one in Eq. (15), in order to achieve faster convergence of the Poisson-NEGF self-consistent calculations. In the above Gummel-type iteration [6], φk and nk3D are the kth step solutions for the potential and charge density, respectively. The boundary conditions for the Poisson’s equation are such that for the gate contact region, φ(x, y, z) = Vg , but for other boundaries including source/drain contacts, the free boundary condition was used: ∂φ = 0, ∂n¯

(41)

where n¯ denotes the direction normal to the boundary. The free boundary condition was used at the source/drain contacts to achieve the thermal equilibrium conditions in the source/drain [5]. The usual interface condition of   ∂φ  ∂φ  ox  = si  (42) ∂n¯ ox ∂n¯ si was imposed at the silicon/oxide interfaces. Eq. (40) was numerically solved by the standard Newton–Raphson method after discretized by the seven-point finite difference method on a rectangular hexahedral mesh. The resultant sparse matrix was solved by the KSP routines provided by the PETSc package [2]. Specifically, the GMRES routine with the Block Jacobi preconditioning was found to be reliable and efficient, among other combinations. The PETSc package offers sparse-matrix solvers in a parallel computing environment. In our 20-node 3.2 GHz Pentium 4 cluster system, it takes about 4 minutes to solve the above nonlinear Poisson’s equation with 6–7 N steps, for the device discretized by 100 × 100 × 100meshes. A linearized version of Eq. (40):    1 q ∇ 2 φk = − ND − nk3D 1 + , (43) (φk − φk−1 )  kB T was also implemented in our simulation, which is computationally efficient because it does not need the Newton–Raphson iteration steps. However, for unusually high drain and/or voltage biases, the above linearized Poisson’s equation either converges poorly or fails, and non-linearized one in Eq. (40) gives more reliable results. The linearized or non-linearized Poisson equation is solved self-consistently with the NEGF transport equations which will be described in the next subsection. Fig. 5 shows the convergence behavior of the Poisson-NEGF iterations where it can be seen that the self-consistent iterations converge at about 20 steps.

1068

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

3.4. Solution of 1D NEGF equations The 1D Green’s function Gm (E) in Eq. (8) was obtained by directly inverting the matrix: ⎡ ⎤ t 0 ... d1,1 (E) ⎢ ⎥ t 0 d2,2 (E) ⎢ t ⎥ ⎢ ⎥ ⎢ ⎥ . . . . −1 .. .. .. .. Gm (E) = ⎢ ⎥, ⎢ ⎥ ⎢ ⎥ t dNx −1,Nx −1 (E) t ⎣ 0 ⎦ ... 0 t dNx ,Nx (E)

(44)

where

⎧ if n = 1 ⎪ ⎨ E − (2t + Em (1) + ΣS (E)) if n = 2, . . . , Nx − 1 , dn,n (E) = E − (2t + Em (n)) ⎪ ⎩ E − (2t + Em (Nx ) + ΣD (E)) if n = Nx

(45)

where Nx is the number of meshes in the x direction and t≡

h ¯2 , 2m∗x x2

(46)

where x is the 1D mesh spacing. The self energies are given by

 1 ΣS (E) = E − 2t − Em (1) − i (4t + Em (1) − E)(E − Em (1)) , (47) 2

 1 ΣD (E) = (48) E − 2t − Em (Nx ) − i (4t + Em (Nx ) − E)(E − Em (Nx )) . 2 Note that only the very first and very last elements of the matrix {Gm } are needed. Numerically, the standard routine such as zgtsv() of LAPACK was employed to solve the tridiagonal linear equations. In the NEGF formalism, integration over energy as in Eqs. (10) and (16) must be performed to obtain the charge densities and drain current. Using simple uniform mesh for the numerical integration is either inefficient in terms of computation time or results in poor evaluation of the integration, because the charge density or transmission probability often yield sharp peaks as a function of energy, especially when the quantum phenomena such as tunneling and resonance come into play importantly. In our simulation, we have devised a simple algorithm for an iterative evaluation of the integration as follows. To be specific, let us consider the energy integration for current in a subband mode m as given by Eq. (16). Let us denote I (k) as the total current (for the mode m) at iteration step k, which is evaluated using a uniform mesh of size E/2k , where E = Emax − Emin and Emin and Emax are the minimum and maximum energies of the integration range, respectively. At the next step k + 1, the mesh size becomes half of the previous step, E/2k+1 , and the total current I (k+1) is given by, using a simple trapezoidal integration method: 1 (k) E  (k+1) I + k+1 I(Ei ), 2 2 k+1

I (k+1) =

(49)

i=1

where I(E) =

2q T (E)(fS (E) − fD (E)). h (k+1)

(50) (k+1)

In Eq. (49), Ei are middle points of the energy points used in the integration of the previous step: Ei = (k) (k) (Ei−1 + Ei )/2. To handle possible sharp peaks in I(E) with respect to E and achieve the computational efficiency at (k+1) ) at step k + 1 only when both the same time, after some pre-assigned step, for example k = 8, we evaluated I(Ei (k) (k) (k) I(Ei−1 ) and I(Ei ) of step k are higher than some fraction of the maximum current value Imax at step k. Otherwise, we

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

1069

Fig. 6. Drain current (Id ) versus drain voltage (Vd ) characteristic of a NWFET. The gate voltage is fixed at 0.6 V. Upper and lower curves represent the total current and the tunneling component of the total current. (k+1)

simply took the average of the two to set the new value in the middle : i.e., I(Ei the iteration stops when the total current converges.

(k)

(k)

) = (I(Ei−1 ) + I(Ei ))/2. And

4. Simulation results In this Section, we briefly show some simulation results that were generated using the simulator that we have developed. Simulated device is a nanowire FET with a square cross-section of width of 5 nm and channel length of 10 nm. The gates are all around the channel, the oxide thickness is 1 nm, and the channel and S/D doping densities are 1010 cm−3 and 1020 cm−3 , respectively. The silicon is (1 0 0) oriented and the mid-gap gate work function is used. Fig. 6 shows the drain current together with the tunneling component of the current when the gate voltage is fixed at 0.6 V. The current-voltage characteristic resembles that of the planar MOSFET, namely, initial rapid rising of the current followed by a saturation, but due to the fact that the channel length is short, the drain current is not as well saturated and rather it increases almost linearly. It can be seen in the figure that the tunneling component of the current contributes to the linear increase considerably. The conduction band edge profiles for different drain voltages are shown in Fig. 7 where the gate voltage is fixed at 0.6 V. The figure shows that as the drain voltage is increased, the height of the source/channel barrier does not change much whereas the width of the barrier becomes shorter, which explains the role of the tunneling component of the current in Fig. 6.

Fig. 7. Conduction band edge profiles related to Fig. 6 as the drain voltage changes from 0 (uppermost) to 0.5 V (lowermost) in a step of 0.1 V.

1070

M. Shin / Mathematics and Computers in Simulation 79 (2008) 1060–1070

5. Conclusions We have described in detail the simulation methods and techniques for the three-dimensional quantum simulation of the multigate nanowire field effect transistors in the ballistic transport regime. In particular, we have given a detailed description of the solution of the two-dimensional Schr¨odinger equation both in the k-space and in the product space, latter of which we have developed for the simulations of nanowire field effect transistors. Numerical methods for the solution of the three-dimensional Poisson’s equation and the one-dimensional transport equation in the non-equilibrium Green’s function framework have been also presented. Acknowledgements This research was supported by the Ministry of Information and Communication, Korea, under the Information Technology Research Center (ITRC) support program supervised by the Institute of Information Technology Assessment (IITA). References [1] A. Abramo, A. Cardin, L. Selmi, E. Sangiorgi, Two-dimensional quantum mechanical simulation of charge distribution in silicon MOSFETs, IEEE Trans. Electron Dev. 27 (2000) 1858–1863. [2] S. Balay, K. Buschelman, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith, H. Zhang, PETSc Web Page (http://www.mcs.anl.gov/petsc), 2001. [3] M. Bescond, K. Nehari, J.L. Autran, N. Cavassilas, D. Munteanu, M. Lannoo, 3D quantum modeling and simulation of multiple-gate nanowire MOSFETs, IEDM Tech. Digest (2004) 617–620. [4] Y. Cui, Z. Zhong, D. Wang, W. Wang, M. Lieber, High performance silicon nanowire field effect transistors, Nano Lett. 3 (2003) 149–152. [5] S. Datta, Nanoscale device modeling: the Greens function method, Superlatt. Microstruct. 28 (2000) 253–278. [6] H.K. Gummel, A self-consistent iterative scheme for one-dimensional steady state transistor calculations, IEEE Trans. Electron Dev. 11 (1964) 455–465. [7] J.T. Park, J.P. Colinge, Multiple-gate SOI MOSFETs: device design guidelines, IEEE Trans. Electron Dev. 49 (2002) 2222–2229. [8] H.-S. Philip Wong, Beyond the conventional transistor, Solid State Electron. 49 (2005) 755–762. [9] J. Wang, E. Polizzi, M. Lundstrom, A three-dimensional quantum simulation of silicon nanowire transistors with the effective-mass approximation, J. Appl. Phys. 96 (2004) 2192–2203.