Ballistic Calculation Of Non Equilibrium Green's Function In Nanoscale Devices Using Finite Element Method

IOP PUBLISHING

JOURNAL OF PHYSICS D: APPLIED PHYSICS

J. Phys. D: Appl. Phys. 42 (2009) 105109 (11pp)

doi:10.1088/0022-3727/42/10/105109

Ballistic calculation of nonequilibrium Green’s function in nanoscale devices using finite element method O Kurniawan, P Bai and E Li Computational Electronics and Photonics Group, Institute of High Performance Computing, Singapore 138632, Singapore E-mail: [email protected]

Received 6 March 2009, in final form 8 April 2009 Published 29 April 2009 Online at stacks.iop.org/JPhysD/42/105109 Abstract A ballistic calculation of a full quantum mechanical system is presented to study 2D nanoscale devices. The simulation uses the nonequilibrium Green’s function (NEGF) approach to calculate the transport properties of the devices. While most available software uses the finite difference discretization technique, our work opts to formulate the NEGF calculation using the finite element method (FEM). In calculating a ballistic device, the FEM gives some advantages. In the FEM, the floating boundary condition for ballistic devices is satisfied naturally. This paper gives a detailed finite element formulation of the NEGF calculation applied to a double-gate MOSFET device with a channel length of 10 nm and a body thickness of 3 nm. The potential, electron density, Fermi functions integrated over the transverse energy, local density of states and the transmission coefficient of the device have been studied. We found that the transmission coefficient is significantly affected by the top of the barrier between the source and the channel, which in turn depends on the gate control. This supports the claim that ballistic devices can be modelled by the transport properties at the top of the barrier. Hence, the full quantum mechanical calculation presented here confirms the theory of ballistic transport in nanoscale devices. (Some figures in this article are in colour only in the electronic version)

compute Green’s function matrices. Once Green’s function matrices are computed, the electron density and the current at the terminal can be obtained immediately. Moreover, it is possible to include phase breaking phenomena such as phonon scattering or photon–electron interaction through the self-energy matrices within this NEGF framework. This makes NEGF an important tool to investigate the thermal and optical properties of nanoscale devices. The NEGF formalism can be applied to finite structures following Caroli et al [10–13]. The domain is divided into the left contact, the right contact and the channel. It was shown that the current can be expressed in terms of the distribution functions in the contacts and the local properties of the channel such as the occupation and the density of states. For ballistic transport, the current can be written in the usual Landauer formula. Hence, the ballistic transport of nanoscale devices can be modelled using the NEGF framework.

1. Introduction New theoretical tools are always needed with the advances in technological progress. The progress of device scaling for transistors has reached the nanometre scale. Intel has begun the production of 45 nm transistors, and in 2009 will start producing chips with 32 nm process technology [1]. As the effort to scale the transistor continues, many predicted that the transistor could operate in the ballistic regime [2–4]. In this nanometre scale and ballistic regime, quantum phenomena occurs. It is interesting, therefore, to study the device behaviour from the perspective of quantum mechanics. A full quantum mechanical simulation of device transport is commonly done using the nonequilibrium Green’s function (NEGF) calculation [5–7]. The formulation of NEGF which we followed was proposed by Keldysh, Kadanoff and Baym [8, 9]. The main calculation of this NEGF method is to 0022-3727/09/105109+11$30.00

1

© 2009 IOP Publishing Ltd

Printed in the UK

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

One way to scale the transistor to nanometre size is to use a double-gate MOSFET. Some works have been done to simulate the double-gate MOSFET within the NEGF framework. Ren et al studied the scaling of double-gate MOSFET with gate length 10 nm [14]. The simulation results were used to identify issues that were needed to be addressed in order to fabricate nanoscale transistors. Svizhenko et al also simulated the device under scattering and compared their results with Medici, a 2D device simulator [15]. Their results show that the quantum calculation differs significantly from the semi-classical calculation implemented in Medici. They developed a recursive algorithm to optimize the simulations. The reason for this optimization is that their approach uses real-space discretization which requires huge matrix size and computational time. To reduce the size of the matrix and the computational time, Ren et al used subband or mode space technique in their software which is named NanoMOS [16]. The same model as NanoMOS has been used in some literature to analyse the double-gate MOSFET transistors [17]. For example, Venugopal et al studied the geometry of the channel access in double-gate MOSFET transistors [18]. In this technique, the Schr¨odinger equation is solved in the plane perpendicular to the transport direction, then the subband energies and eigenfunctions are calculated. The wave function is expanded in the subbands and the resulting transport equation is simply a one-dimensional equation. NanoMOS further simplified the calculation by using the uncoupled mode space technique. This means that the different modes along the transport direction are not coupled. Luisier et al studied the effect of the coupling between the modes and concluded that the effect of correlation between the modes is negligible if the channel is uniform [19]. Instead of using a simple parabolic band structure as the one in NanoMOS, Xia et al used a full complex band structure [20]. The empirical tight-binding Hamiltonian was used in the calculation and the result was compared with the one using effective mass approximation. The results, however, show that the effective mass approximation is sufficient for 10 nm Si MOSFET channel length simulation. The effect of band structure is prominent when the channel length is scaled further to about 5 nm or when the transistor is in the OFF state. Since in this work we are dealing with a 10 nm double-gate MOSFET in the ON state, we use the effective mass approximation in our calculation. Most of the previous works, such as NanoMOS, use the finite difference method in their calculations. This current work, on the other hand, uses the finite element method (FEM) to discretize the domain. The advantage of FEM is that it can be extended to more exotic device geometries and to higher dimensions without reformulating the equations. It is also possible to control the accuracy of approximation by refining the mesh. Besides that, nonuniform mesh spacing can be used without additional difficulties. A closer look also reveals that the self-energy matrix of the semi-infinite lead in the finite element formulation differs from that in the finite difference formulation (cf equation (9) in [18] and equation (44) in this paper). Moreover, FEM formulation

Figure 1. The schematics of a double-gate MOSFET device.

satisfies Neumann boundary condition automatically and implicitly. This boundary condition must be applied at the source and drain contacts of ballistic devices [4, 21]. This boundary condition at the source and drain is usually called the floating boundary condition, which is implemented by taking the gradient of the potential normal to the contacts to be zero.
= 0, xˆ · ∇V

(1)

where V is the electrostatic potential. This floating boundary condition is required in ballistic devices to ensure the charge neutrality in the source and the drain region. The authors feel that there is a lack of available finite element formulation for the NEGF. Havu et al gave a formulation for finite element NEGF [22]. The formulation, however, is given for real space basis. In this paper, we present the finite element formulation of the NEGF calculation in mode space basis. By using mode space basis the 2D transport equation is simplified into 1D and the size of the matrix reduces significantly. Moreover, we use a first-order perturbation theory to speed up the calculation of the potential. This technique has been shown to be sufficiently accurate for uniform potential in the channel. Detailed calculation procedures and mathematical expressions are provided so that the calculation can be performed on a personal computer. The implemented formulation includes the nonlinear Poisson equation, discretized Schr¨odinger equation and the transport equation of the NEGF. The developed method has been used to study the selfconsistent potential, the electron density, the Fermi function integrated over the transverse energy at the source and the drain, the local density of states (LDOS), and the transmission coefficient of a ballistic double-gate MOSFET device. We found that the transmission coefficient is affected by the top of the barrier between the source and the channel. This maximum barrier height depends on the gate control of the device. This result supports the statement in [4] which claims that the top of the barrier is significant in determining the transport of ballistic devices. Hence the full quantum mechanical calculation of ballistic devices presented here has been utilized to support the theory of ballistic transport. Furthermore, the numerical calculation described in this work can also be used for other ballistic nanoscale devices.

2. Double-gate MOSFET simulation structure The structure that we simulate is shown in figure 1. The silicon body thickness is 3 nm while the oxide thickness on the top 2

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

where C is just a normalization constant and is usually taken to be the square root of a length, area or volume. The above equation means that the wave function is a plane wave in the z direction. Substituting this into the 3D Schr¨odinger equation results in the equation for the x–y confined wave function. −

(4)

where El is the longitudinal energy in the x–y direction, and Ekz is the energy in the z direction. In this equation U (x, y) is the potential energy in the x–y plane, m∗x and m∗y are the effective mass in the x and y direction, respectively, and (x, y) is the wave function in the confined x–y plane. Now we expand (x, y) in the subband or mode eigenfunctions, thus comes the term mode space.
φ n (x)ξ n (y; x), (5) (x, y) =

Figure 2. Self-consistent calculation of the Poisson equation and the transport equation.

and at the bottom is 1 nm. The source and drain length is 5 nm and the gate length is 10 nm. The device domain is discretized using the FEM. Gmsh was used to generate the 2D mesh of the device [23]. The Dirichlet boundary nodes are located on the top and at the bottom gate electrodes. The potentials at the source and the drain contacts are left to float freely to ensure the charge neutrality in the device [4, 21]. In other words, it satisfies the Neumann boundary conditions as discussed previously.

n

where n is the index for the modes, and ξ n (y; x) is obtained by solving the Schr¨odinger equation in the plane normal to the transport direction, which in this case is along y. −

3. Finite element formulation for the NEGF

h ¯ 2 d2 ξ(y; x0 ) n + U (y; x0 )ξ(y; x0 ) = Esub (x0 )ξ n (y; x0 ). 2m∗y dy 2 (6)

Substituting equation (5) into equation (6), and then multiplying with ξ m (y; x) and integrating over y, we get 
∂2 ∂ h ¯2 bmn (x) φ n (x) amm (x) 2 φ n (x) + 2 − 2 ∂x ∂x n 
m + cmn (x)φ n (x) + Esub (x)φ m (x) = El φ m (x), (7)

For the NEGF calculation, we use the fast uncoupled mode space technique [24]. In this mode space technique the domain is discretized in its transverse and longitudinal directions. The modes due to the confinement in the transverse direction are calculated by solving the Schr¨odinger equation perpendicular to the transport direction. The main calculation of the NEGF consists of a selfconsistent computation of the electrostatic equation and the transport equation (figure 2). The detailed calculations are presented in the following sections starting with the transport equation. From the transport equation, the electron density can be obtained and substituted into the Poisson equation. The Poisson equation computes the self-consistent potential, which in turn is substituted back into the transport equation.

n

where 1 a mn (x) = ∗ mx bmn (x) =

1 m∗x

1 cmn (x) = ∗ mx

3.1. Fast uncoupled mode space

 ξ m (y; x)ξ n (y; x) dy =  

y

(8)

ξ m (y; x)

∂ n ξ (y; x) dy, ∂x

(9)

ξ m (y; x)

∂2 n ξ (y; x) dy. ∂x 2

(10)

y

y

1 δm,n , m∗x

The Hamiltonian can be written as  2  h ¯ ∂2 m (x) + Esub hmn = δm,n − amm (x) 2 ∂x

In this section, the formulation of the transport equation is presented. The calculation for the transport equation is based on the fast uncoupled mode space technique as described in [24]. We start from the 3D time-independent Schr¨odinger equation, (2) H3D
(x, y, z) = E
(x, y, z),

h ¯2 ∂ cmn (x) − h (11) ¯ 2 bmn (x) . 2 ∂x In the uncoupled mode space technique, we assume that there is no interaction between the modes. The eigenfunctions ξ n (y; x0) are assumed to be the same throughout the device: −

where H3D is the Hamiltonian of the system, E is the eigenvalue and
(x, y, z) is the wave function. Since the z dimension is assumed to be large, we expand the solution for the wave function as
(x, y, z) = C(x, y) exp(ikz z),

h ¯ 2 ∂ 2 (x, y) h ¯ 2 ∂ 2 (x, y) − + U (x, y)(x, y) 2m∗x ∂x 2 2m∗y ∂y 2 = (E − Ekz )(x, y) = El (x, y),

(3) 3

ξ m (y; x0 ) = ξ m (y; x0 ),

(12)

∂ m ξ (y; x0 ) = 0. ∂x

(13)

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

With this assumption, bmn (x) = cmn (x) = 0 for all m and n. Moreover, hmn = 0 for m = n. Hence, it can be seen that the interaction between the modes are zero. The Hamiltonian becomes block diagonal and is given by hmm = −

h ¯2 ∂2 m amm (x) + Esub (x). 2 ∂x

considered similar to the density of states. The spectral density matrix is [7]

(14)

 ∂ h ¯ m + Esub (x) φ m (x) = El φ m (x), − amm (x) 2 ∂x 2

(15)

This eigenvalue problem can be written in a matrix form as shown below: (16)

where S m is the overlap matrix. m (x) is the subband energy along In the Hamiltonian, Esub x which is obtained by solving the Schr¨odinger equation in the confined direction, i.e. equation (6). For the fast uncoupled mode space, this equation is computed only once for the whole device. After we obtain the potential energy, we calculate the average values  U (y) = (1/Lx )

where Lx is the total length of the device. This average potential is then substituted into equation (6) to obtain the m (x0 ) averaged subband eigenvalues and eigenvectors, i.e. Esub m and ξ (y; x0). The subband energies are then calculated from the first-order stationary perturbation theory  m m Esub (x) = Esub (x0 ) + U (x, y)|ξ m (y; x0 )|2 dy −

U (x, y)|ξ m (y; x0 )| dy.

(22)

Dm (E) = i(Dm (E) − Dm† (E))

(23)

m m 2 nm 2D (x, y) = n1D (x)|ξ (y; x)| .

y

2

Sm (E) = i(Sm (E) − Sm† (E)),

is the Fermi function that includes the integration of the k space along the z direction. In the calculation of the density of states, we include the factor of two to account for the symmetry of the unprimed valleys in silicon. In this equation m∗z is the effective mass in the z direction, and F−1/2 is the complete Fermi–Dirac integral of order −1/2. This equation includes the spin degeneracy. The 2D electron density is obtained by multiplying the 1D density with the wave function in the transverse direction:

(17)

0



(21)

where µS and µD are the electrochemical potentials of the source and drain contacts, a is the size of the unit cell and    ∗  2mz kB T 1/2 µ−E F−1/2 (25) f1D (E) = 2 kB T πh ¯2

Lx

U (x, y) dx,

m m m† Am D (E) = G (E)D (E)G (E),

are the broadening matrix due to the source and drain contacts. The electrons in the device come from the injection of electrons in the source, which try to achieve equilibrium with the source Fermi level, as well as the electrons in the drain, which try to achieve equilibrium with the drain Fermi level. The 1D electron density can then be obtained from the diagonal elements of  ∞ 1 m f1D,S (µS − El )Am = nm 1D,x S (El )S 2π a −∞

m + f1D,D (µD − El )Am (24) dEl , D (El )S

2

hmm φ m = El S m φ m ,

(20)

where

When this equation is discretized using the FEM, it will give us the Hamiltonian matrix and the overlap matrix. The Schr¨odinger equation is 

m m m† Am S (E) = G (E)S (E)G (E),

(18)

(26)

This is for one subband. The total electron density is calculated by summing the whole subbands. Usually accurate results can be obtained with less than five subbands [24]. Note that the two partial spectral functions in equation (20) add up to give the total spectral function:

y

The problem has been simplified into a 1D transport equation. Now we are ready to calculate the Green’s function. For a given mode or subband, the retarded Green’s function is calculated from  −1 . (19) Gm (E) = ES m − hmm − Sm (E) − Dm (E)

m Am (E) = Am S (E) + AD (E).

(27)

The LDOS at position r along the channel is obtained from Note that this Green’s function matrix is a function of energy. The superscript m indicates the index for the modes. All the terms in this equation are matrices of size Nx × Nx . The matrix S m is the overlap matrix in equation (16), and the Sm (E) and Dm (E) are the self-energy matrices to account for the interaction with the source and drain contacts. The selfenergy matrices will be discussed later. Once the Green’s function matrix is calculated, the electron density and the terminal currents can be easily obtained. To obtain the electron density, we first need to calculate the spectral density. The spectral density can be

D(r; E) = Am (E)r,r S m /(2π ).

(28)

Now, the current is calculated from  +∞ I = (q/2π h ¯) T (E) [f1D (µS − E) − f1D (µD − E)] dE −∞

(29)

and T (E) =

M
m=1

4

T m (E)

(30)

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

is the transmission function and is obtained by summing the transmission function for each subband:  T m (E) = Trace Sm (E)Gm (E)Dm (E)Gm† (E) . (31)

(a)

In the calculation of electron density and current (cf equations (24) and (29)), an integration over the energy axis must be performed. The integration for the electron density requires further discussion. The reason is that the poles of the Green’s functions are located at the real energy axis. To overcome this, the integration is normally performed by adding a constant imaginary component. However, when the integration is done too close to the real energy axis, a very fine mesh is required. On the other hand, errors are introduced when the integration is performed too far from the real energy axis. Hence, this technique is not recommended particularly for self-consistency calculation. For this reason, most people who use density functional theory with NEGF commonly solve the problem by using contour integration and separating the calculation of the electron density to reduce the errors [25–30].

(b)

Figure 3. (a) Finite element of a 1D geometry and (b) linear shape functions.

when solving equation (32) or (6), but it becomes dx when solving equation (33) or (15). For the Hamiltonian of equation (32), αe = −

h ¯2 , 2m∗y

β e = U, 3.2. Finite element discretization of Schr¨odinger and NEGF equation

HT ξ −

= 0,

HL φ − El SL φ = 0

(36)

while for the Hamiltonian of equation (33),

The previous section presents the calculation of the transport equation. The transport equation is solved in the mode space approach. In this approach, the Schr¨odinger equation is solved in the confinement direction to obtain the subband energies (cf equation (6)). Once the subband energies are obtained, it can be used to build the Hamiltonian matrix of the device (equations (14) and (15)). This Hamiltonian matrix is then substituted into equation (19) to calculate Green’s function. With Green’s function evaluated at every energy, the electron density (cf equation (24)) and the terminal current (cf equation (29)) can be computed. In order to proceed, we need to discretize equations (6) and (15). We discretize the two equations using the FEM. The two equations are eigenvalue problems, which can be written in matrix form as n Esub ST ξ

(35)

αe = −

h ¯2 , 2m∗x

m (x). β e = Esub

(37) (38)

It is important to note that we only solve equation (32) or (6), but not equation (15). This is because, in the longitudinal direction, we solve the transport equation through the NEGF method. Nevertheless, the overlap matrix is required in both equations. The element of the overlap matrix is given by  Sije =

 Nie Nje dk



(39)

and dk can be either dx or dy depending on the equations to be solved. If we use the piecewise continuous linear shape function (figure 3), the integral can be solved analytically to give

(32) (33)

for equations (6) and (15), respectively, where HT is the Hamiltonian in the transverse direction, ST is the overlap matrix in the same direction, and HL is the Hamiltonian equation (14) in the longitudinal direction and SL is its overlap matrix. Since these two equations are basically a 1D problem, the finite element formulation for each element is given simply by [31]    e ∂N e j e e ∂Ni e e e Hij = (34) α + β Ni Nj dk, ∂k ∂k

e = H11

e αe el e + β , = H22 le 3

e =− H12

where Hije can be the elements of either HT or HL . The differential dk in the integral of equation (34) becomes dy

e αe el e = H21 + β , le 6

(40)

(41)

e S11 =

le e = S22 , 3

(42)

e S12 =

le e . = S21 6

(43)

The only term left for us to calculate the Green’s function is the self-energy matrices (cf equation (19)). 5

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

3.3. Self-energy matrix of the infinite lead The self-energy matrix for an infinite lead was available in the literature for finite difference discretization [7, 18]. The expression for finite element discretization using effective mass approximation was given by Polizzi and Datta [32]. It is shown that the self-energy matrix is independent of the discretization step. The expression for the matrix element is given by h ¯2 [L ]ij = − ∗ ikL δi∈s δj ∈s , (44) 2mL

(a)

(b)

where the subscript L indicates the infinite lead which can either be the source (S ) or the drain (D ), and s is the node index at the interface between the contacts and the device. The self-energy matrix is nonzero only in the interface nodes. The wave vector kL is obtained from the dispersion relation E=

h ¯ 2 kL2 . 2m∗L

Figure 4. (a) Triangular finite element of a 2D geometry and (b) linear shape functions.

Note that U is the potential energy and has a unit of one electron volt. The system of equation can be written in the matrix form: KU = B, (48)

(45)

So for a given energy, we can calculate kL and then substitute into equation (44) to get the self-energy matrix. This expression does not depend on the discretization spacing as the one in [24] which uses the same effective mass approximation. We further implemented the technique proposed by Lake et al to include an energy dependent imaginary potential in the self-energy of the contacts, i.e. ıη [33]. Physically, this potential represents the scattering induced broadening and acts as an approximation for the imaginary part of the retarded selfenergy. The imaginary potential is nonzero only in the contacts and zero in the device. Including this potential is crucial to take into account the injection from the source that lies below the continuum [34]. With all this information, the transport equation can be computed to give the electron density. The electron density is then substituted into the Poisson equation to obtain the electrostatic solution.

where the K matrix and the B column vector are assembled from each finite element and U is the solution of potential energy at every node.   e ∂N e e ∂N e  j j e e ∂Ni e ∂Ni Kij = + r dx dy (49)

r ∂x ∂x ∂y ∂y Bie =

To obtain the electrostatic solution, we solve the Poisson equation assuming an n-type semiconductor: q (N + − n),

0 D

q  e e n − ND+ Nie dx dy,

0

(50)

where i, j = 1, 2, 3. The superscript e is to indicate the calculation at the elements, and Nie is the interpolation function or the shape function. In this work we use a piecewise continuous linear triangulation shape function, as shown in figure 4. The shape function is given in [31] and provided in the appendix. Solving the integration gives  K e = ( re /4Ae ) be be T + ce ce T (51)

3.4. Nonlinear Poisson equation

∇( r ∇V ) = −



and

and be and ce are column vectors with the exact formula given in [31, p 80] and in the appendix. The superscript T indicates the transpose operation and the term A is the area of the 2D finite element. For the linear triangulation finite element with three nodes, the size of matrix K e will be 3 × 3 and the size of the column vectors be and ce is 3 × 1. Similarly,

(46)

where V = −U/q is the electrostatic potential, U is the potential energy in joule and q is the unit charge. The righthand side is the net charge, which is computed by taking the difference between the positively charged donor concentration, ND+ and the free carrier electron charge n. The material permittivity, = 0 r , is allowed to vary with the regions. In this case, r is the relative permittivity and 0 is the free space permittivity. For 2D, the Poisson equation can be written as [31]      ∂ ∂U ∂U ∂ q  − n − ND+ . (47)

r −

r = ∂x ∂x ∂y ∂y

0

e

B e = (Ae /3)(q/ 0 )(ne − ND+ )

(52)

and the size of column vector B e is 3 × 1. After imposing the boundary conditions at the gate electrodes, we can then solve the system of equations to obtain the potential energy. The boundary conditions at the gate electrodes are calculated from the gate bias and the bend bending, i.e. UG = −VG + G − ψSi (53) 6

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

where VG is the gate bias, G is the gate metal work function and
Si is the electron affinity in the silicon. The symbols  and ψ should not be confused with the wave functions in the previous sections. As can be seen from figure 2, the potential is used to solve the transport equation, which will give us the electron density. The electron density is then substituted back into the Poisson equation to correct the value of the potential. This is a nonlinear Poisson equation. We solve this equation using the technique presented in [16, 35]. It is shown that this technique improves the convergence of the iteration. The technique basically involves two iterations. The first one is to solve the nonlinear Poisson equation using the Newton–Raphson method, and the other one is to iterate the self-consistent equation as shown in figure 2. The former is the inner iteration while the latter is the outer iteration. The inner iteration is described as follows. Once we obtain the electron density, we do not substitute it immediately into the Poisson equation. On the other hand, we calculate the quasi-Fermi level   kB T −1 n2D , (54) F1/2 Fn = Uold + q NC

assembled from all the elements, and it can be shown that the matrix for one element is given by Jf (Uk )eij = Kije −

Note also that the property of the complete Fermi–Dirac integral as shown below has been utilized [36]: d Fk (x) = Fk−1 (x). dx

4. Simulation results and discussion The simulation structure is shown in figure 1. For simplicity we use a 0.25 nm uniform grid spacing along the transport direction. In the transverse direction, the finite elements are generated automatically with a 0.1 nm characteristic length. The size of the element, then, is computed by linearly interpolating these characteristic lengths on the initial mesh according to the Gmsh manual [23]. Though we will only compute using a uniform mesh, the formulation given in this work is the same even if we use the nonuniform grid spacing along the transport direction. The reason is that the formulation is given for each element of the mesh. And for each element, the length is computed in the formulation. This is the advantage of the FEM. The relative permittivity for the silicon body is set to 11.7, while the permittivity for the oxide is set to 3.9. The conduction band of the oxide is 3.34 eV higher than that in the silicon body. For the boundary conditions in the gate electrode, we use 4.05 eV as the electron affinity in the silicon and 4.25 eV as the metal work function for the two gate electrodes. The source and the drain doping concentrations are set to 1 × 1020 cm−3 , while the region in the channel is assumed to have a doping concentration of 1×1016 cm−3 . As can be noticed, the channel is nearly undoped. It is common to have an intrinsic channel or low doping concentrations for the channel of a doublegate device [16, 20, 37–42]. The reason for this is to reduce scatterings along the channel. In this work, we use the bulk value of silicon effective mass. The longitudinal effective mass is set to ml = 0.98m0 , and the transverse effective mass is set to mt = 0.18m0 . The gate electrodes were biased symmetrically, which means both the top and the bottom gates are biased at the same

where now F1/2 is the complete Fermi–Dirac integral of order 1/2. This nonlinear equation is solved using the Newton– Raphson method by rewriting it in the form of (57)

where f (U ) is obtained by moving the right-hand side of equation (56) to the left-hand side. If we discretize using the FEM, the above equation will become a set of system of equations. And the solution of U for the next iteration is calculated from (58)

where sk is obtained by solving the following linear equation: Jf (Uk )sk = −f (Uk )

(62)

After solving the nonlinear Poisson equation, we obtain the energy potential U . This energy potential is then substituted into the transport equation to obtain the new electron density. This is the outer iteration. The iteration continues until the solution converges to the specified tolerance.

Once the quasi-Fermi level is obtained, it is substituted into the nonlinear Poisson equation (cf equations (46) and (54)):     Fn − Unew qNC N+ F1/2 −∇( r ∇Unew ) = − D , (56)

0 kB T NC

Uk+1 = Uk + sk ,

(60)

where the expression for B e is given in equation (52). Note that B e contains the term electron density (cf equation (55)) which is a function of the potential. Hence, we need the expression for the derivative of the electron density,   Fn − U NC ∂n2D =− F−1/2 . (61) ∂U kB T kB T

where kB T /q = 0.00259 eV at room temperature, n2D is the 2D electron density, and NC is the effective density of states in the conduction band. At room temperature, NC = 2.82 × 1025 m−3 . The function F−1 1/2 is the inverse of the complete Fermi–Dirac integral of order 1/2. In the calculation of Fn , we use the previous solution of the potential U . The above equation can be rearranged as   Fn − Unew n2D = NC F1/2 . (55) kB T

f (U ) = 0,

∂Bie , ∂Uje

(59)

and Jf (Uk ) is the Jacobian matrix of f (U ). The subscript k indicates the number of iterations. The Jacobian matrix is also 7

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

is about −0.15 V. This is due to the charge neutrality being satisfied in this region. The potential is allowed to float to any value it chooses to satisfy the charge neutrality. So in this case, the potential is lowered to increase the electron density back to 1 × 1020 cm−3 . The same conclusion can be observed in figure 6(c) as well. In this figure it can be seen that the electron density in the source and drain region does not change significantly even when they are biased. The electron density near the two contacts are preserved as seen from the plot. Figure 7 gives us the LDOS at the left end and at the right end of the device when the gate and the drain are biased at 0.4 V. The maximum peak of the LDOS at the right end of the device is lower than that at the left end by about VD . It can also be observed that the right end of the device has more available states. It is interesting to note that the LDOS at the left end of the device goes down to zero not at 0 eV, but rather, at energy near −0.2 eV. This is the self-consistent energy potential at the source as can be seen from figure 6(d). Similarly, the LDOS at the right end of the device goes down to zero at energy around −0.55 eV, which is the self-consistent potential energy at the drain (figure 6(d)). Now let us take a look at the current equation. At zero bias, we would expect no current flowing. This is seen from equation (29) where the difference of the two Fermi functions is zero. We plot the curve of these two integrated Fermi functions together with the transmission coefficient in figure 8(a). Since the difference between the two Fermi functions is zero, and there is no overlap region with the transmission coefficient curve, no current flows in the device. Note that the transmission curve will display a staircase and goes up to M, where M is the number of modes we use in the calculation. A staircase transmission coefficient is what is expected in a two-dimensionally confined quantum device. The zero current can also be understood from the band energies in figure 6(b). At zero gate bias, the barrier in the channel is high and blocks electron injection from source to drain. Applying a voltage bias on the gate electrode will lower the barrier. Similarly, applying a voltage bias on the drain will lower the barrier as well, although the effect is much smaller compared with on the gate. The plots of the Fermi functions and the transmission curve for VG = 0.4 V and VD = 0.4 V are shown in figure 8(b). The drain bias shifts both the Fermi function and the transmission coefficient. Since the source bias is fixed, curve F1 does not change. Curve F2 , on the other hand, shifts by VD = 0.4 V as can be seen from the plot. In figure 8(b), the transmission coefficient jumps from zero at around −0.1 eV, in contrast to figure 8(a) which is around 0.15 eV. These values turn out to be the top of the barrier between the source and the channel. Figure 6(b) shows that the top of the barrier of the first subband is around 0.15 eV. On the other hand, figure 6(d) shows that the top of the first subband energy is around −0.1 eV. To prove further, we calculate the first derivative of the transmission curve. If the jump from zero occurs at the top of the barrier, then we will be able to see a maxima in its first derivative at this energy level. Figure 9 shows the two maxima of the first derivative of the transmission curve coincide with

-2

N2D (m )

1e+18

1e+17

1e+16 -5

0

5

(a)

10 X (nm)

15

20

25

20

25

0 -0.1

E (eV)

-0.2 -0.3 -0.4 -0.5 -0.6 -5

(b)

0

5

10 X (nm)

15

Figure 5. Comparison with NanoMOS at bias VG = 0.4 V and VD = 0.4 V suggests a reasonable accuracy. (a) Electron density along the transport direction for single mode calculation and (b) self-consistent potential along the transport direction.

potential. The gates are biased from 0 to 0.4 V. The source electrode is grounded while the drain electrode is also biased from 0 to 0.4 V. The potentials, subband energies, electron density, LDOS and the transmission functions were calculated. The code was written using the Octave scripting language (compatible with Matlab), and can be run in a normal PC. The average time it takes to solve a bias point is about 5 min when it is run in a 3 GHz workstation. The performance is expected to be poor since the Octave interpreter performs quite slow when it consists a lot of for-loops. Porting the code to pure C++ will increase the performance. Before proceeding with the analysis, we verified our calculation results with NanoMOS [16]. In the simulation of NanoMOS, we used 0.3 nm grid spacing along the transport direction, and 0.1 nm grid spacing along the transverse direction. The doping concentration, device dimensions, and the bias steps are similar to our calculation. Though our software and NanoMOS use different discretization schemes and assumptions, the results as shown in figure 5 suggest that our calculation gives a reasonable accuracy. We are now ready to further analyse the double-gate device. Figures 6(a) and (b) show the self-consistent potential and the subband energy at zero bias. We can see that the potential satisfies Neumann boundary condition at the source and drain contacts. Even though we set the source and drain at zero bias potential, the calculated potential in figure 6(b) 8

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al 0.2 0.15

E (eV)

0.1 0.05 0 -0.05 -0.1 -0.15 n=1

-0.2

(a)

(b) 1e+19

VG = 0.4 V, VD = 0.4 V VG = 0.0 V, VD = 0.0 V

5

10 X (nm)

15

20

10 X (nm)

15

20

0

1e+18

-0.1

1e+17

-0.2 E (eV)

N2D (m-2)

0

1e+16

-0.3 -0.4

1e+15 -0.5

1e+14 0

5

10 X (nm)

(c)

15

n=1

-0.6

20

0

(d)

5

m Figure 6. (a) Self-consistent potential of single mode calculation at VG = 0 and VD = 0. (b) First subband energies Esub (x). The potential at the source and drain float to satisfy the charge neutrality. The potential of the first subband inside the channel is about 0.2 eV due to the difference between the metal work function and the electron affinity in the silicon. (c) Electron density in the transport direction of a single mode calculation at zero bias and at VG = 0.4 V and VD = 0.4 V. The electron density around the source and drain is maintained at around m 1 × 1020 cm−3 . (d) First subband energy Esub (x) at VG = 0.4 V and VD = 0.4 V. The energy of the top of the barrier is the point where the transmission curve jumps to unity in figure 8.

0.3

density of states, however, is the constant term in front of the complete Fermi–Dirac integral in equation (25). Thus, for a given device configuration, the only term affecting the integrated Fermi function curve are the source and drain biases. The other parameter in equation (29) is the transmission coefficient. The transmission coefficient is affected by the top of the barrier of the energy band in the channel as discussed previously. This top barrier height is determined by the electrostatic potential of the device which is obtained by solving the self-consistent Poisson equation. This shows that the current–voltage characteristic is determined by only two parameters. The first one is the electrode bias at the drain with respect to the source. The drain bias shifts the integrated Fermi function at the drain side. The second one is the energy level at the top of the barrier in the channel. This energy level shifts the transmission curve accordingly. Hence, we can conclude that the top of the barrier in the channel plays a significant role in the transport of ballistic device. This result supports the previous works which claim that the transport of a ballistic device can be modelled by simply considering the transport at the top of the barrier [4, 17, 38, 43].

Source Drain

0.2 0.1 E (eV)

0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

-1

D(E) (eV ) Figure 7. LDOS at the left end and the right end of the device at VG = 0.4 V and VD = 0.4 V. The LDOS maximum peak at the right end of the device is lower than at the left by about VD .

the top barrier of the energy subbands. In this plot, two modes were used in the calculation. Equation (29) shows that the terminal current depends on the integrated Fermi functions and the transmission coefficient. The integrated Fermi function depends on the source and the drain bias which affects the electrochemical potential µ. This Fermi function also depends on the confinement of the device which determines the density of states in the channel. The

5. Conclusion We have presented the finite element formulation of the NEGF and used it to simulate a double-gate MOSFET device. The 9

J. Phys. D: Appl. Phys. 42 (2009) 105109

1.2e+09

F1 F2

1e+09

1.4

0.8

6e+08

0.6 4e+08

0 -0.2 -0.1

-0.1

0.2

-0.15 -0.2 0

4

F1 F2

1.4

1 0.8 1e+09

0.6

5e+08 0.2 0 -0.2 0 E (eV)

16

20

The advantage of the NEGF is primarily its straightforward way of including the noncoherent transport in the formulation. The phonon and electron–photon interaction is included in the same manner as the contacts. Future works will focus on these phenomena, particularly on the electron– photon interaction. The integration of electronics and optics has shown to be one of the most promising paths in the future.

0.4

-0.4

12

Figure 9. Energy subbands and the first derivative of the transmission curve for two modes calculation. The (blue and green) dotted curves are the energy subbands, and the (red) solid curve is the first derivative of the transmission curve. The plot shows that the peak of the first derivative of the transmission curve coincides with the top of the barrier of the energy subbands. (Colour online.)

1.2

1.5e+09

0 -0.6

8 x

2e+09

(b)

n=1 n=2

E (eV) 2.5e+09

20

-0.05

0.4

0.1 0.2 0.3 0.4 0.5

(a)

16

0

0 0

0

0.05

E (eV)

1

Transmission

8e+08

2e+08

F(E) (1/m)

-4

dT/dE 4 8 12

0.1

1.2

Transmission

F(E) (1/m)

O Kurniawan et al

0.2

0.4

Figure 8. (a) Integrated Fermi function and transmission coefficient of a single mode calculation for VG = 0 and VD = 0, where F1 = f1D (µS − E) and F2 = f1D (µD − E) are the fermi functions integrated over the transverse energy. Since no bias is applied, the values of the two integrated Fermi functions are equal. The terminal current is the overlap region between the difference of the integrated Fermi function with the transmission coefficient, which in this case is zero. (b) Integrated Fermi function and transmission coefficient for VG = 0.4 V and VD = 0.4 V. The integrated Fermi curve at the drain is shifted by VD = 0.4 V. The transmission coefficient curve jumps at around −0.1 eV which is the top of the barrier in the channel.

Acknowledgments The authors would like to thank Dr Polizzi and Professor Datta for helpful discussions on self-energy in finite element formulation.

Appendix. Shape functions

FEM is preferable when considering the geometry as well as the boundary conditions in the source and the drain. Moreover, the accuracy can be controlled through mesh refinement. The developed work uses the piecewise continuous linear shape function for the finite element. This can be extended to higher order for better accuracy. Within this formulation, the selfconsistent potential, subband energies, electron density, LDOS and the transmission function of a nanoscale device can be calculated. The source code is made available and can be downloaded from [44]. The calculated transmission coefficient of the double-gate MOSFET device shows an interesting feature. The jump in the transmission coefficient curve is located at the top of the energy barrier between the source and the channel. The plot of the first derivative of the transmission curve shows that its maxima coincide with the energy level at the top barrier. This supports the claim that the ballistic device can be modelled using a simple transport on top of the barrier.

On the other hand, the shape function for 2D finite element is given by [31] Nie (x, y) = where i = 1, 2, 3 and   e e x2 y3 − y2e x3e   a e = x3e y1e − y3e x1e  , x1e y2e − y1e x2e   e x3 − x2e   ce = x1e − y3e  .

 1  e ai + bie x + cie y , e 2A 

y2e − y3e

(A.1)



  be = y3e − y1e  , y1e − y2e

(A.2)

x2e − y1e Note that xi , yi and zi are the x, y, z locations of the ith node in the element. 10

J. Phys. D: Appl. Phys. 42 (2009) 105109

O Kurniawan et al

[23] Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities http://geuz.org/gmsh [24] Wang J, Polizzi E and Lundstrom M 2004 J. Appl. Phys. 96 2192 [25] Taylor J, Guo H and Wang J 2001 Phys. Rev. B 63 245407 [26] Brandbyge M, Mozos J L, Ordejon P, Taylor J and Stokbro K 2002 Phys. Rev. B 65 165401 [27] Xue Y, Datta S and Ratner M A 2002 Chem. Phys. 281 151 [28] Zhang C, Zhang X-G, Krstic P S, Cheng H P, Butler W H and MacLaren J M 2004 Phys. Rev. B 69 134406 [29] Zhang J, Hou S, Li R, Qian Z, Han R, Shen Z, Zhao X and Xue Z 2005 Nanotechnology 16 3057 [30] Li R, Hou S, Qian Z, Shen Z, Zhao X and Xue Z 2007 Chem. Phys. 336 127 [31] Jin J 1993 The Finite Element Method in Electromagnetics (New York: Wiley) [32] Polizzi E and Datta S 2003 Proc. IEEE-NANO (San Francisco, CA) vol 1, p 40 [33] Lake R, Klimeck G, Bowen R C and Jovanovic D 1997 J. Appl. Phys. 81 7845 [34] Lake R, Klimeck G, Bowen R C, Jovanovic D, Sotirelis P and Frensley W R 1998 VLSI Des. 6 9 [35] Ren Z 2001 PhD Thesis Purdue University [36] Trellakis A, Galick A T, Pacelli A and Ravaioli U 1997 J. Appl. Phys. 81 7880 [37] Venugopal R, Ren Z, Datta S and Lundstrom M S 2002 J. Appl. Phys. 92 3730 [38] Rhew J-H, Ren Z and Lundstrom M S 2002 Solid-State Electron. 46 1899 [39] Guo J and Lundstrom M S 2002 IEEE Trans. Electron Devices 49 1897 [40] Laux S E, Kumar A and Fischetti M V 2004 J. Appl. Phys. 95 5545 [41] Martinez A, Bescond M, Barker J R, Svizhenko A, Anantram M P, Millar C and Asenov A 2007 IEEE Trans. Electron Devices 54 2213 [42] Liu Y, Neophytou N, Low T, Klimeck G and Lundstrom M S 2008 IEEE Trans. Electron Devices 55 866 [43] Lundstrom M 2003 IEDM Technical Digest p 3311 [44] FeMOS Finite Element Simulation of double-gate MOS transistor http://software.ihpc.a-star.edu.sg/projects/ FeMos/FeMos.php

References [1] Intel http://www.intel.com [2] Frank D J, Laux S E and Fischetti M V 1992 IEDM Technical Digest (San Francisco, CA) p 553 [3] Timp G et al 1999 IEDM Technical Digest (Washington, DC) p 55 [4] Rahman A, Guo J, Datta S and Lundstrom M S 2003 IEEE Trans. Electron Devices 50 1853 [5] Datta S 1989 Phys. Rev. B 40 5830 [6] Datta S 2000 Superlatt. Microstruct. 28 253 [7] Datta S 2005 Quantum Transport: Atom to Transistor (Cambridge, UK: Cambridge University Press) [8] Kadanoff L P and Baym G 1962 Quantum Statistical Mechanics (Reading, MA: Benjamin) [9] Keldysh L V 1965 Sov. Phys.—JETP 20 1018 [10] Caroli C, Combescot R, Nozieres P and Saint-James D 1971 J. Phys. C: Solid State Phys. 4 916 [11] Caroli C, Combescot R, Nozieres P and Saint-James D 1972 J. Phys. C: Solid State Phys. 5 21 [12] Meir Y and Wingreen N S 1992 Phys. Rev. Lett. 68 2512 [13] Datta S 1995 Electronic Transport in Mesoscopic Systems (Cambridge, UK: Cambridge University Press) [14] Ren Z, Venugopal R, Datta S and Lundstrom M 2001 IEDM Technical Digest (Washington, DC) p 541 [15] Svizhenko A, Anantram M P, Govindan T R, Biegel B and Venugopal R 2002 J. Appl. Phys. 91 2343 [16] Ren Z, Venugopal R, Goasguen S, Datta S and Lundstrom M S 2003 IEEE Trans. Electron Devices 50 1914 [17] Lundstrom M and Ren Z 2002 IEEE Trans. Electron Devices 49 133 [18] Venugopal R, Goasguen S, Datta S and Lundstrom M S 2004 J. Appl. Phys. 95 292 [19] Luisier M, Schenk A and Fichtner W 2006 J. Appl. Phys. 100 043713 [20] Xia T, Register L F and Banerjee S K 2003 IEEE Trans. Electron Devices 50 1511 [21] Venugopal R, Ren Z and Lundstrom M S 2003 IEEE Trans. Nanotechnol. 2 135 [22] Havu P, Havu V, Puska M J and Nieminen R M 2004 Phys. Rev. B 69 115325

11