Simulation of Quantum Transport in Nanoscale Devices

4 downloads 0 Views 4MB Size Report
p-type regions. The drain is biased at 3.0 V and the gate voltage is slightly below the threshold [Tau98]. ...... from the Brillion zone center along the. > ...... V. I. Kuznetsov, A. J. M. M. van Zutphen, H. R. Kerp, P. G. Vermont,. X. Pages, J. J. van ...
Ain Shams University Faculty of Engineering Electronics and Communication Engineering Department

Simulation of Quantum Transport in Nanoscale Devices

A THESIS Submitted in Partial Fulfillment of the Requirements for the Degree of Master of Science in Electrical Engineering By Yasser Mohammed Sabry Gad B.Sc. in Electrical Engineering (Electronics and Communication Engineering) Ain Shams University - 2005

Supervised by: Prof. Dr. Magdy Mahmoud Mohamed Ibrahim Ain Shams University

Prof. Dr. Hani Fikry Ragai Ain Shams University

Prof. Dr. Wael Fikry Farouk Ain Shams University Cairo, 2009

STATEMENT This thesis is submitted to Ain Shams University in partial fulfillment of the requirements for the degree of Master of Science in Electrical Engineering. The work included in the thesis was carried out by the author at the Electronics and Communication Engineering Department, Faculty of Engineering, Ain Shams University, Cairo, Egypt. No part of this thesis has been submitted for a degree or a qualification at any other university or institute.

Name: Yasser Mohammed Sabry Gad Signature: Date:

CURRICULUM VITAE

Name of Researcher

Yasser Mohammed Sabry Gad

Date of Birth

21/05/1982

Place of Birth

Kuwait

First University Degree

B.Sc. in Electrical Engineering

Name of University

Ain Shams University

Date of Degree

June 2005

ABSTRACT Yasser Mohammed Sabry Gad, “Simulation of Quantum Transport in Nanoscale Devices,” Master of Science dissertation, Ain Shams University, 2009. Rapid device scaling pushes the dimensions of the field effect transistors to the nanometer regime where quantum effects play an important role in determining the transistor characteristics. The Non-equilibrium Green’s function formalism (NEGF) provides a rigorous description of quantum transport in nanoscale devices. The RealSpace (RS) representation is the most accurate yet complex representation used in the NEGF. The geometry of fully-depleted double gate (DG) MOSFETs permits the use of a simpler quasi two-dimensional (2D) representation, the mode-space (MS) which is computationally efficient. This thesis addresses the simulation and design of the nanoscale DG MOSFETs using the NEGF framework in both the RS and MS representations. Different transport models in RS will be implemented in the FETMOSS simulator, whose original version was capable of simulating NEGF using the MS only. In this work, the 2D NEGF simulation methodology and results of DG MOSFETs using the RS approach are presented. A new computationally efficient method for the NEGF in RS is proposed and implemented in the FETMOSS simulator. Other methods existing in the literature were studied, implemented and compared with the proposed one. Moreover, the MS validity is examined by comparison with the RS. Then, a fast method to check the MS validity is investigated and implemented and some ideas are suggested to further reduce the MS computational burden. Exploiting the 2D quantum simulator built in this work, design and simulation of nanoscale DG MOSFET at the end of the International Technology Roadmap for Semiconductors (ITRS) is carried. The design of a 10 nm DG MOSFET is presented using a design optimization procedure. For the first time to be reported, the obtained results show that satisfying both the on-current and the switching speed requirements of the ITRS for the high power device is possible. The use of high-k material gate dielectric for the low standby power 16 nm DG MOSFET is also investigated. Finally, the effect of 2D electrostatics on the nanoscale DG MOSFETs capacitance is studied.

Key Words: DG MOSFETs, Mode-space, NEGF, Quantum transport, Realspace.

ACKNOWLEDGMENT ‫الحمد هلل رب العالمين‬

I would like first to thank my supervisors Prof. Dr. Magdy Mahmoud Ibrahim, Prof. Dr. Hani Fikry Ragai and Prof. Dr. Wael Fikry Farouk for their continuous guidance, encouragement, help and patience. I learned so many valuable things from them. I would be ungrateful if I fail to express my sincere gratitude and appreciation to Prof. Dr. Wael Fikry Farouk for all what he has given to me from his noteworthy experience and knowledge. His patience with me has been unparalleled. He exerted great efforts in supervising my thesis in the very details. I learned a lot from his smart and logical way in thinking. I am sure that the experience I have gained from working under his supervision will be of great help to me in my future professional and social life. I would like also to thank Dr. Tarik Abdolkader for his early support in the subject. I would like also to thank him for his continuous encouragement and his efforts in revising my thesis. Finally I would like to thank my parents. Their patience, care, love and encouragement are what guided me through my whole life. I pray to God that I will always be a good faithful son to them.

CONTENTS LIST OF FIGURES………………………………………………………………..viii LIST OF TABLES…………………………………………………………………xvi LIST OF SYMBOLS……………………………………………………………..xviii

INTRODUCTION………………………………………………………………..….1

DOUBLE-GATE MOSFETS……………………………………………………......2 1.1 BULK MOSFET SCALING CHALLENGES…………………………........3 1.2 DOUBLE-GATE MOSFETS…………………………………………..........6 1.3 QUANTUM EFFECTS IN NANOSCALE DG MOSFETs………....…......10 1.3.1 Quantum Confinement………………………………………..…….....10 1.3.2 Quantum Mechanical Tunneling………………………………...….....14 1.3.3 Quantum Interference…………………………………………….........20 1.4 OTHER MULTI-GATE MOSFET STRUCTURES……………….……....21 QUANTUM TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM……………………………………………………..…24 2.1 SELF-CONSISTENT FIELD METHOD…………………………………..25 2.2 COHERENT TRANSPORT THROUGH A SINGLE ENERGY LEVEL. ..26 2.3 COHERENT TRANSPORT THROUGH A 1D DEVICE…………………29 2.3.1 Discrete Leads E (k ) Relation…………………………………...…....34 2.4 FROM WAVEFUNCTION TO GREEN’S FUNCTION…………...…….......36 2.4.1 The Green’s Function…………………………………………….........37 2.4.2 The Spectral Function…………………………………………………37 2.4.3 The Correlation Function……………………………………………...38 2.4.4 The Scattering Function…………………………………………….....38 2.4.5 Transmission Function…………………………………………….......39 2.5 COMPUTATIONALLY EFFICIENT NEGF-BASED METHODS………39 2.5.1 Recursive Green’s Function Algorithm……………………………….40 2.5.2 Contact Block Reduction Method…………………………………......45 DOUBLE-GATE MOSFETS SIMULATION USING THE NEGF IN REALSPACE……………………………………………………………………….………50 3.1 QUANTUM TRANSPORT MODEL……………………………….…..….....51 3.2 REAL-SPACE APPLIED TO DG MOSFETS……………………….……….53 3.2.1 The Real-Space Hamiltonian Matrix………………………………….......53 3.2.2 The NEGF in Real-Space……………………………………………...….55 3.3 COMPUTATIONALLY EFFICIENT METHODS………………………......57 3.3.1 The Recursive Green’s Function Algorithm…………………………......58 3.3.2 The Contact Block Reduction Method…………………………………....59 3.3.3 Proposed Method Based on Gauss Elimination……………………….....62

3.4 INTEGRATION INTO FETMOSS…………………………………………...65 3.4.1 Initial Guess for Real-Space ………………………………………………66 3.4.2 Numerical Integration …………………………………………………….66 3.4.3 Silicon Anisotropy ………………………………………………………..66 3.5 BENCHMARKING FETMOSS………………………………………………68 3.6 NUMBER OF EIGENSTATES IN THE CBR METHOD……………………72 3.7 SIMULATION TIME COMPARISON……………………………………….84 DOUBLE-GATE MOSFETS SIMULATION USING THE NEGF IN MODESPACE……………………………………………………………………………….88 4.1 MODE-SPACE APPROACH……………………………………………...88 4.1.1 Hamiltonian Matrix in Uncoupled Mode-Space………………….…...91 4.1.2 The NEGF in Uncoupled Mode-Space ………………………………..92 4.1.3 UNCOPULED MODE-SPACE VALIDATION BY COMPARISON WITH REAL-SPACE……………………………………………………………………..93 4.2 COMPUTATIONALLY EFFICIENT METHOD TO CHECK THE UNCOUPLED-MODE SPACE VALIDITY…………………………………….102 4.3 FAST MODE-SPACE SIMULATION…………………………………...109 DESIGN & SIMULATION STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END……………………………………………...…………………...113 5.1 DESIGN OF A HIGH-POWER 10 NM FINFET……...………………....117 5.2 STUDY OF HIGH-K MATERIAL FOR LOW-STANDBY-POWER 16.NM FINFET………………………………………………………...…………..…….123 5.3 COMPARISON OF 1D AND 2D DG MOSFET CAPACITANCE..…..126 CONCLUSIONS AND FUTURE WORK……………………………………….129

REFERENCES…………………………………………………………………….131

LIST OF FIGURES Figure 1.1: Simulated constant potential contours of (a) a long-channel and (b) a short-channel nMOSFET. The contours are labeled by their band bending w.r.t neutral p-type regions. The drain is biased at 3.0 V and the gate voltage is slightly below the threshold [Tau98]. .......................................................................................................... 4 Figure 1.2: Short-channel threshold roll-off: Measured low- and high-drain threshold voltages of n-MOSFET versus channel length [Tau98] ................................................ 4 Figure 1.3: The potential barrier along a MOSFET channel at zero and high drain voltage for (a) long-channel case and (b) short-channel case ........................................ 5 Figure 1.4: Subthreshold characteristics of long- and short-channel devices at low and high drain bias [Tau98]. ................................................................................................. 5 Figure 1.5: Simulated IV curves for11-nm-channel-length bulk MOSFETs. Each gray curve corresponds to a different random placement of the dopants in keeping with the designed average doping profiles. The solid black curve is the geometric average of the curves, while the dashed curve is the expected IV curve based on continuum doping profiles [Fra02]. ................................................................................................. 6 Figure 1.6: Structure of (a) double-gate, (b) ground-plane (or back-gate FET, BG). The upper figures are cross-sections of the devices; the lower figures illustrate their energy band diagrams. The gate work functions of the top and bottom gates can be the same for symmetric DG FET (SDG) or different for asymmetric DG FET (ADG) [Won02]. ........................................................................................................................ 7 Figure 1.7: Electric field lines originating from the drain for (a) bulk MOSFET and (b) DG MOSFET. In bulk MOSFETS, electric field lines can penetrate to the source that causes the DIBL. IN DG MOSFET, in contrast, lines are terminated on either top or bottom gate ................................................................................................................ 7 Figure 1.8: Simulated (using FIELDAY) threshold voltage roll-off for DG MOSFET with tox=1.5nm and 1nm. The silicon channels are undoped. Gate material is assumed to have a mid-gap work function [Won98]. ................................................................... 8 Figure 1.9: Simulated (using FIELDAY) DIBL for double-gate MOSFET with tox=1.5nm and 1nm. The Si channel is undoped and T=300K [Won98]. ..................... 8 Figure 1.10: Subthreshold swing obtained with Medici simulations versus channel length for two different values of silicon thickness (NA=1016 cm-3, tox=1.5 nm, and VDS=0.1 V) [Che03] ...................................................................................................... 9 Figure 1.11: Schematic illustrating band-to-band tunneling in a reverse biased p-n junction .......................................................................................................................... 9

Figure 1.12: Double-gate SOI MOSFET. .................................................................... 10 Figure 1.13: Energy band diagram in the transverse direction of DG-SOI MOSFET of Fig. 1.12. ...................................................................................................................... 11 Figure 1.14: Energy gap widening and the splitting of the conduction band into subbands due to carrier confinement. .......................................................................... 11 Figure 1.15: (a) Simulated electron density dependence on gate voltage. (b) Threshed voltage dependence on silicon layer thickness. Gate oxide layer thickness is 7 nm. Buried oxide (BOX) layer thickness is 80nm. Channel doping concentration (NA) is 1017 cm-3. Temperature is 300 K. [Yas93] ................................................................... 12 Figure 1.16: (a) The carrier density distribution near the interface of a DG device in both classical and quantum-mechanical pictures and (b) The average distance of charge carriers from the interface (charge centroid) [Jan98]. ...................................... 12 Figure 1.17: Comparison of nMOS gate capacitance calculated both classically and quantum-mechanically for oxide thickness = 3.0 nm and substrate doping = 3x1017 cm-3. [Ric01] ............................................................................................................... 13 Figure 1.18: Influence of the semiconductor thickness on the electron concentration distribution at a constant surface potential, according to a quantum model [Bog98] .. 13 Figure 1.19: Electron mobility in the symmetrical double-gate transistor versus the silicon film thickness. .................................................................................................. 14 Figure 1.20: Gate oxide thickness (nm) versus technology node (nm) [Wong05]. .... 15 Figure 1.21: Energy band diagrams illustrating (a) FNT and (b) DT. ......................... 15 Figure 1.22: Measured (dots) and simulated (solid lines) tunneling currents in thinoxide polysilicon-gate MOS devices. The dashed line indicates a tunneling current level of 1A/cm2 [Lo97]. ............................................................................................... 16 Figure 1.23: Graph showing the dielectric constants (and variations) of different high-k gate material candidates as reported in the literature [Ban05]. ........................ 17 Figure 1.24: Electron mobility in a symmetrical DG MOSFET device with different insulator taking into account the effect of remote phonon scattering [Col08]. ........... 18 Figure 1.25: (a) Cross-sectional view of 8-nm channel length EJ-MOSFET. (b) Temperature dependence of the subthreshold swing at various channel lengths [Kaw99]. ...................................................................................................................... 19 Figure 1. 26: Simulated subthreshold characteristics of a double-gate MOSFET with ultrathin (t=2.5nm) intrinsic channel and oxide thickness tox=2.5nm (a) gate length L=10nm. (b) L=2.5nm. Dotted lines show the ideal slope of the curves (60mV/decade) and show also the estimated gate leakage current. [Sve03] ......................................... 19

Figure 1.27: The electron density along the channel of a DG-SOI MOSFET calculated by both the semi-classical BTE (dashed lines) and the quantum transport model (solid lines) [Ren01]............................................................................................................... 20 Figure 1.28: FinFET Device Structure........................................................................ 21 Figure 1.29: Triple-gate (trigate) MOSFET................................................................ 21 Figure 1.30:  -gate device cross section. The thickness and width of the device are t si and W si . The radius of curvature of the top and bottom corners is rtop and rbot respectively. ................................................................................................................. 22 Figure 1.31: dg m / dVG in the  -gate MOSFET of Fig. 1.29. Gate is N+ polysilicon A: rtop  rbot  1nm , B: rtop  rbot  5nm ...................................................................... 23 Figure 2.1: Design sequence to achieve desired customer need ................................ 25 Figure 2. 2: Flow chart illustrating the self-consistent field method used in device simulation using the NEGF .......................................................................................... 26 Figure 2. 3: Single level device sandwiched between two contacts. .......................... 27 Figure 2. 4: A 1D device connected to two contacts. All wave functions in the device can be represented as left incident D(1) , right incident  D( 2 ) and bound state  D( B ) . .... 30 Figure 2. 5: A 1D discretized device connected to two semi-infinite contacts leads. 32 Figure 2. 6: Continuous and discrete E (k) relations plotted on the same graph. The discrete relation is in close agreement with the continues for the range E (k )  U 0  t x ...................................................................................................................................... 35 Figure 2. 7: A layered structure connected to 3 contacts. Contact 1 is coupled to the left most layer, contact 2 is coupled to the right most layer and contact 3 is coupled to all the layers and is coupling them together. The RGF can’t be used as along as contact 3 is coupled to the device. ............................................................................... 40 Figure 2. 8: The overall system can be divided into two parts left L and right R ....... 41 Figure 2. 9: (a) The left connected Green’s function of the left most q layers (q=7 for example) and the isolated Green’s function of the q+1th layer are given. (b) We seek the left connected Green’s function of the overall q+1 layer. ..................................... 42 Figure 2. 10 The forward recursion sequence. Initially, the left connected green’s L1 function of the first layer g1,1 is calculated. Then, the first two layers left connected green’s function block g 2L,22 is obtained and so on until the left connected Green’s function block g NLNx ,xN x of the overall system is obtained. .............................................. 42

Figure 2. 11: The backward recursion sequence. Initially, G N x , N x is used to calculate GN x 1, N x , GN x , N x 1 and GN x , 1N x 1 . The same procedure is repeated until the diagonal the upper and lower off diagonal blocks of the fully connected Green’s function are obtained. ....................................................................................................................... 44

Figure 2. 12: Three terminal device transmission function between terminal 1 and 2 examined in [Mam03], versus energy obtained using the CBR method. The solid line corresponds to using all of the isolated system eigenstates, dotted line for only 7% of the eigenstates using von Neumann boundary condition and dashed line for the use of 7% of the eigenstates using Dirichlet boundary condition. ......................................... 48 Figure 3. 1: A model double-gate MOSFET used in this work ................................. 51 Figure 3. 2 (a) Mesh nodes vertically counted using one index. (b) The simulation domain is enclosed by the bolded black line. .............................................................. 54 Figure 3. 3: (a) The active device is composed of layers adjacent to each other. Each layer has on-site energy α and, the coupling energy between any adjacent layers is β . (b) Each layer is composed of smaller sites (nodes) with on-site energy 2tx+2ty and coupling hoping energy –ty between vertically adjacent nodes and –tx between horizontally adjacent nodes.......................................................................................... 55 Figure 3. 4: The overall system composed of the active device and the contacts. ...... 56 Figure 3. 5: Grid points that connect the active device to the contacts are numbered first then the interior device points are numbered. ...................................................... 59 Figure 3. 6: The flow chart of FETMOSS. .................................................................. 66 Figure 3. 7: Constant-energy surfaces characterizing the conduction-band structure near the conduction band edge in Si. ........................................................................... 67 Figure 3. 8: The IDS-VGS C/CS of device 1 at VDS=25mV. The left axis is the log-scale while the right axis is the linear-scale. ......................................................................... 70 Figure 3. 9: The IDS-VGS C/CS of device 1 at VDS=0.7V. The left axis is the log-scale while the right axis is the linear-scale. ......................................................................... 70 Figure 3. 10: The IDS-VGS C/CS of device 2 at VDS=25mV. The left axis is the log-scale while the right axis is the linear-scale. ......................................................................... 71 Figure 3. 11: The IDS-VGS C/CS of device 2 at VDS=0.9V. The left axis is the log-scale while the right axis is the linear-scale. ......................................................................... 71 Figure 3. 12 : The electron density along the x-direction at midpoint along the ydirection for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation........................................................................................................... 73

Figure 3. 13: The transmission function versus energy for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation. ............................... 74 Figure 3. 14: The energy spectrum of the current for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation. ....................................... 75 Figure 3. 15: The electron density along the x-direction at midpoint along the ydirection for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation. ....................................................................................... 77 Figure 3. 16: The transmission function, on linear scale, versus energy for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation..... 77 Figure 3. 17: The energy spectrum of the current for device 1 in the off-state for 1, 4, 10, 20 and 40 % of the eigenstates considered in the simulation. ............................... 78 Figure 3. 18: The transmission function, on log scale, versus energy for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation..... 78 Figure 3. 19: The electron density along the x-direction at midpoint along the ydirection for device 2 in the on-state for 3, 5, 6 and 40 % of the eigenstates considered in the simulation........................................................................................................... 79 Figure 3. 20: The transmission function, on linear scale, versus energy for device 2 in the off-state for 3, 4, 6 and 40 % of the eigenstates considered in the simulation....... 80 Figure 3. 21: The energy spectrum of the current for device 2 in the on-state for 3, 5, 6, and 40 % of the eigenstates considered in the simulation. ...................................... 80 Figure 3. 22: The electron density along the x-direction at midpoint along the ydirection for device 2 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation. ....................................................................................... 82 Figure 3. 23: The transmission function, on linear scale, versus energy for device 2 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation... 82 Figure 3. 24: The transmission function, on log scale, versus energy for device 1 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation... 83 Figure 3. 25: The energy spectrum of the current for device 2 in the off-state for 3, 20, 40, and 50 % of the eigenstates considered in the simulation. .................................... 83 Figure 3. 26: The self-consistent error versus time using full matrix inversion of the Green’s function in the real-space. .............................................................................. 85 Figure 3. 27: The self-consistent error versus time using proposed method for the Green’s function in the real-space. .............................................................................. 85 Figure 3. 28: The self-consistent error versus time using the RGF in the real-space. . 86

Figure 3. 29: The self-consistent error versus time using the CBR method. ............... 86 Figure 4. 1: A model DG MOSFET divided into vertical slices. The horizontal spacing between the slices is Δx .................................................................................. 89 Figure 4. 2: Percentage error in the on-current I on and the off-current I off versus TSi. The percentage difference is calculated as I (%)  1  IUMS / I RS  100 .................... 95 Figure 4. 3: The off-state 1st order mode eigenfunction versus the y -position at different x -positions and at TSi  1.5nm ....................................................................... 96 Figure 4. 4: The on-state 1st order mode eigenfunction versus the y -position at different x -positions and at TSi  1.5nm ........................................................................ 96 Figure 4. 5: The off-state 3rd order mode eigenfunction versus the y -position at different x -positions and at TSi  1.5nm ........................................................................ 97 Figure 4. 6: The on-state 3rd order mode eigenfunction versus the y -position at different x -positions and at TSi  1.5nm ....................................................................... 97 Figure 4. 7: The off-state 1st order mode eigenfunction versus the y -position at different x -positions and at TSi  7nm ......................................................................... 98 Figure 4. 8: The on-state 1st order mode eigenfunction versus the y -position at different x -positions and at TSi  7nm .......................................................................... 99 Figure 4. 9: The off-state 3rd order mode eigenfunction versus the y -position at different x -positions and at TSi  7nm ........................................................................ 99 Figure 4. 10: The on-state 3rd order mode eigenfunction versus the y -position at different x -positions and at TSi  7nm ....................................................................... 100 Figure 4. 11: Percentage of the mth subband current I m to the total current I versus

Tsi for the off- and on-states. ..................................................................................... 101 Figure 4. 12: 3rd order mode subband energy, for example, versus the x - position at TSi  7nm and for the off- and on-states. ................................................................... 101 Figure 4. 13: R11 as a function of El and x at TSi  7nm in the off-state. The scale on the right side of the graph gives the value of the coupling term versus the intensity where brighter regions mean larger value of the coupling term. Moreover, the lowest 1 order mode subband energy E is indicated by the solid line versus x. .................. 103 Figure 4. 14: The 2D LDOS for the off-state 1st order mode versus El and x -position at Tsi  7nm . E 1 is indicated by the solid line versus x . ........................................ 104

Figure 4.15: R12 as a function of El and x at Tsi  7nm in the off-state. E 1 and E 2  are indicated by the solid line versus x . .................................................................... 105 Figure 4.16: R11 as a function of El and x at Tsi  7nm in the on-state. E 1 is also indicated by the solid line versus x . ........................................................................... 105 Figure 4.17: R12 as a function of El and x at Tsi  7nm in the on-state. E 1 and E 2  are indicated by the solid line versus x . .................................................................... 106 Figure 4.18:

The maximum absolute value of the total coupling effect on the

m th mode Rm where m  1, 2, 3, 4 and 5 versus TSi in the off-state. ........................ 107 Figure 4.19:

The maximum absolute value of the total coupling effect on the

th

m mode Rm where m  1, 2, 3, 4 and 5 versus TSi in the on-state. ........................ 107 Figure 4. 20: Total coupling effect of all the modes weighted by their percentage contribution (as defined in (4-24)) versus TSi in the off and on-state. ...................... 108 Figure 4.21: The 1st order mode eigenfunction versus the vertical position at TSi  1.5nm . Both simulation result (SR) and analytical approximation (AA) are shown in the plot. SR obtained at VGS  0 and VGS  VDD are the identical and only one of them is shown. ................................................................................................ 110 Figure 4.22: The 3rd order mode eigenfunction versus the vertical position at TSi  1.5nm . Both simulation result (SR) and analytical approximation (AA) are shown in the plot. SR obtained at VGS  0 and VGS  VDD are the identical and only one of them is shown ................................................................................................. 110 Figure 5. 1: (a) 3-D view of a FinFET, (b) A-A′ cross- section and, (c) B-B′ crosssection. ....................................................................................................................... 115 Figure 5. 2: The metal gate work function needed to satisfy ITRS off-current (< 0.76A / m ) versus the underlap length at different extension lateral abruptness for the design in table 5.3. ............................................................................................... 118 Figure 5. 3: The on-current versus the underlap length at different extension lateral abruptness for the design in table 5.3. The work function was set to satisfy the ITRS off-current target (< 0.76A / m ). ................................................................................ 119 Figure 5. 4.: The metal gate work function needed to satisfy ITRS off-current (< 0.76A / m ) versus the underlap length at different extension lateral abruptness for the design in table 5.6. ............................................................................................... 121 Figure 5. 5: The on-current versus the underlap length at different extension lateral abruptness for the design in table 5.6. The work function was set to satisfy the ITRS off-current target (< 0.76A / m ). ................................................................................ 122

Figure 5.6. Verification of our model results against that in [Sai07]. The device parameters and dimensions are: Lg= 20 nm, Tsi= 5 nm, Tox= 1 nm, NA= 1017 cm-3. The drain and source were connected together. ................................................................ 124 Figure 5. 7. High-k material study for LSTP 16nm technology node in the year 2015. .................................................................................................................................... 125 Figure 5. 8: The gate capacitance per unit width versus VGS at different Lg obtained from FETMOSS ......................................................................................................... 127 Figure 5. 9. The gate capacitance per unit area versus VGS at different Lg obtained from FETMOSS and compared to Shred. .................................................................. 128

LIST OF TABLES Table 1-1: Technology generation nodes starting from 1994 to 2008 ........................... 3 Table 2- 1: Single level device scalar quantities and their corresponding matrices for a multi-level device......................................................................................................... 29 Table 3-1: Transport, width, and confinement effective masses and degeneracies for three different technologically important semiconductor wafer orientations [Rah05]. ...................................................................................................................................... 67 Table 3-2: The simulated devices dimensions, doping concentration, material parameters, simulator options and FD grid. ................................................................. 69 Table 3-3: CBR method results for device 1 in the on-state ........................................ 73 Table 3-4: CBR method results for device 1 in the off-state. ...................................... 76 Table 3-5: CBR method results for device 2 in the on-state ........................................ 79 Table 3- 6: CBR method results for device 2 in the off-state. ..................................... 81 Table 3-7: CBR method results summary. ................................................................... 84 Table 3-8: Simulation time summary. These simulations were carried out on a home PC: Intel® Pentium 4 CPU 2.4GHz, 768 MB RAM. .................................................. 87 Table 4-1: List of the device parameters used in the simulation together with the model device shown in Figure 3.1. .............................................................................. 94 Table 4-2: The maximum absolute value of the coupling terms Rnm at TSi  7nm , VGS  0 and VDS  VDD . The dominant terms are highlighted. ................................... 111 Table 4-3: The maximum absolute value of the coupling terms Rnm at TSi  7nm , VGS  VDD and VDS  VDD . The dominant terms are highlighted. ................................. 111 Table 4-4: The maximum absolute value of the coupling terms Rnm at TSi  5nm , VGS  0 and VDS  VDD . The dominant terms are highlighted. ..................................... 112 Table 4-5: The maximum absolute value of the coupling terms Rnm at TSi  5nm , VGS  VDD and VDS  VDD . The dominant terms are highlighted................................... 112 Table 5-1: ITRS the technology requirements for HP, LSTP and LOP DG MOSFETs. .................................................................................................................................... 113

Table 5-2: Design parameters and simulation results reported in [Has08]. ............... 116 Table 5-3: Design parameters for our study .............................................................. 117 Table 5-4: Design parameters for optimal on-current and the resulting switching speed. The rest of parameters are given in table 5.3. ................................................. 120 Table 5-5: Design parameters for near mid-gap work function and the resulting oncurrent and switching speed. The rest of parameters are given in table 5.3. ............. 120 Table 5-6: Design parameters for better manufacturability design. ......................... 121 Table 5-7: Design parameters for optimal on-current and the resulting switching speed using 6 nm fin thickness. The rest of parameters are given in table 5.6. ................... 122 Table 5-8: Design parameters for near mid-gap work function and the resulting oncurrent and switching speed. The rest of parameters are given in table 5.6. ............. 123 Table 5-9: Relative dielectric material constant and conduction band offset for dielectric materials used in our simulation taken from [Nag05]................................ 125 Table 5-10: Simulation parameters for 2D capacitance study .................................. 126

LIST OF SYMBOLS A

Spectral function

AD

Spectral function due to drain contact

AS

Spectral function due to source contact

α

On-site (layer) energy matrix

β

Coupling energy matrix between adjacent sites

Cgg

Total gate capacitance calculated using 2D quantum simulation

Cg,total

Total gate capacitance as defined in the ITRS



Eigenfunction in the y-direction



Density of states of a single level device

y

Spatial grid points spacing in y-direction

x

Spatial grid points spacing in x-direction

E

Total electron energy

EC

Conduction band edge energy

EF

Fermi energy

El

Longitudinal electron energy due to motion in x- and ydirection

Ekz

Transverse electron energy due to motion in z-direction

EOT

Equivalent oxide thickness

(n)

E

Nth mode Subband energy

εI

Dielectric constant of insulator

εSi

Dielectric constant of silicon

f

Fermi-Dirac distribution

F

Modified Fermi function taking transverse modes into account

g

Extension lateral abruptness

gD

Surface Green’s function of the drain contact

gS

Surface Green’s function of the source contact

gL

Left connected Green’s function

G

Retarded Green’s function

Gn

Electron correlation function

Go

Retarded Green’s function of the isolated device

Gp

Hole correlation function

(m)

G

Mth Green’s function

γ

Coupling strength between a contact and a single level device

Γ

Broadening function

ΓD

Broadening function due to drain contact

ΓS ΓS

Broadening function due to source contact (m)

Broadening function of the mth mode due to source contact

ΓD (m)

Broadening function of the mth mode due to drain contact

H

Hamiltonian matrix

Hm

Hamiltonian of the mth mode

ħ

Dirac constant

ISD,leak

MOSFETs subthreshold off-state leakage current

ION

MOSFETs drive current in the on-state

I (m)

Electrical current of the mth mode from source to drain

J

Electrical current density

Jg,limit

MOSFETs maximum gate leakage current density

KB

Boltzmann constant

L

MOSFETs Channel length

LD

MOSFETs drain length

Lg

MOSFETs Physical gate length

LS

MOSFETs source length

LU

MOSFETs underlap region length

Λ

Scale length for 2-D effects in MOSFETs

ml

Longitudinal effective mass in Si

mt

Transverse effective mass in Si

m*

Electron effective mass

mx*

Electron effective mass in x-direction

my*

Electron effective mass in y-direction

mz*

Electron effective mass in z-direction

n

Electron concentration

n(m)

Electron density of the mth mode

N

Total number of spatial grid points

NA

Acceptor concentration

ND

Donor concentration

NE

Number of energy grid points

Neigen

Number of eigenstates used in the CBR method

Niter

Number of self-consistent iterations

Nm

Number of modes used in the mode-space representation

Nmode

Number of propagating modes used in the CBR method

NS

Supply function

Nx

Number of spatial grid points in the x-direction

Ny

Number of spatial grid points in the y-direction

p

Holes concentration

φm

Metal work function

(n)

Nth mode expansion coefficient

ψ

2D wavefunction

Ψ

3D wavefunction

Rm

Total coupling effect on the mth mode

Rmm

Coupling effect on the mth mode due to the mth mode itself

Rnm

Coupling effect on the mth mode due to the nth mode

ρ

Density matrix



Contact Self energy

∑D

Drain contact self energy

∑S

Source contact self energy

in



In-scattering self energy

∑out

Out-scattering self energy

∑D(m)

Drain contact self energy of the mth mode

∑S(m)

Source contact self energy of the mth mode

t

1-D coupling energy between adjacent sites in the spatial grid

tbias

Average simulation time per bias point

titeration

Average simulation time per self-consistent iteration

ttotal

Total simulation time

T

Temperature in Kelvin

T(E)

Transmission function

TC

Transmission coefficient

TI

Insulator thickness

TOX

Oxide thickness

TSD (m)

Mth mode transmission function from source to drain

TSi

Silicon film thickness



MOSFET intrinsic delay

Wdm

Maximum depletion depth

INTRODUCTION Rapid device scaling pushes the dimensions of the field effect transistors to the nanometer regime. In this regime of operation quantum effects play an important role in determining the transistor characteristics. These effects can be accurately predicted only using quantum mechanical based device simulation. Non-equilibrium Green’s function formalism (NEGF) provides a rigorous description of quantum transport in nanoscale devices. Computational efficiency is needed to make the device simulation suitable for device design and characteristic prediction. The NEGF method, unfortunately, has the disadvantage of being heavy in computations. The Real-space (RS) representation is the most accurate yet complex representation used in the NEGF. The geometry of fully-depleted double gate (DG) MOSFETs, however, permits the use of a simpler representation, the mode-space (MS) which is computationally efficient. An existing software tool, named as FETMOSS, for quasi two-dimensional (2D) simulation of DG MOSFETs was built before by Dr. Tarik Abdolkader during his PhD thesis. FETMOSS is working under MATLAB environment and is based on the numerical solution of Poisson and Schrödinger equations self-consistently. In this work, we build on FETMOSS and increase its capabilities to carry out true 2D simulation of DG MOSFETs. In chapter 1, bulk MOSFET scaling challenges are highlighted and the DG MOSFET structure is presented. The principal quantum mechanical effects encountered in DG MOSFETs are investigated showing there critical role in determining the device electrical characteristics. Finally, a brief overview on other Multiple-Gate MOSFETs is given. In chapter 2, the important properties of coherent transport and the NEGF for nanoscale device simulation are reviewed. Since the NEGF is computationally intensive, two methods commonly used in the literature that make the NEGF less computationally demanding are also reviewed. These are the Recursive Green’s Function (RGF) and, the more recent one, the Contact Block Reduction (CBR) method. In chapter 3, the 2D simulation methodology and results of DG MOSFETs using the RS approach are presented. The quantum transport model used for the DG MOSFET including our assumptions is given. Then, the application of the RGF and CBR methods to the DG MOSFET is presented. Moreover, a new computationally efficient method is proposed. Finally, all these methods were implemented and integrated into the FETMOSS simulator and compared. In chapter 4, a review of MS approach is given. Then, the uncoupled mode-space (UMS) validity is examined by comparing against the RS. A fast method to check the UMS validity is then investigated and implemented. Finally, some ideas are suggested to reduce the MS computational burden. In chapter 5, design of a 10 nm DG MOSFET is presented using the developed 2D simulator. Design optimization procedure is applied, that results in satisfying both the oncurrent and the switching speed requirements of the ITRS for the high power device. Then, the use of high-k material gate dielectric for the low standby power 16 nm gate length DG MOSFET is investigated. Finally, the effect of 2D electrostatics on the device capacitance is investigated. Finally, the thesis ends by extracting conclusions and stating future work that might be done based on this work.

CHAPTER 1

DOUBLE-GATE MOSFETS

Field Effect Transistor (FET) invention was a breakthrough in the field of integrated circuits. The basic idea of the FET was patented by Lilienfeld in 1930 [Lil30]. The idea was extensively explored and experimentally studied by numerous researchers at different universities and companies [Loj07]. In 1960, the Metal Oxide Semiconductor Filed Effect Transistor (MOSFET) was practically demonstrated by Dawn Kahng and Mohamed Attalla [Kah60]. The demonstrated transistor was pchannel MOSFET with 25 µm drawn channel length and 250 µm channel width. Technology generation is denoted by a figure, for example 0.5 µm. The figure is called a node and refers to the minimum metal line width that can be afforded by the technology. This restricts the minimum drawn MOSFET channel length that can be fabricated using the technology. Each new technology generation node is about 0.7 smaller than the previous node. This periodic size reduction is called scaling. Smaller MOSFETs are desirable for several reasons: 1) Increasing the packing density and, as a result, the number of functions per chip. This, in turn, reduces the cost per function, 2) Smaller transistors switch faster because smaller gate capacitance leads to smaller RC delay and 3) Shorter transistors have higher transconductance that provide a higher gain-bandwidth product. Transistor scaling was very aggressive in the past 40 years. In 1971, the technology generation node was the 10 μm. At that time, Intel® introduced its first processor (the “4004” microprocessor). The “4004” employed a 10 μm silicon-gate PMOS technology and was built of 2,300 transistors [Wik01]. The current technology generation node is the 45 nm. Recently, Intel® announced that production for the “Penryn” processor will be started in 2008. The “Penryn” employs a 45 nm high-k metal-gate CMOS technology with about 800 million transistors [Boh07]. In 1965 Gordon Moore empirically predicted that the number of transistors per chip will be doubled every two years (every technology generation) [Moo65]. In Table1.1 the technology generations are given staring from 1994 to 2008. It is interesting to note that Moore’s law is satisfied across the entire technology nodes. It worth mentioning that, this agreement between technology nodes and Moore’s law isn’t by chance. The semiconductor industry is very keen to continue its exponential advance.

Table 1-1: Technology generation nodes starting from 1994 to 2008 Year Node ( nm)

1994 500

1997 350

1999 250

2000 180

2002 130

2004 90

2006 65

2008 45

It is predicted that conventional CMOS technology will come to its end at the 22 nm node [Tho05]; scaling bulk MOSFET beyond the 45 nm technology is accompanied with great challenges. Challenges include Short-Channel Effects (SCEs), Drain Induced Barrier Lowering (DIBL), random dopant fluctuation and high leakage current [Sko05]. Industrial and academic efforts are now being exerted to find a replacement for conventional bulk MOSFET. There are two main directions: 1) adopting new transistor structures or 2) adopting new materials other than Silicon (Si) [Sko05]. Searching for new device structure that circumvents the forgoing challenges while, in the same time, scalable to 22 nm, Double-Gate (DG) MOSFET is found to be the most promising candidate [Bans05, Kim05, Lee05, Wu05, Xio05, Yin05]. This is because, for a given gate length and oxide thickness, DG MOSFETs show better short-channel behavior, lower leakage current, higher drive current, and smaller subthreshold swing than other MOSFET structures [Xio05]. In the same time, they have better power consumption [Kim05] and scaling capabilities [Yin05]. In this chapter, an overview of the bulk MOSFET scaling challenges is given in section 1.1. In section 2.1 the DG MOSFET structure is presented and then, the advantages of the DG MOSFET over bulk MOSFET are explained. The principal Quantum Mechanical Effects encountered in DG MOSFETs are illustrated in section 1.3 showing how these effects play a critical role in determining the device electrical characteristics. Finally, a brief overview on other Multiple-Gate (MG) MOSFETs is given in section 1.4

1.1

BULK MOSFET SCALING CHALLENGES

Scaling the MOSFET beyond the 45 nm technology is accompanied with great challenges. For digital applications, scaling challenges include controlling leakage currents and SCEs, increasing drain saturation current while reducing the power supply voltage, and maintaining control of device threshold voltage across the chip and from chip to chip. For analog/mixed-signal/RF applications, the challenges additionally include sustaining linearity, low noise figure, high power-added efficiency, and good transistor matching [Sko05] The origin of SCEs can be understood by analyzing the potential distribution inside the MOSFET. For long-channel MOSFETs, the effects of the source- and drain-body junctions on the electrostatic potential in the depletion region below the channel are insignificant. On the other hand, for the short-channel MOSFETs, the effects of the source- and drain-body junctions are very important and results in a two-dimensional (2-D) electric field pattern in the depletion region as shown in Fig. 1.1. The constant potential contours in Fig. 1.1 (a) are almost along the y-direction indicating that a 1D vertical field is dominating in the long MOSFET. On the other hand, the constant potential contours in Fig. 1.1 (b) are curved indicating the 2D nature of the electric field. The 2D field pattern arises from the proximity of the drain and source regions [Tau98].

4

CHAPTER 1: DOUBLE-GATE MOSFETS

Figure 1.1: Simulated constant potential contours of (a) a long-channel and (b) a short-channel nMOSFET. The contours are labeled by their band bending w.r.t neutral p-type regions. The drain is biased at 3.0 V and the gate voltage is slightly below the threshold [Tau98].

An important short-channel phenomenon is the dependence of the threshold voltage on the transistor channel length [Tau98]. The classical explanation of this phenomenon is based on the following: the depletion charges associated with the source- and drain-body junctions reduce the charges required for the gate to obtain strong inversion. The relative importance of this reduction increases as the channel length decreases. Therefore, the threshold voltage decreases as the channel length decrease as shown in Fig. 1.2. In a CMOS VLSI technology, channel length varies statistically from chip to chip, wafer to wafer, and lot to lot due to process tolerances [Tau98]. The SCE is therefore an important consideration in device design; one must ensure that the threshold voltage doesn’t become too low for the minimum-channel-length device on the chip.

Figure 1.2: Short-channel threshold roll-off: Measured low- and high-drain threshold voltages of nMOSFET versus channel length [Tau98]

Another important short-channel phenomenon is the Drain Induced Barrier Lowering (DIBL) [Tau98]. The DIBL can be illustrated considering the potential barrier to electrons for an n-channel MOSFET seen from the source. This barrier prevents electron current from flowing to the drain. When the gate voltage is below the threshold voltage, there are only a limited number of electrons injected from the source over the barrier and collected by the drain (subthreshold current). In the

long-channel case, the potential barrier is flat over the most part of the channel. The drain electric field only affects the end of the channel potential near the drain as shown in Fig. 1.3 (a). As the channel length is shortened, the drain field penetrates deeply into the middle of the channel, which lowers the potential barrier near the source as shown in Fig. 1.3 (b). The lowering of the barrier seen at the source due to the drain voltage is referred to as DIBL. The DIBL causes the threshold voltage decrease as the drain voltage is increased as shown in Fig. 1.2. Thus, the DIBL causes a substantial increase of the subthreshold current when the drain voltage is increased [Tau98]. ΔEB

(a)

VDS=0

VDS=0

VDS=VDD

VDS=VDD (b)

Figure 1.3: The potential barrier along a MOSFET channel at zero and high drain voltage for (a) longchannel case and (b) short-channel case

Fig. 1.4 shows the subthreshold characteristics of long and short channel devices at different drain bias voltages. For long devices, the subthreshold current is independent of drain voltage. For short-channel devices, however, there is a parallel shift of the curve to a lower threshold voltage for high drain bias conditions. At even shorter channel lengths, the subthreshold slope starts to degrade as the barrier is controlled by the drain and not only by the gate.

Figure 1.4: Subthreshold characteristics of long- and short-channel devices at low and high drain bias [Tau98].

To find the minimum gate length at each generation of technology, it is needed to analyze the 2-D potential distribution in modern FETs. The most common analysis way is to use a 2-D simulation tools. However, it was shown that the full 2-D potential distribution for a MOSFET can be solved for in a slightly simplified geometry [Fra98]. This analysis yields a scale length for the 2-D effects in MOSFET’s. To the first order, the scale length is given by:

6

CHAPTER 1: DOUBLE-GATE MOSFETS

  Wdm  TI  Si /  I

(1-1)

where  is the scale length, Wdm is the maximum channel depletion depth, TI is the insulator thickness, and  Si /  I is the ratio of dielectric constants of silicon and the insulator. The scaling limits for Si MOSFETs have been extensively studied for different applications [Fra01]. The channel length limit, in order to have good electrostatic control of the gate on the transistor channel, was estimated to be a multiple of the scaling length; typically from  to 2 . The scaling length can be reduced by reducing the depletion width. Reducing the depletion width in bulk MOSFETs is carried out by increasing the channel doping. Increasing the channel doping faces a serious problem in nano-scale transistors, namely, randomness in the exact location of dopant atoms. Average concentration of doping is defined by ion implantation and diffusion process. However, the stochastic nature of ion implantation and diffusion process leads to fluctuations in the local doping concentration [Key75] .This in turn causes device-todevice variation in MOSFET threshold voltages [Drag98]. Moreover, the presence of random dopants in the source and drain regions causes the boundary between source/channel and drain/channel to be randomly varying depending upon the donors’ locations [Fer06]. Fig. 1.5 depicts an example of the statistical variation expected to occur in an 11-nm bulk MOSFET due to random dopant placement [Fra02]. It is clear that the threshold voltage variations are too much to be accepted.

Figure 1.5: Simulated IV curves for11-nm-channel-length bulk MOSFETs. Each gray curve corresponds to a different random placement of the dopants in keeping with the designed average doping profiles. The solid black curve is the geometric average of the curves, while the dashed curve is the expected IV curve based on continuum doping profiles [Fra02].

1.2

DOUBLE-GATE MOSFETS

Due to the great challenges facing the scalability of bulk MOSFETs, industrial and academic efforts are now being exerted to find a replacement. There are two main directions, adopting new transistor structures or adopting new materials [Sko05]. New transistor structures seek to improve the electrostatics of the MOSFET (i.e. to increase the gate voltage control on the transistor) to reduce the SCEs. Among these structures is the Double-Gate (DG) MOSFET. The DG FET shown in Fig. 1.6 (a) was proposed in 1984 [Sek84]. It was shown that one can improve the immunity to SCEs using a top as well as a bottom gate. The usage of two gates, instead of one, showed excellent

gate control on the channel potential. Moreover, the usage of two gates increase the device current capability due to larger gate capacitance per unit area.

Figure 1.6: Structure of (a) double-gate, (b) ground-plane (or back-gate FET, BG). The upper figures are cross-sections of the devices; the lower figures illustrate their energy band diagrams. The gate work functions of the top and bottom gates can be the same for symmetric DG FET (SDG) or different for asymmetric DG FET (ADG) [Won02]. Gate

Gate

Source

Drain

Source

Drain Gate

(a) (b) Figure 1.7: Electric field lines originating from the drain for (a) bulk MOSFET and (b) DG MOSFET. In bulk MOSFETS, electric field lines can penetrate to the source that causes the DIBL. IN DG MOSFET, in contrast, lines are terminated on either top or bottom gate

The reason of the excellent gate control on the channel potential can be explained as follows (see Fig. 1.7). Electric field lines originating from the drain can penetrate to points near the source or even to the source itself in bulk MOSFETs as shown in Fig. 1.7 (a). In contrast, when another gate is added at the bottom, as shown in Fig. 1.7 (b), the electric field lines will properly be terminated at either the top or bottom gate. Thus the drain control on the channel potential is reduced leading to more immunity to DIBL. When the silicon channel is thin enough, two potential advantages for DG MOSFETs over the bulk MOSFETs appear [Won02]: 1. The channel potential becomes tightly coupled to the gate potential; the 2D short-channel effects are reduced. This leads to shorter allowable channel lengths compared to bulk MOSFET without need to increase the channel

8

CHAPTER 1: DOUBLE-GATE MOSFETS

doping. We can say that the SCEs in DG MOSFET are controlled by device geometry compared to doping in bulk MOSFETs. Fig. 1.8 depicts the threshold voltage roll-off for DG MOSFET versus channel length at different channel thickness. The roll-off for the thinnest channel occurs at small channel length compared to other thicknesses. The same is observed for the DIBL from Fig. 1.9. Therefore, the transistor with the thinnest channel is the most immune to SCEs.

Figure 1.8: Simulated (using FIELDAY) threshold voltage roll-off for DG MOSFET with tox=1.5nm and 1nm. The silicon channels are undoped. Gate material is assumed to have a mid-gap work function [Won98].

2. The channel becomes fully depleted eliminating depletion capacitance that deteriorates the subthreshold characteristics. Fig. 1.10 depicts the subthreshold swing versus channel length for tsi=10nm and tsi=20nm [Che03]. As the channel becomes thinner, the subthreshold characteristics become steeper and the subthreshold swing becomes the smaller. This gives the DG MOSFET the advantage of larger gate overdrive for the same power supply and the same off current. Moreover, as a result of lightly-doped/undoped channel, the DG MOSFETs has two more advantages over bulk MOSFET, namely, better carrier transport and reduced scaling limitation due to drain-to-body band-to-band tunneling leakage.

Figure 1.9: Simulated (using FIELDAY) DIBL for double-gate MOSFET with tox=1.5nm and 1nm. The Si channel is undoped and T=300K [Won98].

Figure 1.10: Subthreshold swing obtained with Medici simulations versus channel length for two different values of silicon thickness (NA=1016 cm-3, tox=1.5 nm, and VDS=0.1 V) [Che03]

The better carrier transport property is due to two reasons: reduced Coulomb scattering due to fewer ionized dopants in the lightly-doped/undoped channel, and reduced surface roughness scattering due to a lower surface electric field [Won02]. In bulk MOSFETs, channel doping is employed to set the threshold voltage, and halo or pocket doping is employed to control the short-channel effects. The ionized depletion charges contribute appreciably to the surface electric field. In a DG MOSFET with an undoped channel, there is no ionized depletion charge and, as a result, the surface electric field is contributed only from inversion carriers. When the electric field across a reverse-biased p-n junction approaches 106 V/cm (see Fig. 1.11), significant current flow can occur due to tunneling of electrons from the valence band of the p-region into the conduction band of the n-region [Tau98]. Therefore, band-to-band tunneling (BTBT) leakage current can occur between the body and drain of a MOSFET. BTBT is estimated to be about 1A/cm 2 for body doping of 5x1018cm-3 VDS=1V [Tau98]. The BTBT leakage is proportional to the channel doping and, thus, undoped channel DG MOSFET can reduce this leakage too much.

qVapp EC EV p

n+

Figure 1.11: Schematic illustrating band-to-band tunneling in a reverse biased p-n junction

10

CHAPTER 1: DOUBLE-GATE MOSFETS

1.3

QUANTUM EFFECTS IN NANOSCALE DG MOSFETs

Despite the interesting advantages of DG MOSFET devices, some challenges arise through its design and characterization. The transport of charge carriers in a silicon film of a thickness in the order of few nanometers , sandwiched between two oxide layers each of thickness ~ 1-3 nm, results in significant Quantum Mechanical Effects (QMEs). The QMEs represent a great challenge for the DG MOSFET modeling and simulation. These QMEs can be grouped into three categories: 1. Effects arising due to carrier confinement along the transverse direction (the direction normal to the Si-SiO2 interfaces) 2. Quantum-mechanical tunneling Effects 3. Quantum Interference The details of the preceding effects will be presented in the following three subsections.

1.3.1 Quantum Confinement The cross section of a DG-SOI device is shown in Fig. 1.12 with the energy band diagram along the confinement direction (y-direction) given in Fig. 1.13. There are two confining potential wells that can be identified from the energy band diagram shown in Fig. 1.13: 1) structural confinement and 2) electrical confinement. Structural confinement potential well arises from the conduction band offset between the Si and the gate insulator region (the conduction band offset between Si and SiO2 is nearly 3.2 eV [Han89]). On the other hand, the electrical confinement potential well arises from the steep band bending at the Si surface which is encountered only at large gate voltages or equivalently large electric field in the confinement direction. Indeed, electrical confinement of carriers dominates in higher doping concentration [Col07]. One can also predict that electrical confinement dominates when the silicon thickness is large such that the quantization energy due to structural confinement is small.

Top Gate Tox

Source

Oxide

Drain x

n++

Si

Tsi Tox

n++

Oxide

Bottom Gate y

Figure 1.12: Double-gate SOI MOSFET.

Tox

Tox

SiO2

SiO2

Si Well (1)

Well (2)

E FN

Ei

VG

 ( y)

E FN

VG

y TSi y 0 Figure 1.13: Energy band diagram in the transverse direction of DG-SOI MOSFET of Fig. 1.12.

Due to carrier confinement, the energy levels available for motion in the ydirection are quantized with a continuum for motion in the x-z plane [Strn67]. The total energy of carriers can be written as, E   n  (2k x2 / 2mx  2k z2 / 2mz )

n  1,2,3,...

(1-2)

where kx, and kz are the wave numbers and m*x and m*z are the effective masses in x and z directions, respectively. The first term εn are discrete values of energy corresponds to the eigenenergies in the y-direction. The eigenenergies in the ydirection are inversely proportional to the square of the effective width of the potential well. The second term is a continuous sum of kinetic energies of carriers for motion in x and z directions. As a result of energy quantization, the conduction band splits into subbands each of which has a minimum energy εn, as illustrated in Fig. 1.14. The quantization of energy is stronger (or equally, the separation of energy levels is higher) for smaller Si-film thicknesses and for larger surface electric fields. Subbands

Conduction band

3 2 1 EG (DG )

EG (bulk )

Figure 1.14: Energy gap widening and the splitting of the conduction band into subbands due to carrier confinement.

Due to the presence of the subbands, we can equivalently say that the quantum confinement raises the conduction band edge EC to the lowest order eigenenergy ε1. This shift has a direct influence on the device threshold voltage because it needs more band bending (potential energy lowering) to create the inversion layer. The threshold voltage calculated both classically and quantum-mechanically is given in Fig. 1.15 [Yas93]. The threshold voltage, calculated using the quantum mechanics, increases

12

CHAPTER 1: DOUBLE-GATE MOSFETS

dramatically for ultra-thin devices. It also noticed that this increase is accompanied with the increased sensitivity of threshold voltage with the silicon film thickness. This poses an additional design challenge as fluctuations in silicon thickness will lead to unpredictable values of threshold voltage [Maj01, Tsu05].

(a)

(b)

Figure 1.15: (a) Simulated electron density dependence on gate voltage. (b) Threshed voltage dependence on silicon layer thickness. Gate oxide layer thickness is 7 nm. Buried oxide (BOX) layer thickness is 80nm. Channel doping concentration (NA) is 1017 cm-3. Temperature is 300 K. [Yas93]

The carrier density is peaked at the interface in the classical view. In the view of quantum mechanics, infinite potential barrier at the interface necessitates zero probability density function, and thus, zero carrier density at the interface (see Fig. 1.16 (a)). This quantum effect is called “repulsive boundary condition”. Due to the repulsive boundary condition shown in Fig. 1.16 (a), the average distance of charge carriers from the interface (called charge centroid) is increased approximately by a factor of 2. This can be seen from the quantum-mechanical based simulation results given in Fig. 1.16 (b) [Jan98]. The increased charge centroid means larger effective oxide thickness and smaller gate capacitance. A comparison of gate capacitance calculated both classically and quantum-mechanically is given in Fig. 1.17 [Ric01]. The classical based simulation overestimates the gate capacitance in both inversion and accumulation regimes. Overestimating the gate capacitance leads to over estimating the device driving current capability.

electron cocentration

semi-classical quantum-mechanical

Distance from interface

(a)

(b)

Figure 1.16: (a) The carrier density distribution near the interface of a DG device in both classical and quantum-mechanical pictures and (b) The average distance of charge carriers from the interface (charge centroid) [Jan98].

One more phenomenon occurs due to repulsive boundary condition is the volume inversion. For the DG MOSFET shown in Fig. 1.12, the electron concentration profile along the y-direction with different silicon film thickness has been studied in [Bog98]. The result is shown in Fig. 1.18 where one can observe high electron concentration in the center of the silicon film when thin enough (5nm). This high concentration, corresponding to volume inversion, has a direct impact on the electrons mobility as shown in Fig. 1.19. In thick films there is no interaction between the front and the back channel and there is no volume inversion. If the film gets thinner, volume inversion appears and the mobility is increased because of reduced Si-SiO2 interface scattering [Col07]. In very thin silicon films, however, the carriers are in proximity to the interface and mobility drops with thinner thickness.

Figure 1.17: Comparison of nMOS gate capacitance calculated both classically and quantummechanically for oxide thickness = 3.0 nm and substrate doping = 3x1017 cm-3. [Ric01]

Figure 1.18: Influence of the semiconductor thickness on the electron concentration distribution at a constant surface potential, according to a quantum model [Bog98]

14

CHAPTER 1: DOUBLE-GATE MOSFETS

Mobility

10

20

30

40

Silicon film thickness (nm) Figure 1.19: Electron mobility in the symmetrical double-gate transistor versus the silicon film thickness.

1.3.2 Quantum Mechanical Tunneling One of the most important effects that limit device scaling is the quantummechanical tunneling of carriers through the energy barriers. This tunneling results in leakage current, which increases the static power dissipation and decreases the logic operating margins. Two mechanisms of this leakage current of particular importance [Fran02] for nanoscale DG MOSFETs: tunneling through the gate insulator and source-to-drain tunneling through the channel barrier.

1.3.2.1 Gate insulator tunneling Gate oxide thickness scaling has been the magic tool in controlling short channel effects as MOSFETs gate length have been reduced from one technology node to another. As depicted in Fig. 1.20 the gate oxide thickness is approximately linearly scaled with channel length to maintain the same amount of gate control over the channel to ensure good short channel behavior [Wong05]. With gate-oxides of thickness in the range of few nanometers (< 7 nm), quantum gate-oxide tunneling can affect the electrostatics within the device and generates undesirable gate tunneling leakage current [Ren01]. Tunneling through the gate oxide isn’t pronounced in DG MOSFETs only, but also in conventional bulk MOSFETs. Indeed, the thickness limit for SiO2 is set by gate-to-channel tunneling leakage for a given application [Fra01]. Two tunneling mechanisms can lead to gate current: Fowler Nordheim Tunneling (FNT) [Len69] and Direct Tunneling [Shu94]. The energy band diagrams in Fig. 1.21 illustrate the difference between the two tunneling mechanisms.

Figure 1.20: Gate oxide thickness (nm) versus technology node (nm) [Wong05].

Tox

Tox

Tox

Metal

EFM

EFM

FNT

SiO2

Si

SiO2

Si

FNT Metal

Metal

VG

Metal

Tox

SiO2

SiO2

EFM E

Si

DT

FM

Si

DT

VGV

VG

G

EFS

EEFSFS

(a)

EFS

(b)

Figure 1.21: Energy band diagrams illustrating (a) FNT and (b) DT.

FNT tunneling is effective for Tox > 4 nm and negligible for normal device operation [Tau98]. On the other hand, DT can be very large for thin gate insulator thickness and increases exponentially with reducing the thickness as shown in Figure1.22. The DT is the dominate mechanism for Tox < 2 nm.

16

CHAPTER 1: DOUBLE-GATE MOSFETS

Figure 1.22: Measured (dots) and simulated (solid lines) tunneling currents in thin-oxide polysilicongate MOS devices. The dashed line indicates a tunneling current level of 1A/cm2 [Lo97].

In order to decrease the tunneling leakage current, insulator thickness must increase while equivalent oxide thickness (calculated from the insulator capacitance) must continue to decrease or at least remains the same for good control of the gate on the channel. If a suitable insulator can be found, it would be characterized by three thicknesses [Fran01]: 1) physical insulator thickness: TI , 2) equivalent oxide tunneling thickness: ToxTeq and 3) equivalent oxide capacitive thickness (referred to as EOT): ToxCeq . The goal is to find an insulator with the property that when its ToxTeq is equal

to the minimum SiO2 thickness, its ToxCeq is significantly less than the minimum SiO2 thickness that would enable further scaling. This is possible only with the introduction of CMOS-compatible, high dielectric constant (high-k) materials. This is largely a material issue since its success depends upon achieving high layer uniformity, integration with other Si processes, minimal/controlled reactions with Si and the gate electrode, low fixed-charge, defect, and trap densities in the insulator and at the interface between the insulator and the Si substrate [Fra01]. Fig. 1.23 shows the different dielectric constants of different high-k materials as reported in the literature [Ban05]. They have been deposited under different conditions, different precursors (biochemical substances which are intermediate compounds in a chain of reactions from which a more stable or definitive product is formed) and different thicknesses.

Figure 1.23: Graph showing the dielectric constants (and variations) of different high-k gate material candidates as reported in the literature [Ban05].

As seen in Fig. 1.23, many high-k materials could replace SiO2. Unfortunately, there is a lot of technological challenges appear for the integration of high-k materials with CMOS process [Sha01]: 1) some of these materials are not thermally stable on silicon, forming silicide during CMOS processing, 2) most of them are also not compatible with high-temperature processing and 3) most of these materials are fabricated using techniques such as atomic layer deposition and therefore process integration is also another issue to be considered for use of these materials for ULSI technology. A recent study discussing metal/high-k gate stack issues from both materials science and device physics viewpoints has been reported [07]. There is an emphasis on the critical importance of digging deeper into related material science. The effective electron mobility in inversion layers of CMOS devices employing high-k dielectric was studied and the results are somehow disappointing. Fig. 1.24 depicts the expected electron mobility based on Monte Carlo simulation for a DG MOSFET taking into account the effect of polar phonon scattering [Col08]. It is clear that there is a degradation in the electron mobility when insulators with high dielectric constant are used. It is interesting to note that volume inversion in DG MOSFET doesn’t help to improve the mobility degradation produced by remote polar phonon scattering when high-k dielectrics are used [Col08].

18

CHAPTER 1: DOUBLE-GATE MOSFETS

Mobility (cm2 /Vs)

103

102 1010

1013 1011 1012 Sheet Density ( cm-2 )

1014

Figure 1.24: Electron mobility in a symmetrical DG MOSFET device with different insulator taking into account the effect of remote phonon scattering [Col08].

Other limitations on the use of the high-k gate dielectrics or ultra thin effective gate are the polysilicon gate depletion effect which adds about 0.3–0.7 nm to the total SiO2 gate thickness [Hu96]. Therefore, it is desirable to replace the polysilicon gate with metal gates to minimize the gate depletion effect in addition to reducing the gateline sheet resistance [Che99]. Elimination of the polysilicon depletion reduces the vertical electric field at a given gate voltage, which in turn increases the carrier mobility in the inversion layer [Col08]. 1.3.2.2 Source Barrier Tunneling Direct source-to-drain tunneling through the channel potential barrier is another possible mechanism of tunneling. The contribution of this tunneling current to the MOSFETs subthreshold current was observed for channel lengths shorter than 20 nm at low temperature as shown in Fig. 1.25 [Kaw99]. The cross-sectional diagram in Fig. 1.25 (a) shows the electrically variable shallow junction nMOSFETs which use a second gate over the top of the first gate to induce inversion layers in the source and drain regions. Such inversion layers are much shallower than the usual implanted source–drain extensions, which reduces the influence of the drain on the channel of the FET. As indicated by Fig. 1.25 (b), the subthreshold swing (S-factor) is saturated at low temperatures for the shortest device. This saturation is an indication that the current is dominated by direct tunneling through the channel barrier since such tunneling is not very sensitive to temperature.

Figure 1.25: (a) Cross-sectional view of 8-nm channel length EJ-MOSFET. (b) Temperature dependence of the subthreshold swing at various channel lengths [Kaw99].

Recent theoretical analyses show that source-to-drain tunneling becomes important at room temperature only for channel lengths below ~5 nm [Sve03]. As depicted in Fig. 1.26 source-to-drain tunneling along the channel becomes visible as an upward bend of the curves at negative gate voltages. For L=10 nm, the slope becomes lower than the thermally determined value (dashed line) whereas for L=2.5 nm, the tunneling current dominates the subthreshold current at almost all negative, so that the device works essentially as a tunnel transistor.

VGS (a)

VGS (b)

Figure 1. 26: Simulated subthreshold characteristics of a double-gate MOSFET with ultrathin (t=2.5nm) intrinsic channel and oxide thickness tox=2.5nm (a) gate length L=10nm. (b) L=2.5nm. Dotted lines show the ideal slope of the curves (60mV/decade) and show also the estimated gate leakage current. [Sve03]

20

CHAPTER 1: DOUBLE-GATE MOSFETS

1.3.3 Quantum Interference The phase breaking length is defined as the distance over which the electron wave’s phase is destroyed by some process. It was estimated to be in the range of 100 nm for GaAS and 50 nm for Si [Fer97]. The phase breaking process arises from the interaction of one electron with phonons (lattice vibrations), photons (electromagnetic vibrations) or other electrons. For macroscopic devices, phase randomizing scattering dominates, and electrons wave phase is randomized by collisions and lose their phase during the transport process. Therefore, quantum interference effects can be neglected in macroscopic devices and a semi-classical approach based on Boltzman transport equation BTE can be used to describe the transport [Lun00]. Small devices whose dimensions are much larger than microscopic objects like atoms but not large enough for many scattering events, are called mesoscopic devices [Dat95]. In mesoscopic devices electrons may transport from one side of the device to the other side with little or no scattering at all. In this case, the phase of the electrons wave nature plays an important role in the transport process because electrons can interfere instructively or destructively. Thus, we can’t describe the transport without including the wave nature of the electrons. Figure.1.27 depicts the electron density along the channel of a DG MOSFET calculated by both the semi-classical BTE and the Schrödinger equation. Due to the interference between the incident and reflected electron waves, carrier density oscillations can be seen near the channel barrier edges. At high temperatures, the interference effect is somewhat washed out by the statistics. At low temperatures, however, the oscillation patterns become sharper, indicating strong interference. The BTE simulations do not show any such effect [Ren01].

Figure 1.27: The electron density along the channel of a DG-SOI MOSFET calculated by both the semi-classical BTE (dashed lines) and the quantum transport model (solid lines) [Ren01].

1.4

OTHER MULTI-GATE MOSFET STRUCTURES

In an ever increasing need for higher current drive and better SCEs immunity, MOSFETs are evolving from classical, planar, single-gate devices into three dimensional devices with multiple gates (double-, triple- or multi-gate devices). A brief review on some of these structures and their main features is given in this section. The first fabricated double-gate MOSFET was made in a tall and narrow Si island called fin [Dig90]. The FinFET structure (shown in Fig. 1.28) is very similar to that MOSFET expect for the presence of a thick dielectric layer called hard mask on the top of the Si fin [Xue99]. The hard mask is used to prevent the formation of parasitic inversion channels at the top corners of the device. The undesired effect of corners will be addressed later in this chapter.

Figure 1.28: FinFET Device Structure

One can expect that increasing the number of gates will increase the current-drive capability of the device and suppress the SCEs. When a 3rd gate is present in the MOSFET structure, we reach the triple-gate (or tri-gate) MOSFET [Lem04] shown in Fig. 1.29. The SCEs suppression can be improved by extending the sidewall portions of the gate electrode to some depth in the buried oxide (not shown in figure) which is called  -gate MOSFET [Jon01] or underneath the channel region which is called  gate [Lia02].

Figure 1.29: Triple-gate (trigate) MOSFET

22

CHAPTER 1: DOUBLE-GATE MOSFETS

The best possible control of SCEs can be achieved using the surrounding-gate MOSFET [Col08]. The first surrounding-gate MOSFET was fabricated by wrapping a gate electrode around a vertical Si pillar [Tak05]. More recently planar surroundinggate devices with square or circular cross sections have been reported [Pas06]. Multiple surrounding-gate channels can be stacked on the top of one another, while sharing common gate, source and drain [Sun04]. Such devices are called MultiBridged Channel MOSFET (MBCFET) and are used to increase the device currentdrive capability per unit area. Devices with a triple- or quadruple-gate structure present a non-planar Si/SiO2 interface with corners where premature inversion may form in the corners due to charge-sharing effects between two adjacent gates [Col08]. Consider the  -gate device shown in Fig. 1.30 with t si  Wsi  30 nm and oxide thickness of 2nm. Fig. 1.31 depicts the simulated dg m / dVG characteristic that is used to identify the threshold voltage. The humps of this curve correspond to the formation of channels in the device. For device A with rtop  rbot  1 nm , the lowest doping curve has a single hump which indicates that both corners and edges build up channels at the same time. With increasing the doping, two humps can be seen in the curve. The first hump corresponds to channel formation in the top corner while the second one corresponds to channel formation in bottom corners and edges. In contrast, when rtop  rbot  5nm (device B) only a single hump is identified in the curve even with higher doping. From the previous observations, one can deduce that corner effect can be suppressed by either lowering the doping concentration or using corners with larger radius [Col08]. In deed, the corner effect is not expected to pose any problem because the technology trend is to use undoped channels with metal gates and controlling the device threshold by adjustable metal gate work function [Jun04].

Figure 1.30:

 -gate device cross section. The thickness and width of the device are t si and W si .

The radius of curvature of the top and bottom corners is

rtop and rbot respectively.

Figure 1.31: dg m / dVG in the

 -gate MOSFET of Fig. 1.29. Gate is N+ polysilicon A:

rtop  rbot  1nm , B: rtop  rbot  5nm Development efforts into multi-gate transistors have been reported by AMD, Hitachi, IBM, Infineon, Intel, TSMC, Freescale, UC Berkeley and others and the ITRS predicts that such devices will be the cornerstone of sub-32 nm technologies [ITRS]. The primary roadblock to widespread implementation is manufacturability, as both planar and non-planar designs present significant challenges, especially with respect to lithography and patterning.

CHAPTER 2

QUANTUM TRANSPORT USING THE NONEQUILIBRIUM GREEN’S FUNCTION FORMALISM

Numerical simulation of semiconductor devices is an important complement to analysis and experiment. Device simulation provides essential capabilities absent in both analysis and experiment. The prediction of device behavior before fabrication is the most advantageous feature of device simulation over experimental measurement. Running a computer-based simulation costs nothing in comparison with fabrication and different types of measurements. Device simulation allows a spectrum of conditions (under which fabricated device may be damaged) to be tested inexpensively. Moreover, it is possible to simulate devices which aren’t yet manufacturabe. It can also explore the physical phenomena occurring in a semiconductor device and, thus, it improves our understanding of device operation and enables us to develop analytical models. Device simulation is one component of technology computer-aided design (TCAD), which provides the data needed by device compact models used by SPICE [Ant88] simulation. The relationship between various simulation design steps that have to be followed to achieve certain customer need is illustrated in Fig. 2.1 [Dra06]. Rapid device scaling pushes the dimensions of the field effect transistors to the nanometer regime [ITRS]. Scaling bulk MOSFETs to and beyond the 32 nm node technology will face fundamental limits [Fra01] and are expected to be replaced with ultra-thin body (UTB) SOI MOSFETs and Double-Gate (DG) MOSFETs. The International Technology Roadmap for Semiconductors 2007 (ITRS) projection for the DG MOSFETs physical gate length for high performance logic is 4.5 nm for the year 2022 [ITRS]. For extremely scaled dimensions, the quantum effects (illustrated in chapter 1) including carrier confinement along the direction normal to the Si-SiO2 interface, tunneling through the gate oxides, source to drain tunneling and quantum interference play an important role in determining the DG MOSFETs characteristics. These effects can be accurately predicted only using quantum mechanical based device simulation [Sve07].

Figure 2.1: Design sequence to achieve desired customer need

The non-equilibrium Green’s function formalism (NEGF) provides a rigorous description of quantum transport in nanoscale devices [Dat00]. For ballistic transport, NEGF is equivalent to solving Schrödinger equation. The NEGF, however, provides a powerful way for treating the boundary of 2D and 3D problems. Moreover, it provides a way for treating scattering within non-ballistic devices. The rigorous description for this formalism can be found in the literature [Mart59, Kad62, Kyl65, Lan76, Dan84, Ram86, Mah87, Khan87, Jau94, and Hau96] where it is described using an advanced language in the quantum mechanics world, namely, the second quantization language. A simpler description can be found also in the literature [Dat95, Dat00, Dat05, and Ana06]. This chapter reviews important properties of coherent transport, nano-scale device simulation involving this type of transport, the NEGF and finally computationally efficient methods used within the NEGF framework. The chapter begins by introducing the self-consistent field method in section 2.1. It is an iterative method used for solving the non-linear equations of the transport problem together with the electrostatic equation. In section 2.2, coherent transport through a single level conductor is introduced. It helps to understand coherent transport concepts applied to a simple device of one energy level. Then a more realistic situation is given in section 2.3 where coherent transport is applied to a one-dimensional (1D) device using the wavefunction approach. Section 2.4 generalizes this approach, by induction, and introduces the NEGF. Since the NEGF is computationally intensive, two computationally efficient methods are reviewed in section 2.5.

2.1

SELF-CONSISTENT FIELD METHOD

Device simulation process is based on solving two main problems, electromagnetic problem and carrier transport problem. At steady state, electromagnetic becomes electrostatic and treated by solving Poisson’s equation. Poisson’s equation describes the dependence of the electrostatic potential in the device on the distribution of fixed charges (dopant ions, for example) and free carriers

26

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

(electrons and holes) [Tau98]. For inhomogeneous material device, like MOSFETs for example, Poisson’s equation can be written as [Hay00]:



  V   q p  n  N D  N A



(2-1)

where V is the electrostatic potential, q is the electronic charge,  is the permittivity of the medium, p and n are the hole and electron concentrations respectively and

N D and N A are the ionized donor and acceptor doping concentrations respectively. Electrostatic and carrier transport problems are coupled together. To determine the electrostatic potential by Poisson’s equation, we need the carrier distribution. At the same time, the carrier distribution is obtained from transport equations which depend on the electrostatic potential. The solution of the NEGF and Poisson’s equation is carried out by the self-consistent field method [Str72, Dat00]. It is an iterative method starts by assuming an initial guess for the potential distribution in the device (see Fig. 2.2). According to this potential, the NEGF is used to calculate the electron and hole concentrations in the device. With the electron and hole concentrations are known, Poisson’s equation can be solved yielding a new potential distribution. The new potential is compared to the old potential and the solution cycle is repeated until self consistent solution for the potential is obtained. Self-consistency criterion is that the difference in potential between two successive iterations drops below a certain tolerance. Assume initial guess for the potential

V V old  V new

NEGF

V

old

 n, p

Poisson’s equation

n, p  V new No

V new  V old   Yes

Figure 2. 2: Flow chart illustrating the self-consistent field method used in device simulation using the NEGF

2.2

COHERENT TRANSPORT THROUGH A SINGLE ENERGY LEVEL

Coherent transport is a special type of carrier transport arises when device dimensions become smaller than the phase breaking length [Fer97]. In this case, the carrier energy and phase is kept constant and carriers’ wave nature plays an important role in the transport process. For example, electrons can interfere instructively or

destructively. Therefore, we can’t describe the transport in this case without including the wave nature of the carriers in our description. This nature can be included by solving the Schrödinger wave equation. Now we are going to discuss electrons coherent transport through a tiny device with single energy level. In this case, the problem is simple and can be described without going into Schrödinger equation [Dat05]. In section 2.3, transport through a 1D device will be described using Schrödinger equation. Contact 1

Single Level Device Contact 2 2 1  

EF 1

EF 2

qVb  EF1  EF 2

Figure 2. 3: Single level device sandwiched between two contacts.

Consider the simple system shown in Fig. 2.3 consisting of a single level device sandwiched between two contacts. The contacts are assumed to have constant electrostatic potential, relatively large size (reservoirs), extremely small levels spacing (almost continuous), strong enough scattering to maintain thermal equilibrium with constant Fermi energy E F . If no battery is connected between the contacts, the overall system will be in equilibrium with constant Fermi energy. Now assume the contacts are connected by a battery Vb . The battery disturbs the equilibrium situation by lowering the Fermi level in contact 2 by a qVb with respect to contact 1. The FermiDirac distribution function in contact 1 is given by 1 (2-2) f1 ( E )  1  exp[( E  EF 1 ) / K BT ] and in contact 2 1 (2-3) f2 (E)  1  exp[( E  EF 2 ) / K BT ] where K B is Boltzman constant and T is contact temperature in Kelvin. Due to the mentioned difference in Fermi energy, contact 1 pumps electrons into the device and contact 2 pulls out electrons from the device; a situation leads to a current flow under condition that E F 2    E F1 if the device energy level is denoted by  . The situation can be explained as follow [Dat05]: contact 1 would like to see f1 ( ) electrons in the device, while contact 2 would like to see f 2 ( ) electrons with f 2 smaller than f 1 . At steady state, the number of electrons in the device N will be something in between f 1 and f 2 . One can expect that the current flowing from contact 1 to the device is proportional to f1  N and the current flowing out of the

28

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

device to contact 2 is proportional to N  f 2 . We can write the current equations as [Dat05]: I1device  q I device2  q

1 

2 

[ f1 ( )  N ]

(2-4)

[ N  f 2 ( )]

(2-5)

where q is the electronic charge,  is Dirac’s constant,  1 and  2 are contact 1/device and contact 2/device coupling strength respectively. The factor

 is 

interpreted as the rate at which an electron placed in the device will escape into the contact [Dat05]. At steady state, I1device  I device2 and from which one can deduce that

 1 f1 ( )   2 f 2 ( ) 1   2 q  1 2 I [ f1 ( )  f 2 ( )]  ( 1   2 ) N

(2-6) (2-7)

One important thing is missing in the previous equations, the device energy level broadening due to coupling with the contacts [Dat05]. From uncertainty principle, there must be broadening (uncertainty) in the electron energy level as long as the electron has a finite life time in the device due to being able to escape to the contacts. When the device is isolated, it has a sharp density of states (DOS) given by the Dirac delta function  ( E   ) . The broadened DOS could have, in the simplest case, a Lorentzian function centered around E   with a unity area under the curve [Dat05]:

D ( E ) 

 / 2 (E   )2   2 / 4

(2-8)

where    1   2 . Therefore eq. (2-6) and (2-7) should be modified to 

N

 D ( E )

 

I

 1 f1 ( E )   2 f 2 ( E ) dE 1   2

q  1 2 D ( E ) [ f1 ( E )  f 2 ( E )]dE    ( 1   2 )

(2-9) (2-10)

The effect of the potential wasn’t yet included in the simple transport model mentioned above. The potential in this tiny device modulates the density of states by shifting the energy level location. The DOS should be modified to D ( E  U ) where U is the potential energy given by  qV and V is the potential at the device. The potential can be obtained by solving Poisson’s equation given the number of electrons

N using eq. (2-9). The DOS, and as a result the number of electrons, is potential dependent and therefore an iterative solution is needed to find a self-consistent potential value. Based on the single level device description, one may think that the single level DOS D ( E  U ) should be, only, replaced by another DOS in case of a multilevel device and that’s it. Let’s think of it step by step and see if modifying the DOS only will serve it or not. Since we have more than one level, the level energy  should be replaced by a column vector that contains the energy for each level. The transfer rate (or the coupling strength) may not be identical for each energy level also and thus, we need another column vector to describe  . By the same way we can conclude that all the scalar quantities have to be replaced by column vectors. A more general picture may turn all the scalar quantities into matrices with the mentioned vectors on the diagonal of these matrices. In fact, these matrices may contain also off-diagonal terms. This comes to the heart of quantum mechanics; there is no single representation that diagonalizes all matrices [Dat05]. Therefore, for a multi-level device with l levels, each of the scalar quantities should be replaced by matrices of size l  l as given in Table 2-1. A simple derivation of these matrices will be given in section 2.4 with their physical meaning. Table 2- 1: Single level device scalar quantities and their corresponding matrices for a multi-level device Single level device



 1, 2 2D( E ) 2n( E )

N

2.3

COHERENT

Multi-level device

H : Hamiltonian Matrix Γ1, 2 : Broadening Function

A(E ) : Spectral Function

G n (E ) : Correlation Function dE n ρ G ( E ) : Density Matrix 2

TRANSPORT

THROUGH

A

1D

DEVICE Consider a 1D device connected to two contacts as shown in Fig. 2.4. Again the contacts are assumed to have constant electrostatic potential, relatively large size, extremely small levels spacing, strong enough scattering to maintain thermal equilibrium with constant Fermi energy E F . As shown in Fig. 2.4, contact 1 injects a flux of carriers into the device. Some carriers are reflected from potential barriers within the device with a reflection coefficient r . The rest of carriers are transmitted across with a transmission coefficient t and enter contact 2. The contacts are assumed to absorb the carriers with no reflection [Dat95]. The same scenario occurs with the carriers injected from contact 2. That is to say two sources of carriers exist, contact 1 and contact 2.

30

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM Contact1

Contact2

1D Device L

0

eik1x reik1x

te ik1 ( x  L )

 D(1)

te ik2 x

 D( 2)

eik2 ( x L) reik2 ( xL)

 D( B ) Figure 2. 4: A 1D device connected to two contacts. All wave functions in the device can be represented as left incident D , right incident  D and bound state  D . (1)

( B)

( 2)

Now according to this situation, there are three kinds of the wavefunctions in the device [Ana06]: 1.  D(1) : Wavefunction due to injection from contact 1. 2.  D( 2 ) : Wavefunction due to injection from contact 2. 3.  D( B ) : Bound and quasi-bound states which are filled by scattering mechanisms. These states will be neglected as we are considering coherent transport of carriers It is not correct to think that  D(1) and  D( 2 ) should be added to obtain the total wavefunction within the device. One can’t add two wavefunctions originating from two different sources because they have no phase relation connecting them [Dat05]. The Landauer-Buttiker approach expresses the expectation value of an operator in terms of  D(1) ,  D( 2 ) and the contacts’ distribution functions. The expectation value of 

operator O is given by [Ana06]: 



O  2[( D(1) )* O D(1) ] f1 ( E )  2[( D( 2) )* O D( 2) ] f 2 ( E ) (2-11) k1

k2

where the factor 2 accounts for spin. The electron density within the device as a function of position x is given by [Ana06] n( x)  2|  D(1) ( x; k1 ) | 2 f1 ( E )  2|  D( 2) ( x; k 2 ) | 2 f 2 ( E ) k1

The summation over

(2-12)

k2

k can be transformed to integration using periodic boundary

conditions in the contacts

1

  2  dk [Pie03]: k

dk1 1 dk 1 | D(1) ( x; E ) |2 f1 ( E )dE  2 2 | D( 2) ( x; E ) |2 f 2 ( E )dE dE 2 dE 2 (1)  2 LDOS ( x; E ) f1 ( E )dE  2 LDOS( 2) ( x; E ) f 2 ( E )dE (2-13)

n( x )  2 

where LDOS( x; E ) is the local density of states. It expresses how the density of states change both with energy and position and hence the name local. At equilibrium f1 ( E )  f 2 ( E ) and the LDOS is filled according to one Fermi function. Under bias, we fill up part of it, LDOS (1), according to the Fermi function in contact 1 and the other part of it, LDOS (2), according to the Fermi function in contact 2. This concept allows us to treat a non-equilibrium problem using equilibrium statistical mechanics [Dat00]. The electrical current density is given by [Ana06]: J ( x)  2

 q  d (1) ( x; k1 ) d (1) ( x; k2 ) [ D(1) ( x; k1 )]* D f1 ( E )  [ D( 2) ( x; k2 )]* D f 2 ( E ) (2-14)   2im  k1 dx dx k2 

The current is continuous and, hence, we can substitute at x  0 for example to obtain that:  q  2 2 (2-15) J 2  ik1 (1  r ) f1 ( E )   ik 2 t f 2 ( E )   2im  k1 k2  Current continuity (current entering the device is equal to that is getting out of it) imposes that [Lib80] v2 t

2

integrations and using v 

 v1 (1  r ) , transforming summations to 2

k 1 dE , eq. (2.15) can be rewritten as:  m  dk

2q v 2 2 t  f1 ( E )  f 2 ( E )dE h  v1 2q  2  T ( E ) f1 ( E )  f 2 ( E )dE h

J 2

(2-16)

where T (E ) is called the transmission function. Now to be able to calculate the electron density and electrical current, D(1) ( x; E ) ,  D( 2) ( x; E ) and T (E ) must be calculated first. They are obtained by solving Schrödinger equation, one time for incident waves from contact 1 to find  D(1) ( x; E ) and another time for incident waves from contact 2 to find  D( 2) ( x; E ) . The 1D steady state Schrödinger equation for a single parabolic conduction band is given by [Dat89]:  2 d 2   EC ( x)  E (2-17) 2m* dx 2 where m * is the effective mass and EC (x) is the conduction band edge along the device

32

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

Generally, the conduction band edge along the device is non-uniform and a numerical solution of Schrödinger equation is needed. Applying a central finite difference approximation for 2nd order derivative in (2-17) and denoting  ( xi ) by  i as shown in Fig. 2.5, then:

 2 i  i 1  2 i  i 1  x 2 x 2

1  2 1





2

(2-18)

x

i



N

N

x 1

x

Nx 1 N x

i

x =L

x =0

Figure 2. 5: A 1D discretized device connected to two semi-infinite contacts leads.

Substituting with eq. (2-18) into eq. (2-17), we obtain:

 t x i 1  2t x i  t x i 1  EC i i  E i where t x 

(2-19)

2 is the coupling energy between adjacent sites on the grid. 2m * (x) 2

The boundary condition at nodes i  1 and i  N x should be treated carefully to account for the scenarios shown in Fig. 2.4. Let’s ignore the boundary condition for the moment by assuming closed boundary condition; i.e. 0   N x 1  0 , and correcting this assumption after just few lines. Writing eq. (2-19) at each node along 1D grid in the x-direction, a set of linear equations are obtained and can be cast in the matrix form: Hψ  Eψ  tx 0  0  2t x  EC1  t 2t x  EC 2  t x 0   x  where    0  H  0   0  t x 2t x  ECNx 1  tx     0  0  tx 2t x  ECNx  

(2-20)

N x N x

 1     2    ψ      N x 1     Nx  

Note the structure of the Hamiltonian matrix; it consists of diagonal and off-diagonal terms. Each element on the diagonal is related to a certain point in the grid. Offdiagonal terms represents coupling between the points. It makes it easier to remember this structure when dealing with the quantum mechanics matrices. Eq. (2-20) is an eigenvalue problem. The solution of such a problem yields the eigenstates and eigenvalues of the device when it is isolated from the outside world; i.e. when the system is closed. In contrast, the transport problem involves a transfer of

electrons between the device and the contacts. Now let’s see how these matrices should be modified to account for injection from either contacts as shown in Fig. 2.4.

Injection from contact 1 ( D ) (1)

For this case, we can write that  ( x)  e ik1x  re ik1x x  0 Then  D(11)   (0)  1  r  r   D(11)  1 and  D(10)   (x)  e ik1x  reik1x  e ik1x  ( D(11)  1)eik1x Thus, when i  1 in eq. (2-19),  D(10) should be substituted by eik1x  ( D(11)  1)eik1x (1) instead of zero. In a similar way when i  N x , it can be deduced that  DN should be x 1

(1) e ik2x instead of zero in eq. (2-19). substituted by  DN x

Therefore, in order to account for the open boundary, eq. (2-20) must be modified as follows:  EI - H - Σ1  Σ 2 ψ (1) D  iΣ 1  Σ 1 

(2-21)

where I is the identity matrix, Σ 1 is called the self-energy of contact 1 and given by:  t x eik1x 0   0 0  Σ1           0

0     0 N  N x x

(2-22)

and Σ 2 is called the self-energy of contact 2 and given by:

0  Σ2     0

      0       t x eik2x  N  N x x  

0

(2-23)

Injection from contact 2 ( D ) ( 2)

The analysis is very similar to the case of injection from contact 1 and yields that:  EI - H - Σ1  Σ 2 ψ (2) D  iΣ 2  Σ 2 

(2-24)

34

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

Now it is time to mention very important properties of the self-energy. The selfenergy Σ has two great differences that distinguish it from normal potential energy encountered in quantum mechanics [Dat95]. The first difference is that Σ is energy dependent. Because of that, it is more convenient to think of the energy as an independent variable and we want to know the response of the device to the incident electrons with different energies. This is always the case for an open system rather than for the isolated one which has certain eigenenergies. The second difference between Σ and normal potential is that Σ is not Hermitian. This makes the effective Hamiltonian H  Σ  has complex eigenvalues and the imaginary part of these eigenvalues broadens the density of states and gives the eigenstates a finite life time. Finite life time and energy level broadening are conceptually related to each other through uncertainty principle. The energy dependence of Σ 1 and Σ 2 comes from the wave numbers k1 and k 2 . The wave numbers are obtained at a given energy from the E (k ) relation of contact 1 and contact 2. The E (k ) relation of contacts will be mentioned hereinafter. Eq. (2-21) and (2-24) have their solution in the form:

ψ (1) D ( E )  -EI - H - Σ1  Σ 2  Γ 1

(2-25)

ψ (2) D ( E )  -EI - H - Σ1  Σ 2  Γ 2

(2-26)

-1

-1

where the broadening functions Γ 1 and Γ 2 are given by:





(2-27)





(2-28)

Γ 1  i Σ1  Σ1 Γ 2  i Σ 2  Σ 2

Physically the broadening function determines the electron exchange rates between the contact and the device [Dat05]. It is responsible for the broadening of the energy levels and the states finite life time. If the self-energy matrices are purely real then there is no exchange rate between the device and the contacts and the effective Hamiltonian is Hermitian with real eigenvalues. In this case, the real self-energy just shifts the eigenvalues of the device as if it was a normal potential.

2.3.1 Discrete Leads E (k ) Relation For a continuous 1D lead, the potential is constant and, thus, the lead eigenfunctions are plane waves with a parabolic dispersion relation [Dat95]:

2k 2  k ( x)  e , E ( k )  U 0  2m* ikx

(2-29)

where U 0 is the potential in the lead. Now for a discrete lead one can expect that the lead eigenfunctions are also plane waves given by  n  eiknx . Substituting back into

the discretized Schrödinger equation (same as eq. (2-19) but with EC replaced by U 0 ), the dispersion relation is obtained [Dat95]: E(k )  U 0  2t x (1  cos(kx))

(2-30)

The velocity can be also obtained v(k ) 

1 E  2xt x sin(kx)  k

(2-31)

Actually, there is a range of energies in which the discrete representation is accurate. This range can be defined by measuring how the E (k ) relation of the discrete representation is close to the continuous limit. Fig. 2.6 depicts E (k )  U 0 normalized to t x , obtained from both using eq. (2-29) and (2-30), versus the wave number k normalized to  / x . It is noticeable that, both equations are in close agreement for the range ( E (k )  U 0 ) / t x  1 . Thus it is important to make sure that the energy range of interest falls within this range. If not, we have to decrease the grid spacing to increase the factor t x which, in turn, allows the simulation of larger range of energies with acceptable accuracy.

Figure 2. 6: Continuous and discrete E (k) relations plotted on the same graph. The discrete relation is in close agreement with the continues for the range E (k )  U 0  t x

Now let’s summarize the solution steps for the transport problem shown in Fig. 2.4: 1. Guess the conduction band edge profile ECi along the device 2. For each energy E :

36

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

a. Calculate k1 and k 2 using the lead dispersion relation (2.30) as follows:

k1 

 E  EC 0   E  ECNx 1 1  , k2  cos 1 1  cos 1 1  x 2t x x 2t x   

  

b. Calculate Σ1 and Σ 2 using eq. (2-22) and (2-23) then Γ 1 and Γ 2 using eq. (2-27) and (2-28) (2) c. Calculate ψ (1) D (E) and ψ D (E) using eq. (2-25) and (2-26) 3. Calculate ni using eq. (2-13) where LDOS is given by: LDOS(1) (i; E )  2t x  sin(k1x) | D(1) (i; E ) |2

LDOS( 2) (i; E)  2t x  sin(k2 x) |  D( 2) (i; E) |2 4. Solve Poisson’s equation to find ECi 5. Go to step 2 and repeat until self consistency is achieved (follow chart shown in Fig. 2.2) 6. Calculate the current using eq. (2-16) where T (E ) is given by T (E) 

v2 sin(k 2 x) 2 1) t (E)   ( E ) (DN x sin(k1x) v1

2

(2-32)

Note: the LDOS given above are obtained from eq. (2-13) divided by x . In a discrete representation the wavefunction is normalized to the number of points N x . The real wavefunction in a continuous representation, however, is normalized to the length of the domain N x x .

2.4

FROM

WAVEFUNCTION

TO

GREEN’S

FUNCTION In this section we describe how one can reach the NEGF equations starting from the wavefunction approach described in the previous section. This isn’t a formal derivation of the non-equilibrium Green’s function formalism. The rigorous description for this formalism can be found in the literature [Mart59, Kad62, Kyl65, Lan76, Dan84, Ram86, Mah87, Khan87, Jau94, and Hau96]. However it is described there using an advanced language in the quantum mechanics world, namely, the second quantization language. A simpler description can be also found in the literature [Dat95, Dat00, Dat05, and Ana06]. Here we are following the latter description.

2.4.1 The Green’s Function The Green's function is named after the British mathematician George Green, who first developed the concept in the 1830s [Wik08].The concept of Green’s function appears in many physical contexts; whenever the response of a given system R is related to the excitation S by a differential operator Dop such that : Dop R  S . The Green’s function in 1D is the defined by the following equation [Wik08]:

DopG   ( x  x' ) It is interesting to notice that the Green’s function G( x, x' ) represents the response of the system at point x due to a unit impulse excitation at a point x' . It is important to highlight that the inverse of a differential operator, or equivalently the Green’s function, is not uniquely specified till the boundary condition is specified [Dat95]. As a result, two different Green’s functions exist, the retarded and the advanced Green’s function. When solving Schrödinger equation, the retarded Green’s function corresponds to outgoing waves from the source while advanced function corresponds to incoming waves far from the source [Dat95]. The retarded Green’s function of eq. (2.21) and (2.24) is given by: -1 G(E)  EI - H - Σ1  Σ 2  (2-33) Thus eq. (2-25) and (2-26) can be rewritten as: ψ (1) D ( E )  G ( E ) Γ1 ( E )

(2-34)

ψ (2) D ( E )  G ( E ) Γ 2 ( E )

(2-35)

2.4.2 The Spectral Function To proceed with the NEGF, we have to consider some generalizations. First of them is the spectral function which is the generalization (matrix version) of the LDOS [Dat05]. Consider the following generalization: 1 dk1 1 dk1  D (i; E ) D* (i; E )  ψ D ( E )ψ D ( E ) 2 dE 2 dE

where

  D1 D* 1  D1 D* 2   D 2 D* 1  D 2 D* 2   ψ D ( E )ψ D ( E )      * *  DN x D1  DN x D 2

*    D1 DN x  *   D 2 DN x      *   DN x DN  x 

dk ψ D ( E )ψ D ( E ) which leads, with dE some mathematical manipulation, to the formal expression [Dat05]:

The spectral function can be defined as: A(E) 

38

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

A(E)  G(E) Γ (E)G  (E)

(2-36)

The LDOS is given by the diagonal elements of the spectral function divided by 2

2.4.3 The Correlation Function The density matrix is the generalization (matrix version) of the electron density [Dat05] and obtained by replacing the LDOS with the spectral function in eq. (2-13): ρ  2

1 1 1 n A1 f1 ( E )dE  2 A2 f 2 ( E )dE  2 G ( E )dE (2-37) 2 2 2

where G n (E ) is called the correlation function and given by: G n ( E)  A1 ( E) f1 ( E)  A2 ( E) f 2 ( E)

(2-38)

The density matrix, in real-space (RS), ρi, j, t  is a generalization of the electron density, because it accounts for the correlation between different points, in space, at certain time t. It is more interesting to know that the correlation function G n i, j; t , t ' takes into account the correlation between different states and at two different times t and t’. At steady state the electron correlation function is Fourier transformed to its energy dependent version G n i, j; E  which is equivalent to the one in eq. (2-38). Then, the energy resolved electron density is given by the diagonal of the electron correlation function divided by 2 . A closely related quantity to G n is the hole correlation function G P . The hole correlation function is related to the empty states by the same way the electron correlation function is related to the filled states. Interested reader may consult [Dat95] for an insightful discussion on the difference between the carrier density, the density matrix and the correlation functions.

2.4.4 The Scattering Function Substituting for the spectral function using eq. (2-36) into eq. (2-38) yields the kinetic equation [Dat95]: G n ( E )  G( E ) Σ in ( E )G  ( E )

(2-39)

where Σ in E   Γ1 E  f1 E   Γ 2 E  f 2 E  is called the in-scattering self-energy. It can be viewed as a sum of two terms: Σ in  Σ1in  Σ 2in . The in-scattering self-energy in the device arises from coupling of the device to the external contacts. The term Σ1in represents in-scattering rate of electrons from contact 1 to the device assuming that the first point of the device is empty [Ana06]. A similar statement applies to Σ 2in . One more kinetic equation exists in the form:

G p ( E )  G( E ) Σ out ( E )G  ( E )

(2-40)

where Σ out is the out-scattering self-energy. The out-scattering self-energy is closely related to the in-scattering one, but the first represents out-scattering rate of electrons from the device to contact 1 (or 2) assuming that the first (or last) point of the device is occupied [Ana06]. One may wonder what the importance of eq. (2-39) and (2-40) is; they seem useless. In fact, these two equations are central equations in the NEGF [Dat95]. In the coherent transport case, however, they are not used. They came into play when non-coherent transport is considered. They play the role of the Boltzman Transport Equation (BTE) in the semiclassical treatment [Lun00]. For Different scattering mechanisms, in-scattering and out-scattering self-energies exist and must be taken into account to solve eq. (2-39) and eq. (2-40) for the correlation functions [Dat95]. Non-coherent scattering is out of the scope of this thesis. Interest reader can find more details in [Dat95, Dat00, Dat05, and Ana06].

2.4.5 Transmission Function Now we want to express the transmission function T (E ) in terms of the Green’s function. The quantities in the T (E ) expression given by eq. (2-32) can be replaced as follows [Ven03]: sin(k1x)  Γ1 (1,1) / 2t x , sin(k2 x)  Γ 2 ( N x , N x ) / 2t x 1)  ( E )(DN

2 x

 G ( N ,1; E ) Γ1 (1,1; E ) Γ1 (1,1; E )G  ( N ,1; E )

Accordingly, we can deduce that : T ( E)  Γ1 (1,1; E )G( N ,1; E ) Γ 2 ( N x , N x ; E)G  ( N ,1; E ) which can be rewritten as:



T ( E )  Trace Γ1 GΓ 2 G 

2.5



(2-41)

COMPUTATIONALLY EFFICIENT NEGF-BASED METHODS

The retarded Green’s function is a central quantity in the NEGF. As seen from eq. (2-33), Green’s function is calculated by matrix inversion for effective Hamiltonian matrix. Hamiltonian matrix size, in real-space (RS), is the same as the number of points in the grid. Number of gird points N  N x N y for 2D and N  N x N y N z for 3D devices. Numerical matrix inversion consumes a large number

of operations that is proportional to N 3 [Car95]. Moreover, Green’s function should be calculated for each energy point considered in the simulation. This makes total number of operations scales as N E N 3 where N E is the number of energy points. Therefore, efforts have been exerted to develop computationally efficient methods to reduce the computational burden. Based on Dyson’s equation [Hau96, Ale03], two methods exist: the recursive Green’s function (RGF) algorithm [Kli98, Kli02, Sve02,

40

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

Ana06 and Ahm07] and the contact block reduction (CBR) method [Mam03, Sab03, Sab04, Hack02, Mam04, Mam05, Kha06 and Mam07]. The RGF algorithm builds up the Green’s function recursively without full inversion of the Hamiltonian matrix. It can be used in the presence of scattering as well as in the ballistic limit for a layered device [Kli98, Sve02, and Ana06]. The drawback of the RGF is that it can be used only if the effective Hamiltonian matrix EI  H  Σ  is block tri-diagonal. This means it allows only nearest neighbor layers coupling in RS. Unfortunately, a lead couples all the layers connected to it [Mam03] and, therefore, the RGF works when the device has no more than two contacts (see Fig. 2.7 for illustration). NEMO 1-D [Kli98, NEM] is the first NEGF-based TCAD tool and it uses the RGF algorithm to reduce the computation and memory needed to carry out 1D simulation of the resonant tunneling devices. Recently, NEMO 1-D has been extended to the Quantum Dot Lab [Kli02] simulator available on nanohub: www.nanohub.org . A newly developed tool, NanoFET [Ahm07] uses also the RGF algorithm to simulate 2D planar DG MOSFETs devices. The CBR method calculates the coherent transport transmission function of an arbitrarily shaped, multi-terminal 2D or 3D structure in an efficient way [Mam03]. The method has been extended to calculate the electron density in [Mam04] and the overall CBR method has been applied to 10nm gate length bulk n-MOSFET. The ability to calculate the electron density in an efficient way is necessary for efficient fully quantum mechanical self-consistent calculation. Recently, the CBR method has been used to analyze and optimize the behavior of 10nm FinFET device in the quasiballistic regime of operation in [Kha06] and the results were compared to experimental ones [Mam07]. Contact 3

Contact 1

Contact 2

Figure 2. 7: A layered structure connected to 3 contacts. Contact 1 is coupled to the left most layer, contact 2 is coupled to the right most layer and contact 3 is coupled to all the layers and is coupling them together. The RGF can’t be used as along as contact 3 is coupled to the device.

2.5.1 Recursive Green’s Function Algorithm Within the NEGF framework, charge density is obtained by integrating the main diagonal of the electron correlation function G n and current density by integrating the diagonals adjacent to the main diagonal [Ana06]. As a result, the entire G n matrix isn’t needed to calculate the charge and the current density. Suppose that the simulated system is composed of N x layers and divided into left subsystem L and right subsystem R as shown in Fig. 2.8. Dyson’s equation relates the Green’s function of the full system L +R in terms of the Green’s function of the subsystems L,

R and the coupling between them. Let the full system Green’s function denoted by G , it satisfies the following equation: DG  I

(2-42)

 DLL DLR  G G  where G and D have the forms  LL LR  and   respectively.  DRL DRR  GRL GRR 

Figure 2. 8: The overall system can be divided into two parts left L and right R

Imagine that the left and right subsystems were isolated (i.e. there is no coupling between them) then, the overall system Green’s function is denoted by G 0 and satisfies the following equation: D0G 0  I

(2-43)

G 0 LL 0   DLL 0  where G and D have the forms  and  0 D  respectively. 0 RR   0 G RR  0

0

Dyson’s equation relates G in terms of the G 0 and the coupling matrix D0  D and is given by [Ana06] G  G 0  G 0UG  G 0  GUG0

where

(2-44)

 0  DLR   U    D 0 RL  

A similar equation exists for the electron correlation function [Ana06] G n  G n0  G n0U G   G 0UGn  G n0  GUGn0  G nU G 0 (2-45)

where G n0  G 0 Σ inG 0 The RGF algorithm consists of forward recursion and backward recursion steps

42

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

2.5.1.1 Forward Recusrion Left connected Green’s function g L of a group of layers is defined as the Green’s function taking into account all the layers to the left but assuming no layers exist on the right; i.e. it takes into account the coupling to the left layers only. The forward recursion algorithm build up the left connected Green’s function of the full system recursively. The forward recursion formula is derived by assuming that we already have the left Lq 1

Lq

connected Green’s function matrix block g q,q and we want to obtain gq 1,q 1 (see Fig. 2.9). The overall system composed of the q+1 layers is divided into left subsystem and right subsystem. The left subsystem is composed of q layers on the left and the right subsystem is the right most layer. Using eq. (2-42), (2-43) and (2-44), it can be shown that [Sve02, Ana06]



gqLq1,1q 1  Dq 1,q 1  Dq 1,q gqLq,q Dq ,q 1

g Lq  D1:q1,1:q

Dq11,q 1



1

(2-46)

g Lq1

(a) (b) Figure 2. 9: (a) The left connected Green’s function of the left most q layers (q=7 for example) and the isolated Green’s function of the q+1th layer are given. (b) We seek the left connected Green’s function of the overall q+1 layer.

g1L,11

g 2L,22

g 3, 3

L3

g NLNx  GN x , N x x ,Nx

Figure 2. 10 The forward recursion sequence. Initially, the left connected green’s function of the first L1

layer g1,1 is calculated. Then, the first two layers left connected green’s function block and so on until the left connected Green’s function block

g 2L,22 is obtained

g NLNx ,xN x of the overall system is obtained.

Initially, the left connected green’s function of the first layer g1L,11 is calculated by

g1L,11  D1,11 . Then, using the above recursion formula, the first two layers left connected green’s function block g 2L,22 is obtained in terms of g1L,11 . By the same way

g3L,33 is obtained in terms of g 2L,22 and so on until the left connected Green’s function block g NLNx ,xN x of the overall system is obtained. It is equal to the last block of the fullyconnected Green’s function of the full system GN x , N x . The forward recursion sequence is shown schematically in Fig. 2.10. In a similar way, using eq. (2-42), (2-43), (2-44) 1 and (2-45), gqnLq 1, q 1 is recursively obtained as [Ana06]







1 Lq 1 in nLq  Lq 1 gqnLq (2-47) 1, q 1  gq 1, q 1 Σ q 1, q 1  Dq 1, q gq , q Dq , q 1 gq 1, q 1

2.5.1.2 Backward recursion The backward recursion algorithm builds up the fully connected Green’s function of the full system recursively. It starts from the right most corner block of the fully connected Green’s function matrix G N x , N x which was obtained from the forward recursion. The backward recursion formula is derived by assuming that we already have the block Gq 1,q 1 of the fully connected Green’s function matrix and we want to obtain Gq,q and Gq 1,q . The overall system is divided into left subsystem, composed of q layers on the left, and right subsystem , composed of the rest of layers from q+1 to N x . Using eq. (2-42), (2-43) and (2-44), it can be shown that the diagonal and adjacent diagonals elements of the fully connected Green’s function are given by [Sv202, Ana06]:

Gq ,q  gqLq,q  gq ,q Dq ,q 1Gq 1,q

(2-48)

Gq 1,q  Gq 1,q 1 Dq 1,q g Lq q,q

(2-49)

Gq,q 1   g Lq q,q Dq,q 1Gq 1,q 1

(2-50)

Lq

The right most corner block of the fully connected Green’s function GN x , N x is, already, calculated from the forward recursion. The backward recursion starts by calculating GN x 1, N x using eq. (2-50). Then GN x , N x 1 is calculated using eq. (2-51). The next diagonal block GN x , 1N x 1 is obtained using eq. (2-52). The same procedure is repeated until the diagonal, upper and lower off diagonal blocks of the fully connected Green’s function are obtained. The backward recursion sequence is shown schematically in Fig. 2.11. Similarly, using eq. (2-42), (2-43), (2-44) and (2-45), Gqn1,q and Gq,n q can be written in terms of Gqn1,q 1 as [Ana06]: n   Lq Gqn1,q  Gq 1,q1 Dq 1,q gqnLq , q  Gq 1, q 1 Dq 1, q gq , q

(2-51)

44

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM







Lq n   Lq nLq   nLq Gqn,q  gqnLq , q  gq , q Dq ,q 1Gq 1, q 1 Dq 1, q gq , q  gq , q Dq , q 1Gq 1,q  Gq , q 1 Dq 1,q gq , q



(2-52)

G1,1

G     



Gq1,q



2 ,1

Gq ,q  Gq1,q

 

GN x 3, N x 2 GN x  2 , N x  2

GN x 2, N x 1

GN x 1, N x 2

GN x 1, N x 1 G N 1, N x x G N x , N x 1

GN x , N x

     

Figure 2. 11: The backward recursion sequence. Initially, G N x , N x is used to calculate GN x 1, N x ,

GN x , N x 1 and GN x , 1N x 1 . The same procedure is repeated until the diagonal the upper and lower off diagonal blocks of the fully connected Green’s function are obtained.

To summarize, the algorithm to obtain the three diagonals of G , G n and G p [Ana06]:

is

For G 1

1) g1L,11  D11 2) For q=1,2,…, Nx-1, compute eq. (2-46) 3) For q=1,2,…, Nx-1, compute gq,Lq q 4) GN x , N x  g NLNx ,xN x 5) For q= Nx-1, Nx-2,…, 1, compute eq. (2-50), (2-49) and (2-48) in this order 6) For q=1,2,…, Nx-1, compute Gq,q 1 and Gq1,q For G 1) 2) 3) 4) 5)

n

and

Gp 

nL1 g11  g11L1 Σ11in g11L1 For q=1,2,…, Nx-1, compute eq. (2-47) x GNn x , N x  g NnLN x ,Nx For q= Nx-1, Nx-2,…, 1, compute eq. (2-52)and (2-51) Use Gqn,q 1  Gqn1,q





6) Use G p  i G  G   G n Then, the electron density is given by the diagonal elements of the integrated correlation function n

1 Diag (G n (E ))dE 2 

(2-53)

The current density flowing from the left contacts into layer 1 of the device is given by [Ana06] I 2





q 1 Trace i[G1,1 ( E )  G1,1 ( E )] 1in ( E )  G1n,1 ( E )1 ( E ) dE (2-54)   2

For 2D devices, the operation count of this algorithm scales approximately as N y3 N x [Ana06]. The dependence on N y3 arises because matrices of the sub Hamiltonian of layers should be inverted (recall eq. (2-46)), and the dependence on N x corresponds to one such inversion for each of the layers. These operations are carried out for each energy step and, therefore, the total number of operations is estimated as N E N y3 N x .

2.5.2 Contact Block Reduction Method The CBR method allows one to calculate 1D, 2D or 3D ballistic transport properties in a device that may have any number of leads. There are three key points in the CBR method: 1. Dyson’s equation [Hau96, Ale03] is used together with a clever splitting of the simulation domain. This makes elements of the Green’s function can be calculated with inversion of a relatively small matrix [Mam03]. Moreover, the electron density and terminal current can be calculated from only a small part of the Green’s function matrix. Thus, there is no need to find all the elements of the Green’s function matrix [Mam04]. 2. Decoupled (isolated) device eigenstates are used as a basis for the open boundary problem. The number of eigenstates needed to maintain acceptable accuracy is greatly reduced by applying von Neumann boundary condition for the decoupled device Hamiltonian [Mam03]. 3. Use of the leads propagating modes as a basis instead of the real-space basis in single band case [Mam05]. The basic idea is that only propagating modes contribute to the current and, neglecting non-propagating modes has a neglected effect on the Green’s function. In this section, we are going to briefly explain the method stressing on the key points mentioned above.

2.5.2.1

The Green’s Function

The Greens’ function of the isolated device G 0 is given by the spectral representation [Dat95]: G 0 (i, j , E )   

ψ (i )ψ* ( j ) ,  0 E     i

(2-55)

46

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM

where ψ are the eigenfunctions of the isolated device,   are the corresponding eigenenergies, i and j are the grid point’s indices in real-space. The Green’s function of the device can be calculated using Dyson’s equation [Hau96, Ale03]: G  X 1G 0

(2-56)

where the matrix X is given by X  I  G0 Σ

(2-57)

The simulation domain is spitted into two sub-domains, D: the interior part of the device and C: the boundary region that connects the interior parts to the contacts. Accordingly, the isolated device Hamiltonian matrix can be written as [Mam03]:  H0 H 0   0C  H DC

0  H CD  H D0  N  N

(2-58)

where H C0 NC  NC is a relatively small matrix corresponds to N C boundary points and H D0 N D  N D is a huge matrix corresponds to the N D interior points such that N C  N D  N where N the total number of points in the grid. Similarly, the isolated device Green’s function can be written as [Mam03]:

 G0 G 0   0C GDC

0  GCD  GD0  N  N

(2-59)

As a result of the above splitting, the self energy matrix is given by [Mam03] Σ Σ  C 0

0 0  N  N

(2-60)

Substituting with eq. (2-60) into eq. (2-57), the matrix X is given by

 X X  0C  GDC Σ C



where X C  I  GC0 Σ C



NC  NC

0 I 

(2-61)

is a relatively small matrix.

Substituting into eq. (2-59), the Green’s function of the coupled device is obtained in terms X C-1 which is an inverse of a small matrix:  X C1GC0 G 0 1 0 0 GDC Σ C X C GC  GDC

0   GC GCD  X C1GCD   (2-62) 0 1 0 0  GDC Σ C X C GCD  GD  GDC GD 

2.5.2.2

The Transmission Function

For ballistic transport, the energy resolved terminal current flowing between contact 1 and contact 2 is given by: I 12 ( E )  2





2q f 1  f 2 T12 ( E ) h

where f 1 and f 2 are the Fermi Dirac distribution function inside contacts and

T12 E  is the transmission function given by





Tλ1 λ 2  Trace Γ λ1 GΓ λ2 G  , 1  2

(2-63)

where the broadening function is given by



Γλ i Σλ  Σλ





(2-64)

Due to the sparse nature of the broadening function matrices, only a small part ( GC ) of the coupled Green’s function is needed in order to calculate the transmission function. Using eq. (2-62) and (2-63), the transmission function is given by [Mam03]



Tλ1 λ 2  Trace Γ Cλ1 GC Γ Cλ2 GC







(2-65)





where

Γ Cλ  i Σ Cλ  Σ Cλ

2.5.2.3

The Spectral Function

(2-66)

The spectral function filled by the λth contact, A is given by

Aλ  GΓ λG  (2-67) Substituting eq. (2-62) into eq. (2-67), the general form for the spectral function matrix element is obtained [Mam04]: Aλ i, j  

G

l,kC

0 i,l

B

1 C



Γ Cλ BC1

  G  

l,k

0  C k, j

(2-68)

where BC  I  Σ CGC0 Substituting eq. (2-55) into eq. (2-68) we obtain [Mam04]: Aλ (i, j; E)   ψ α (i)ψ β (j) *

α,β



 

Trace ψ β ψ α BC1 Γ Cλ BC1



(E -    i )(E -    i )

 ,η  0 (2-69)

The double summation on the eigenstates is composed of two terms, the 1st one is energy independent but position dependent and the 2nd one is the opposite. Therefore,

48

CHAPTER 2: QUANTUM

TRANSPORT USING THE NON-EQUILIBRIUM GREEN’S FUNCTION FORMALISM 2 2 N E  Neigen N [Mam05] the number of operations can be estimated as Nop  Neigen

where N eigen is the number of eigenstates to be used. 2.5.2.4

Von Neumann Boundary Condition

Isolated device Hamiltonian is constructed using closed boundary condition (Dirichlet boundary condition). The closed boundary condition assumes zero wavefunction at virtual points just outside the device; the boundary points have very small non-zero wavefunction. Therefore, it is possible to expand physical quantities in the open boundary system in terms of the closed boundary eigenstates [Mam08], but with a large number of eigenstates. It has been shown in [Mam03] that the use of the von Neumann boundary condition; i.e. zero derivative for the wavefunction at the boundary, result in a dramatic reduction in the needed number of eigenstates. Conceptually, we can choose whatever Hamiltonian for the isolated device as along as we use the appropriate self energy that compensate for the chosen Hamiltonian. The best way to mimic open boundary systems with an energy-independent eigenvalue problem is to use von Neumann boundary condition. Fig. 2.12 shows the transmission function between terminal 1 and 2 versus energy for three terminal device examined in [Mam03]. It is clear that von Neumann boundary condition allows one to use only a small percentage (less than 10%) of the isolated system eigenstates to find the open boundary system transmission function with acceptable accuracy.

Figure 2. 12: Three terminal device transmission function between terminal 1 and 2 examined in [Mam03], versus energy obtained using the CBR method. The solid line corresponds to using all of the isolated system eigenstates, dotted line for only 7% of the eigenstates using von Neumann boundary condition and dashed line for the use of 7% of the eigenstates using Dirichlet boundary condition.

2.5.2.5

Mode-Space Reduction in Single-Band Case

Numerical calculation effort of the transmission function and the spectral function of the open system can be further reduced in the single-band case. This is done by transforming the Green’s function and the self-energy into a basis of the leads mode-space [Mam05]. The lead modes are orthogonal and complete and, therefore, can serve as a basis for our matrices representation instead of the real-space basis. The key idea behind choosing these modes as a basis is that since the potential is constant

inside a given lead, then its modes are uncoupled which results in diagonal selfenergy matrices [Mam05]. Besides diagonal self-energy matrices, only few modes contribute to the transmission at a given energy [Mam05] which is physically meaningful. If the energy, at which the transmission to be found, is larger than the mode eigenenergy, then this mode will be propagating and will contribute to the current otherwise the mode will be decaying and won’t contribute to current. The self energy matrix in mode-space, taking into account von Neumann boundary condition, is given by [Mam08]:

Σ mm '

t (1  eikm a ) , m  m'   0 , m  m' 

(2-70)

2 is the coupling between 2m* a 2 the points along the propagation direction, a is grid point separation in that direction and mode dispersion relation is given by E   m  2t (1  cos(km a)) . Due to the diagonal nature of the self-energy matrix in the mode-space representation; the broadening function is given only by the imaginary part of the self-energy matrix: where m is the mode index (not grid point index), t 

Γ C  modespace  2Im( Σ C modespace )

(2-71)

This makes the broadening function matrix size N mod e  N mod e instead of N C  N C where N mod e is the number of propagating modes and N mod e  N C . In contrast to the transmission, purely real elements in the self-energy matrix are needed to determine the Green’s function of the open device G . However, for the case that both the device and the leads have the same cross section, then the lead decaying modes can be, safely, neglected without any loss of accuracy [Mam05]. The diagonal elements of the spectral function in real-space, making use of the leads mode-space reduction, can be obtained from the following expression [Kha07]:

A(i, i, E)   Greal-mode- space(i, m; E) Γ (m,m; E) 2

(2-72)

m

where Greal-mode- space(i, m; E ) is given by: Greal-mode- space(i, m; E)   m'

α

ψα (i)ψα (m' ) 1 BC (m' , m) ,  0 E  εα  iη

(2-73)

Now, the number of operations can be estimated as Nop  N E Neigen N mod es N [Kha07]. This makes a great reduction in the computational burden due to the absence of quadratic terms in N or N eigen .

CHAPTER 3

DOUBLE-GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

The term real-space (RS) means that the Schrödinger equation is solved in the coordinate space r. There is more than one approach that solves Schrödinger equation in real-space, the finite difference method (FDM) is one of them. The FDM results from a Taylor expansion of the involved wave function [Kik05]. In this chapter, the finite difference method was adopted to solve the two dimensional (2D) Schrödinger equation in real-space within the NEGF. The RS representation is the most accurate representation used in the NEGF [Ven02, Rah05 and Kha07] due to the following reasons 1) it can be used to accurately predict the electrical characteristics of double-gate (DG) MOSFETs whether the Si body is ultra-thin or not, 2) it can be used to accurately predict the electrical characteristics of electronic devices with arbitrary oriented wafer orientation, 3) non-coherent scattering can be accurately included in the simulation and finally 4) accurate calculation of the gate leakage current can be achieved using the RS representation. We are going to make use of the first and second advantages of RS in this chapter. The third and fourth advantages are left as a future work. The RS representation, however, suffers from being heavy in computation and, therefore, another representation called the uncoupled mode-space (UMS) is used instead. Unfortunately, the UMS is limited to simulate structures where the potential distribution is separable. The validity of the UMS for DG MOSFET simulation will be addressed in chapter 4. In this chapter, the computationally efficient, Recursive Green’s Function (RGF) and Contact Block Reduction (CBR) methods (previously discussed in chapter 2) are applied to DG MOSFETs and implemented in FETMOSS. Then, we use our implementation to study the accuracy of the CBR method for DG MOSFETs simulation. Simulations of DG MOSFETs with channel lengths of 15 and 5 nm shows that, the needed percentage of eigenstates (for acceptable accuracy) is bias dependent, and can vary from 6% in the on-state to 40% in the off-state. In addition to the RGF and the CBR methods, a new computationally efficient method, based on the gauss elimination, is proposed and implemented. Finally, the three methods are used to simulate a 5 nm DG MOSFET and their simulation times are compared versus the traditional NEGF with full matrix inversion.

The rest of this chapter is organized as follows: section 3.1 presents the quantum transport model used for the DG MOSFET, lists the assumptions made by the model and shows how 3D devices with one infinitely large dimension can be treated as 2D devices. Section 3.2 shows how the RS approach can be applied to the DG MOSFETs under the assumptions discussed in section 3.1. Section 3.3 deals with our implementation of the RGF and the CBR methods to DG MOSFETs and propose a new method based on gauss elimination. In section 3.4, the integration of the RS solution into the FETMOSS simulator is presented. The new transport model is then benchmarked in section 3.5. The CBR method accuracy is studied in section 3.6. Finally, a comparison of the simulation time taken by the various methods (full matrix inversion, the proposed method, the RGF and the CBR method) is given in section 3.7.

3.1 QUANTUM TRANSPORT MODEL The DG MOSFET model device geometry is shown in Figure 3.1. We make the following assumptions: Z

Top Metal Gate X Tox

Y n+ Source

Tsi

Oxide p Channel

n+ Drain

Oxide

Bottom Metal Gate

LS

L

LD

Figure 3. 1: A model double-gate MOSFET used in this work

1) The channel length in x-direction is shorter than any characteristic scattering length such that the device is operating in the ballistic limit [Lel00]. The mean free path for undoped 2.5nm thick SOI layers was estimated to be ~ 10nm [Tho04]. Thus ballistic transport assumption may be valid for Si DG MOSFETs devices with channel length ~ 10nm . 2) The width of the device in the z-direction is so large compared to other dimensions of the active device, such that the potential along that direction is rendered constant. 3) The metal contacts are so large such that thermal equilibrium is maintained and the Fermi level in these regions is determined by the applied voltage [Dat00]. The interface between the active device and these contacts are the actual sites for the energy dissipation [Muk06]. This dissipation is the

52 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

mechanism responsible for maintaining thermal equilibrium in the metal contacts 4) N-channel transistor where holes contribution, to both the transport and the electrostatic problems, can be neglected. In other words, only electrons are considered in this model. 5) No electron penetration in the insulator region. This is equivalent to say that the envelope wave function inside the insulator region is zero and that the active device is decoupled from upper and lower gate contacts. By this assumption, we implicitly set the gate current I G  0 during the selfconsistent solution. 6) An effective mass Hamiltonian is used to model the electron transport. The effective mass approximation is valid as long as the potential varies slowly relative to the atomic scale [Dat89]. 7) A single band Hamiltonian is used to model the electron transport. The validity of the single band assumption emerges from the fact that the Si conduction-band minimum occurs at k  0.82 / a  from the Brillion zone center along the 100  direction and other minima in the Si conduction-band structure occur at considerably higher energies and typically ignored [Pie03]. Generally, the 3D single band effective mass Schrödinger equation for a homogenous medium can be written as:  2   2

  1 2 1 2 1 2   *   EC x, y, z  x, y, z   Ex, y, z    2 * 2 * 2  m y y m z z   m x x 

where ( x, y, z ) is the envelope wave function, m*x , m*y and m*z are the electron effective mass in x- , y- and z-direction respectively, EC is the conduction band edge and E is the total energy. The eigenstates in the transverse (z- direction) are plane waves as a consequence of the assumption of invariant potential in that direction; i.e. Ec  Ec x, y  . The envelope wave function can be expanded in terms of the orthonormal basis  ( x, y) exp( jk z z ) / w

where the quantum number

kz

corresponds to the transverse eigenenergy Ekz   k z / 2m and w is the transistor 2

2

* z

channel width. The 2D wave function  ( x, y) is obtained from the solution of the 2D Schrödinger equation:

 2   2

  1 2 1 2   *   Ec x, y  x, y   El x, y  (3-1)  2 * 2  m y y   m x x 

where El  E  Ek z is the longitudinal energy due to motion in x- and y- direction. As long as we neglect elastic and inelastic scattering process that couple different transverse modes, we can think of the transverse modes as separate 2D

devices connected in parallel [Dat00]. Each transverse mode k z has an extra transverse energy Ek z that should be added to the longitudinal energy whenever the total energy is needed. The Fermi Dirac function f needs the total energy as an argument and, therefore, it should be replaced by another function F that takes the transverse modes into account and given by [Dam03]:  Ef  E  2m* zk BT  (3-1) F ( E , E f )   f El  Ek z , E f  1 / 2  2  kz  k BT 



where



1/ 2 is the Fermi-Dirac integral of order  1 / 2 and given by: 1 / 2 x  

1







0

t 1 / 2 dt 1  e t  x 

The summation was carried out by transforming it to integration over the transverse energy Ek z using the relation connecting it to the quantum number k z . The above function accounts for all the transverse modes exactly as well as the sum over spins. Now the device can be treated as if it were purely 2D [Dam03].

3.2 REAL-SPACE APPLIED TO DG MOSFETS In this section, direct descritization of eq. (3-1) is applied to solve the transport problem in the RS. In section 3.2.1, the RS Hamiltonian matrix is obtained. Then, section 3.2.2 presents the NEGF treatment of the problem.

3.2.1 The Real-Space Hamiltonian Matrix Using the central difference approximation for the second derivative of and denoting  ( xi , y j ) by  i, j then,

 2 i , j x

2

  i, j



2

y

2





 i 1, j  2 i , j   i 1, j

x 2

(3-2)

 i , j 1  2 i , j   i , j 1

y 2

Substituting with eq. (3-2) into eq. (3-1), we obtain :  t y i , j 1  t y i , j 1  2t y  2t x  i , j  t x i 1, j  t x i 1, j  EC i , j i , j  El i , j (3-3)

where t x 

2 2 t  and y 2 2 2m* x x  2m* y y 

Counting the 2D mesh nodes vertically as shown in Figure 3.2 (a), eq. (3-3) is rewritten as:

54 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

 t y i 1  t y i 1  2t y  2t x  i  t x i  Ny  t x i  Ny  Eci i  El i (3-4)

where the counter i runs from i=1 (top left corner node) to i  N x N y (bottom right corner node). The simulation domain is the region enclosed by the bold black line shown in Figure 3.2 (b). There are three remarks that have to be mentioned and clarified here: 1) As a result of the assumption that the insulator region is impenetrable barrier, the envelope wave function inside the insulator region is zero. 2) The nodes that belong to metal contacts are outside the simulation domain. 3) Closed boundary condition will be temporarily assumed at the source/drain contact. The correct open boundary condition will be incorporated later using the self energy concept According to the preceding remarks, node 1 is that at the top left corner of the source region and node N x N y is that at the bottom right corner of the drain region. X Y

 1

Ny  2

2 

Ny  1  Ny

 ( Nx  1) Ny  1

Ny  1



 



 







 

NxNy  1



NxNy

2 Ny  1 2 Ny

Oxide

( Nx  1) Ny  2

  



   

(a)

Source

Drain

Channel

 Oxide



(b)

Figure 3. 2 (a) Mesh nodes vertically counted using one index. (b) The simulation domain is enclosed by the bolded black line.

Writing eq. (3-4) at each node in the simulation domain, a set of linear equations are obtained and cast in the matrix form: H l ψ  ECψ  El Iψ

where  t x  1  α β 0 0   0    β α β 0  2    β   Hl   ψ     0            N N   0  0 0  β α  N x N y N x N y  x y  N x N y 1

0   t x  0          t x  N N y y 0



 ty 0  0  2t x  2t y  t 2t x  2t y  t y  0  y α  0        0   t y 2t x  2t y   0

N y N y

  EC     

       0  Ec N x N y  N N  N N x y x y

Ec1 0   0 Ec2 0  



0

0

0 0 

The Hamiltonian matrix H l has a physical meaning that worth to be mentioned. The active device can be thought as (see Figure 3.3 (a)) being composed of vertical layers adjacent to each other. Each layer represents a big site with on-site energy α and, the coupling energy between any adjacent sites is β .The matrix β is diagonal since it couples the horizontally adjacent nodes with energy  t x .The first site (left most one) has no adjacent site to the left due to the closed boundary condition applied when deriving the Hamiltonian matrix. Also the last site (right most one) has no adjacent site to the right due to the same reason. Each layer is composed of (see Figure 3.3 (b)) a set of nodes (smaller sites), each with on-site energy 2t x  2t y , and the coupling energy between any two vertically adjacent nodes is  t y . By recognition of this physical meaning of the Hamiltonian matrix, one can deduce it for other device structures without going into mathematics. Actually, the previously mentioned physical meaning is very similar to that adopted by the tight binding model (TBM) which is typically used for calculations of electronic band structure and energy gap [Gor97]. 

X





Y





2t x  2t y





   



 

 









 









 

 

  

 tx

2t x  2t y  ty



(a)

 tx



2t x  2t y  ty 2t x  2t y



2t x  2t y  ty 2t x  2t y   

 tx

 tx



2t x  2t y  ty 2t x  2t y

(b)

Figure 3. 3: (a) The active device is composed of layers adjacent to each other. Each layer has on-site energy α and, the coupling energy between any adjacent layers is β . (b) Each layer is composed of smaller sites (nodes) with on-site energy 2tx+2ty and coupling hoping energy –ty between vertically adjacent nodes and –tx between horizontally adjacent nodes.

3.2.2 The NEGF in Real-Space Closed boundary condition assumes that the active device is isolated from the outside world. Open boundary condition, in contrast, means that the active device is coupled to the drain and source contacts and in turn to the outside world. Figure 3.4 shows the overall system composed of the active device and the contacts. The Green’s function of the active device can be calculated taking into account the effect of coupling to contacts using the self energy concept [Dat00].

56 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

The retarded Green’s function of the active device is given by:

G  El I  H d   S   D 

1

(3-5)

where H d  H l  EC ;  S and  D are the source and drain contact self energy respectively. For the system shown in Figure 3.4, the self energies are given by:  βg S β 0  0   0      S           0  N N N N  0 x y x y

(3-6)

0  D     0

(3-7)

0        0    0 βg D β  N N  N N x y x y  

where g S and g D are the surface Green’s functions [Dat05] of the source and drain contacts respectively and given by [Dam03]:



α  E  i0



 α  E  i0  qVS I 

gS 



 qVS I  4t x2 I



(3-8)



(3-9)

2

2t x2

gD 



α  E  i0



 α  E  i0  qVD I 



2

 qVD I  4t x2 I

2t x2 









































    



  







 

 







  



 

 







 



   



 

 







  





 

 







  



source contact





 



   



   

      

source

channel

drain

drain contact

active device

Figure 3. 4: The overall system composed of the active device and the contacts.

Once g S and g D are obtained using (3-8) and (3-9), the self energy matrices  S and

 D are computed from (3-6) and (3-7) and the broadening functions, Γ S and Γ D are obtained as:

  i 

Γ S  i  S   S ΓD

D



 D

 

(3-10) (3-11)

The device Green’s function G is obtained using (3-5) and the spectral functions due to the source/drain contacts can be afterwards obtained as:

AS  GΓ S G 

(3-12)

AD  GΓ DG 

(3-13)

The total spectral function is the sum of the spectral functions due to the source/drain contacts: A  AS  AD The longitudinal-energy-resolved electron density at a grid point i is obtained by [Ven03]: n(i; El ) 





1 AS i, i; El F ( El , E f S )  AD i, i; El F ( El , E f D ) (3-14) 2xy

where E f S and E f D are the Fermi level of the source and drain contacts respectively. The factor

1 comes from the process of 2D discretization [Dat00]. In a discrete xy

representation the wavefunction is normalized to the number of points N x N y . The real wavefunction in a continuous representation, however, is normalized to the area of the domain N y N x  xy . The longitudinal energy resolved electron density n(i, El ) is, further, summed over the Si six conduction band valleys [Pie03] and, finally, the total electron density n(i) is obtained by integration over the longitudinal energy. The transmission coefficient from the source contact to the drain contact is defined in terms of the Green’s function and the broadening function as:





TSD  Trace Γ SGΓ DG  (3-15) The longitudinal-energy-resolved terminal current in the ballistic limit is, afterwards, obtained as: q I ( El )  TSD F ( El , E f S )  F ( El , E f D ) (3-16) 2 The terminal current is, further, summed over the six conduction band valleys and, finally, integrated over the longitudinal energy.





3.3 COMPUTATIONALLY EFFICIENT METHODS The self-consistent solution involves electron density calculation for many times during its iterations. Therefore, electron density calculation has a great impact on the total simulation time and must be as efficient as possible. However, the Green’s function (needed for the spectral function calculation) is calculated using an inversion of a huge matrix as given by eq. (3.5). In chapter 2, the RGF and the CBR method has been reviewed and their efficient computations were emphasized. The application of both of them to the DG MOSFET and their implementation into

58 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

FETMOSS will be detailed hereinafter. We are going also to propose another simple and efficient way to calculate the spectral function and, hence, the electron density using Gauss elimination method.

3.3.1 The Recursive Green’s Function Algorithm Eq. (3-5) can be rewritten as: GD  I

(3-17)

where D  El I  H d   S   D  . Then, the blocks of the matrix D , needed for the RGF algorithm are given by:

D1,1

E  C1  El I  α  βg s β    0 

DN x , N x

Dq , q



 EC  N x N y  N y 1  El I  α  βg D β   0  0 

 EC  qN y  N y 1  El I  α    0 

0    EC N  y  0

0



0

0

EC N

xN y

(3-18)

    

(3-19)

0   , q  2,3,4, N  1 (3-20) x    ECqN  y 

Dq 1,q  Dq,q 1  β , q  1,2,3, N x  1

(3-21)

The in-scattering self-energy is given by:







Σ in El   Γ S El F El , E f S  Γ D El F El , E f D







where Γ S  i Σ S  Σ S and Γ D  i Σ D  Σ D



(3-22)



To evaluate eq. (2-47) , we need only the diagonal blocks of the in-scattering selfenergy matrix. The diagonal blocks are given by:



in Σ11  i[ β(gS - gS ) β ]F El , E f S





Σ Ninx N x  i [ β ( gD - gD ) β ] F El , E f D

(3-23)



(3-24)

inqq  0 , q  2,3,4, N x  1 (3-25) Once the quantities in eq. (3-18) to eq. (3-25) are ready, the RGF algorithm is applied in the sequence summarized in chapter 2. Then longitudinal-energy-resolved electron

density at a grid point i is obtained from the diagonal elements of the correlation function:

n(i; El ) 

1 G n i, i; El  2xy

(3-26)

Using eq. (2-54), the longitudinal energy resolved terminal current is obtained as: I





q Trace [G1,1 ( E )  G1,1 ( E )]F ( El , E fs )  iG1n,1 ( E ) β(g S - gS ) β 2



(3-27)

n Note that G n (1,1; E ) isn’t equivalent to G1,1 (E ) as the first one denotes the first element of the matrix while the latter one denotes the first block of the matrix which corresponds to layer 1 in the device

3.3.2 The Contact Block Reduction Method The CBR method splits the simulation domain into two sub-domains, D: the interior part of the device and C: the boundary region that connects the interior parts to the contacts. For this reason, the grid points that connect the active device to the contacts are numbered first then the interior device points are numbered as shown in Figure 3.5. X Y

1

2 Ny  1



( Nx  1) Ny  1

Ny  1

2

2 Ny  2



( Nx  1) Ny  2

Ny  2

Ny  1

 

 

 

 



 











 

NxNy  1



NxNy



Ny

3Ny  1

3Ny

1st

3rd



 2 Ny  1

2 Ny

last

2nd

Figure 3. 5: Grid points that connect the active device to the contacts are numbered first then the interior device points are numbered.

Note that the 1st layer is the left most layer, the 2nd layer is the right most layer, and the 3rd layer is the 2nd one on the left. As a result of this numbering scheme we have 1) the 1st and the 3rd layers are coupled together and 2) the 2nd and the last layers are coupled together. Thus the isolated device Hamiltonian matrix, with Dirichlet boundary condition, will be no more tri-diagonal matrix as given in the matrix from of eq. (3-4), there will be some modifications in the 1st, 2nd, 3rd and last rows of the matrix. The new Hamiltonian matrix is given by:

60 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

α 0 β 0  0 α 0    β 0 α β 0  0 0 β α β Hl    0 0 β                   0 β 0   

   0    0 

       β 0

      β α β

0 β  0  0    0  β α 

(3-28)

The self-energy matrix accounts for the coupling to the source and drain contacts. Now, the source contact is coupled to the 1st layer while the drain contact is coupled to the 2nd layer. As a result, the self-energy matrix should be modified and is given by:   S     0 

0   D  N  N C C 0

 0  0  N N

(3-29)

where N C  2 N y , N  N x N y , Σ S  βg S β , Σ D  βg D β , g S and g D are the surface Green’s function of the source and drain contact given by eq. (3-8) and eq.(3-9) respectively. Applying von Neumann boundary condition, instead of Dirichlet boundary condition, needs further modifications in the Hamiltonian and self-energy matrices in the blocks related to the 1st and the last layers. The block matrix β should be added to the marked blocks of the Hamiltonian matrix as given by (3-30):

0 α  β  0 α β   β 0  0  0 Hl    0             0 β 

β 0 α β 0    0

0  β α     

  0 β     

   0    0 

       β 0

       α β

0 β  0  0    0  β α 

(3-30)

The self-energy matrix must be also modified such that the overall problem isn’t changed; i.e. we must subtract the terms added to Hamiltonian matrix. Therefore, we need to subtract the block β which leads to the following expressions for the source and drain self energy: Σ S  βg S β - β

(3-31)

Σ D  βg D β - β

(3-32)

The eigenstates of the isolated device are obtained from the eigenvalue problem:

Hl  EC ψ   ψ

(3-33)

Then, the isolated device Green’s function matrix in the contact region is obtained by: GC0 El    

ψ (1  N C )ψ (1  N C ) El     i

(3-34)

Now to make use of the mode-space reduction technique, the propagating modes of the source and drain contacts must be obtained. The source and drain contacts’ modes can be obtained from the eigenvalue problems: H S ψ S m   S mψ S m

(3-35)

H Dψ Dm   Dmψ Dm

(3-36)

where  tx 0  0  2t x  EC1  t 2 t  E  t 0   x x C2 x     0  HS   0   0  t x 2t x  ECN y 1  tx     0  0  tx 2t x  ECNy   N y N y

 tx 0  0 2t x  EC Ny 1   t  2 t  E  t 0  x x C Ny  2 x    0    0 HD      0  t 2 t  E  t x x C 2 N y 1 x     0  0  t 2 t  E x x C 2 Ny   N y N y

The propagating modes satisfy the inequality    El  max . These modes are used to transform the real-space matrices to their mode-space version. The transformation matrix is given by:





 ψ S 1 ψ S 2   ψ S  N mode Rr tom   0 



ψ

0

D 1  ψ D  2 

 (3-37)  ψ D  N mode 





Mode-space self-energy can be obtained by transforming eq. (3-31) and (3-32) to their mode-space version. Alternatively, the self-energy matrix can be calculated in modespace, directly, with von Neumann boundary condition according to eq. (2-68) as: ΣC mod espace  Σ S mod espace  Σ Dmod espace (3-38) where

62 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

t x  t x eik S 1x 0   ik S 2 x 0 t  t e  x x  Σ S  mod e  space       0  0   0 0 0  t x  t x eik D1x 0  0 t x  t x eik D 2 x Σ D  mod e  space  0     0  

  0  0  ik SN mode x t x  t xe  0     0       0 ik DN mode x  0 tx  txe  0 

The broadening functions, Γ S mod espace and Γ Dmod espace are calculated using eq. (269). The source- and drain-contact spectral functions are then calculated using eq. (270). The longitudinal energy resolved electron density is, afterwards, calculated by eq. (3-14). Finally, the transmission function is given by:





TSD  Trace Γ S -mode- spaceGC -mode- spaceΓ D-mode- spaceGC -mode- space



(3-39)

The longitudinal-energy-resolved terminal current in the ballistic limit is, afterwards, obtained by eq. (3-16).

3.3.3 Proposed Method Based on Gauss Elimination The idea in the proposed method is based on the sparse nature of the broadening function in the coherent transport case. This sparse nature can be seen from eq. (3.6), (3.7), (3.10) and (3.11). The broadening functions are given by:





 β( g S - g S ) β N  N 0  0  y y   0     ΓS  i        0   0  N N  N N  x y x y 0  Γ D  i   0 

       0

     0   β( g D - g D ) β N y  N y   N x N y N x N y 0





Due to sparse broadening functions, we don’t need the entire Green’s function to calculate the spectral functions in eq. (3-12) and (3-13). Schematically eq. (3-12) and (3-13) can be seen as:

  (3  12)        (3  13)     

                     

           

           

           

where the shaded area is the non-zero elements needed in (or resulted from) the multiplication process. As a result, we need only the following parts the Green’s function: G (1,2)  G (1,1)  G (2,1) G (2,2) GS       G ( N x N y ,1) G ( N x N y ,2)

G (1, N y )   G (2, N y )       G ( N x N y , N y ) N x N y N y 

G (1, ( N x  1) N y  2)  G (1, ( N x  1) N y  1)  G (2, ( N  1) N  1) G (2, ( N x  1) N y  2) x y GD       G ( N x N y , ( N x  1) N y  1) G ( N x N y , ( N x  1) N y  2)

(3-40)

G (1, N y )   G (2, N y )  (3-41)      G ( N x N y , N y ) N x N y N y 

This can be understood from the physical meaning of the Green’s function; it is the impulse response of the system. More precisely, G (i, j ) is the response at point i due to an impulse excitation at point j . In our case, excitation sources are the source and drain contacts. Thus we have excitations only at (see Figure 3.2(a)) points 1, 2… N y and at points ( N x  1) N y  1 , ( N x  1) N y  2 … N x N y . The matrices in eq. (3-40) and (3-41) can be obtained using Gauss elimination method by the following equations: GS  El I  H d  Σ S  Σ D  \ I S GD  El I  H d  Σ S  Σ D \ I D

(3-42) (3-43)

where the A \ B denotes division of the B by A using Gauss elimination method, I S and I D are given by:

64 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

1 0  0 IS      0

0



1



0 

 









0 0 0    0 1 ID    1      0 N N  N 0 x y y

0     0   1 N N N x y y

    0    

1

 

It is important to note that, using eq. (3-40) and eq. (3-41), we calculate only N x N y  2 N y elements of the Green’s function instead of calculating N x N y  N x N y elements. This is a great reduction in simulation time from the operations counts point of view. The spectral functions are calculated using the following equations

AS  GS [ β( gs - gS ) β ]GS

(3-44)

AD  GD [ β( gD - gD ) β ]GD

(3-45)

It is interesting to note that transmission function can also be efficiently calculated without the full Green’s function matrix. Schematically eq. (3-15) can be seen as:    (3  15)  Trace          Trace     

           

     

      

     

       Trace         

      

where the shaded area is the non-zero elements needed in (or resulted from) the multiplication process and the dotted area the elements that, after multiplication, results in the diagonal elements needed to perform the trace operation. Then, the transmission can be efficiently calculated using the following equation:



TSD  Trace [ β( g s - gS ) β ]GDS [ β( g D - g D ) β ]GDS





(3-27)

where GDS is a subset of G D and given by

GDS

G (1, ( N x  1) N y  2)  G (1, ( N x  1) N y  1)  G (2, ( N  1) N  1) G (2, ( N  1) N  2) x y x y      G ( N y , ( N x  1) N y  1) G ( N y , ( N x  1) N y  2)

G (1, N y )   G (2, N y )       G ( N y , N y ) N y N y 

(3-28)

3.4 INTEGRATION INTO FETMOSS FETMOSS is a software tool for the 2D simulation of DG MOSFETs [Abd06]. The tool works under MATLAB environment. In FETMOSS, Poisson’s equation is solved using the finite elements method (FEM) while the UMS NEGF approach is used to simulate the ballistic quantum transport. The confined modes are obtained using the transfer matrix method (TMM) [Abd04]. The use of the UMS makes FETMOSS computationally efficient and can be used for extensive device simulation and design. However, the UMS imposes some restrictions on the simulation capabilities of FETMOSS: 1) the simulated DG MOSFET must be of ultrathin body, 2) closed boundary condition must be applied at either Si/SiO2 interface or SiO2/metal contact interface to be able obtain the confined modes. This prevents the accurate calculation of the gate leakage current and finally 3) scattering can’t be included accurately in the simulation since the scattering potential may not be satisfying the UMS assumption. Moreover, FETMOSS is capable of simulation only DG MOSFET with (001) oriented wafers.

The RS NEGF model discussed previously has been integrated into FETMOSS. By this integration, the 1st restriction has been overcome automatically. The accurate inclusion of gate leakage current or scattering mechanisms becomes possible but left as a future work. Moreover, accurate simulation of (001), (110) and (111) oriented wafers become possible. The difference between these wafer orientations will be mentioned later.

The flow chart of the program is given in Figure 3.6. The user specifies whether to use UMS or RS as a NEGF transport model. Other inputs specified by the user are device dimensions, doping, wafer orientation and the externally applied bias. The simulation starts by assuming an initial guess for the potential distribution in the device. According to this potential, the NEGF is used to calculate the electron concentration in the device. With the electron concentration is known, Poisson’s equation can be solved yielding a new potential distribution. The new potential is compared to the old potential and the solution cycle is repeated until self consistent solution for the potential is obtained (until the difference in potential between two successive iterations is below a certain tolerance.)

66 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Specify NEGF transport model (whether RS or UMS), geometry, doping and bias

Assume initial guess for the potential EC

NEGF ECold  n

ECold  ECnew

Poisson’s equation

n  ECnew No

ECnew  ECold   Yes NEGF Calc. terminal current and other outputs

Figure 3. 6: The flow chart of FETMOSS.

3.4.1 Initial Guess for Real-Space A good initial guess has a great impact on the number of iterations of the selfconsistent solution loop. The UMS approach is less accurate than the RS approach but more efficient from the number-of-computations point of view. This fact motivated us to use the UMS results as initial guess for the RS. If the simulated DG MOSFET body is thin enough, the RS approach may add only one iteration to the UMS iterations, to achieve a self consistent solution. If the simulated DG MOSFET body is relatively thick, the RS approach will add a few iterations to achieve self-consistency. Even in this case, we found that taking the UMS results as initial guess for RS also reduces the number of iterations done by the RS and consequently, the overall simulation time needed.

3.4.2 Numerical Integration It is also important to note that a good choice of the energy grid strongly affect the convergence of the self-consistent scheme. The electron density is obtained by integrating eq. (3-14). The Numerical integration is done by summing the values for each energy step multiplied by the energy grid spacing ΔE. Using a uniform energy grid makes the convergence of the self-consistent loop very poor since peaks in the density of states shift from one iteration to another. This shift makes the error due to the numerical integration varies from one iteration to another. Therefore, FETMOSS uses a non-uniform energy grid and the intelligent integration algorithm of the “quad” MATLAB function [Mat02] to keep the error due to numerical integration constant over all the self-consistent iterations.

3.4.3 Silicon Anisotropy Figure 3.7 displays the constant-energy surfaces characterizing the conduction-band structure near the conduction band edge of Si. The Si conduction-

band minimum occurs along six equivalent directions 100 and accordingly we have six valleys. The effective mass in Si has an anisotropic nature and; the constant energy surface for each valley is ellipsoid with longitudinal and transverse effective masses [Pie03]: ml  0.91m0 and mt  0.19m0 . Z

010

X

5 Y

4 1

100

2 3

6

001 Figure 3. 7: Constant-energy surfaces characterizing the conduction-band structure near the conduction band edge in Si.

The original FETMOSS assumes the transistor shown in Figure 3.1 is fabricated on wafer with (001) surface orientation; the x-, y- and z-directions are aligned with 100 , 001 and 010 respectively. We’ve generalized FETMOSS to simulate wafers with (001), (111) and (110) surface orientations. The effective masses for each wafer orientation are given in Table 3.1 [Rah05]. Table 3-1: Transport, width, and confinement effective masses and degeneracies for three different technologically important semiconductor wafer orientations [Rah05]. Wafer

(001)

Length

[100]

Width

[010]

m *x

m *y

m *z

Degeneracy

ml

mt

mt

2

mt

mt

ml

2

mt

ml

mt

2

(2ml  mt ) / 3

3ml mt

mt

2

4

/( 2ml  mt ) (111)

 

[2 11]



[0 11]

(2 / 3)mt

3ml mt /

(ml  mt )

[(2ml  mt )

(2ml  mt )

/2

2ml mt

(ml  mt )

/( ml  mt )

/2

mt

mt

/( ml  mt )]

mt (110)

[001]

 

[0 1 0]

ml

4

2

68 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

For (001) wafer orientation, the six valleys are reduced to three with a degeneracy factor of 2. Moreover, Valleys can be categorized as primed or unprimed [Lun06]. Unprimed valleys are those with heavy effective mass in the gate confinement direction which are valley 1 and 2 for the case in hand. Heavier effective mass means smaller quantization energy and larger occupation probability. Therefore, for ultrathin bodies only unprimed valleys can be taken into consideration without perceptible reduction in accuracy. All valleys, however, must be considered when simulating thick bodies. For (111) and (110) orientations, the six valleys are reduced to two with a degeneracy factors of 2 and 4. The two valleys of the (111) orientations have the same effective masses and therefore, none of them dominate the other. The (110) orientation, however, has one valley with heavier effective mass (a factor of 1.5) than the other. This factor isn’t large enough to make one valley dominate the other one. Therefore, only (001) orientation has primed and unprimed valleys.

3.5 BENCHMARKING FETMOSS FETMOSS has been benchmarked using an online available simulator NanoMOS 3.0 [Ren03] on nanohub: www.nanohub.org. NanoMOS3.0, is a 2-D simulator for thin body (less than 5 nm), fully depleted, DG n-MOSFETs. A choice of five transport models is available (drift-diffusion, classical ballistic, energy transport, quantum ballistic, and quantum diffusive). The transport models treat quantum effects in the confinement direction exactly and the names indicate the technique used to account for carrier transport along the channel. Each of these transport models is solved self-consistently with Poisson's equation. All the models in NanoMOS are based on the UMS and that’s why accuracy isn’t guaranteed for body thickness larger than 5 nm. Moreover, nanoMOS3.0 is capable of simulating (001) wafer orientation only. The quantum ballistic transport model that uses UMS NEGF has been used to bench mark the quantum ballistic transport model that uses RS NEGF in FETMOSS. It is expected to find well agreement, between UMS NEGF in NanoMOS and RS NEGF in FETMOSS, for ultra-thin body thickness where the UMS assumptions are valid. Two samples of nanoscale DG MOSFETs have been simulated using NanoMOS 3.0 and FETMOSS. The simulated devices dimensions, doping concentration, material parameters, simulator options, and FD grid and supply voltage are given in Table 3-2. The simulated wafer orientation is (001). The simulation results are shown in Figures 3.8, 3.9, 3.10 and 3.11.

Table 3-2: The simulated devices dimensions, doping concentration, material parameters, simulator options and FD grid. Parameter

Dimensions

Doping

Material

1

2

channel length L

5 nm

15 nm

source and drain length LS , LD

5 nm

5 nm

oxide thickness Tox

1 nm

1.2 nm

silicon (body) thickness TSi

2 nm

4.6 nm

channel doping

1010 cm-3

0 cm-3

source and drain doping

2x1020 cm-3

1020 cm-3

junction doping profile

step

step

silicon relative permittivity oxide relative permittivity

Options

Grid

Supply Voltage

Device

 Si

 ox

11.7 3.9

0

0

11.7 3.9

0

0

top and bottom gate work function  m

4.5 eV

longitudinal electron effective mass ml*

0.91

m0

0.91

m0

transverse electron effective mass mt*

0.19

m0

0.19

m0

oxide penetration flag

no

No

valleys

unprimed

All

self-consistence tolerance 

10-3 V

10-3 V

Poisson’s tolerance

10-6 V

10-6 V

vertical node spacing y

0.1 nm

0.2 nm

horizontal node spacing x

0.2 nm

0.2 nm

VDD

0.7 V

0.9 V

4.4 eV

IDS (μA/μm)

IDS (μA/μm)

70 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

VGS (V)

IDS (μA/μm)

IDS (μA/μm)

Figure 3. 8: The IDS-VGS C/CS of device 1 at VDS=25mV. The left axis is the log-scale while the right axis is the linear-scale.

VGS (V)

Figure 3. 9: The IDS-VGS C/CS of device 1 at VDS=0.7V. The left axis is the log-scale while the right axis is the linear-scale.

IDS (μA/μm)

IDS (μA/μm)

VGS (V)

IDS (μA/μm)

IDS (μA/μm)

Figure 3. 10: The IDS-VGS C/CS of device 2 at VDS=25mV. The left axis is the log-scale while the right axis is the linear-scale.

VGS (V)

Figure 3. 11: The IDS-VGS C/CS of device 2 at VDS=0.9V. The left axis is the log-scale while the right axis is the linear-scale.

72 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

The figures show good agreement between the two simulators. Device 2 (Figure 3.10 and 3.11) shows less agreement than device 1 (Figure 3.8 and 3.9). This is because of the thicker body of device 2 where the UMS, used by NanoMOS, accuracy is reduced. Detailed study of the UMS validity will be given in chapter 4.

3.6 NUMBER OF EIGENSTATES IN THE CBR METHOD The CBR method doesn’t provide a clear way to determine the number of eigenstates needed for acceptable accuracy in terminal current calculations. It is expected that this number of states may vary from one structure to another or even from one bias point to another for the same device structure. In this section we are going to study the number of eigenstates needed for acceptable accuracy for the DG MOSFETs device. This study is carried out for the two most important states of the device, the on- and off-states. Seeking for the minimum number of eigenstates needed in order to have an acceptable accuracy, the simulation is repeated many times with increasing the number of eigenstates in each time. An x% of eigenstates indicates that, x N eigen  N T where N T is the total number of eigenstates which is equal to the 100 number of FD grid points. The eigenstates percentage was chosen to start from 1% and increased up to 10 %. For each time, the max energy ( Eeigenmax ) corresponding to the eigenstates percentage x%, the number of self-consistent iterations ( N iter ) and the drain current ( I D ) are recorded. The percentage difference in drain current between the given percentage of eigenstates and the 40% of eigenstates is calculated such that  (%) given by:

 (%) 

I Neigen%  I 40% I 40%

 100

The 40% eigenstates was taken as a reference, because it was found that it gives almost the same results as the 100%. The simulated DG MOSFET devices are the same as those given in Table 3-2. The initial guess was taken to be the UMS simulation results. Table 3-3 gives the results obtained for device 1 in the on-state. It is observed that starting from a 2% eigenstates, the error percentage in current calculation has dropped below a 5%. This observation indicates that taking only 2% of the eigenstates yields an acceptable accuracy. The number of iterations, however, has dropped to a single iteration staring from an 8% of the eigenstates. This indicates that taking only an 8% of the eigenstates results in a good accuracy for the charge calculation.

Table 3-3: CBR method results for device 1 in the on-state

N eigen

E eigen max

N iter

I (mA/μm)

1%

0.149

5

2.5287

6232

2%

0.7193

4

2.9415

9348

3%

1.2933

3

2.9753

12464

4%

1.8153

2

3.0164

15580

5%

2.3102

2

3.0327

18696

6%

2.8497

2

3.0428

21812

7%

3.3727

2

3.0629

24928

8%

3.8226

1

3.0741

28044

9%

4.3533

1

3.0591

31160

10 %

4.8486

1

3.0599

124640

40 %

15.4653

1

3.0571

17.3% 3.8% 2.7% 1.3% 0.8% 0.5% 0.19% 0.56% 0.07% 0.09 --

n (m-3)

3116

| I (%) |

Figure 3. 12 : The electron density along the x-direction at midpoint along the y-direction for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation.

Figure 3.12 depicts the electron density along the x-direction at midpoint along the y-direction for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation. There are strange oscillations for 1% of

74 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

eigenstates. These oscillations are suppressed by increasing the eigenstates percentage. The 8% eigenstates electron density is very close to that of the 40%. That’s why the number of iterations remains fixed starting from 8% of eigenstates. It is important to analyze the transmission function and drain current versus energy and the electron density along the device. This analysis will help to understand the effect of increasing the number of eigenstates. Figure 3.13 depicts the transmission function versus energy for 1%, 2%, 8% and 40% of eigenstates.

Figure 3. 13: The transmission function versus energy for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation.

The transmission has a zero value below -0.14 eV which indicates that the first lead propagating mode has a cut off energy of about -0.14 eV. For a 1% of eigenstates, the transmission decays incorrectly above 0.05 eV. Moreover, it has strange oscillations between -0.1 eV and 0.1eV. The 2% eigenstates has a much better transmission where the oscillations have been fairly suppressed and the decays starts at about 0.25 eV. The 8 % eigenstates transmission is very close to the 40% case. In order to appreciate the effect of inaccurate calculation of the transmission function on the current, Figure 3.14 depicts the energy spectrum of the current. It is interesting to observe that the higher energy spectrum of the transmission has almost no effect on the current. As given by eq. (3.16), the current depends on two factors 1) the transmission function and 2) the difference between Fermi Dirac integrals of order -1/2 evaluated using the source and drain Fermi-levels. The energy spectrum of the drain current peak occurs when both factors have a considerable value. For smaller

ID (µA/µm.eV)

energies than the peak range, the transmission function is very small, and for larger energies the Fermi Dirac integral is very small. Therefore, the current value is dominated by the low non-zero energy spectrum of the transmission, which is fairly accurate for eigenstates expect the 1%. That’s why only a 2% of eigenstates is enough to calculate the current with an error percentage less than 5%.

Figure 3. 14: The energy spectrum of the current for device 1 in the on-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation.

Now, let’s turn to the off-state. Table 3-4 gives the results obtained for device 1 (parameters and dimensions are given in Table 3-1) in the off-state. We can observe that only 4 % of eigenstates is adequate to calculate the charge accurately. Figure 3.15 depicts the electron density along the x-direction at midpoint along the y-direction for device 1 in the off-state for 1, 2, 8 and 40 % of the eigenstates considered in the simulation. There electron density, in the channel, is overestimated for 1% of eigenstates. The electron density is decreased by increasing the eigenstates percentage. The 4% eigenstates electron density is close from the 20% and 40%. That’s why the number of iterations remains fixed starting from 4% of eigenstates.

76 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Table 3-4: CBR method results for device 1 in the off-state.

N eigen

E eigen max

N iter

I (μA/μm)

| I (%) |

3116

1%

0.2804

3

1.8189

246

6232

2%

0.8596

2

0.5899

12.2

9348

3%

1.4824

2

1.1721

123

12464

4%

2.0162

2

0.9593

82.5

15580

5%

2.5165

1

0.8346

58.8

18696

6%

3.0355

1

0.7551

43.6

21812

7%

3.5642

1

0.4326

17.7

24928

8%

4.0095

1

0.4462

15.1

28044

9%

4.5148

1

0.6646

26.4

31160

10 %

5.0025

1

0.4584

12.8

62320

20 %

9.3549

1

0.5036

4.2

124640

40 %

15.6291

1

0.5257

--

The error percentage in current calculations has dropped below 5 % using a 20% of the eigenstates. This is a surprising result in comparison with that in Table 3.3 for the on-state. The on-state current could be calculated using only 2% of eigenstates with acceptable accuracy where the off-states current requires a much larger percentage of eigenstates. To understand why the off-state’s current needs such a large percentage, Figure 5.16 depicts the transmission function, on linear scale, versus energy for 1%, 4%, 20% and 40% of eigenstates. The transmission function in off-state seems to be similar to that in on-state except for an upward shift with energy and being smoother. This upward shift is due to the larger potential barrier in off-state than that in on-states. The current, however, needs a larger percentage of eigenstates for accurate calculation. Figure 3.17 depicts the energy spectrum of the current for device 1 in the off-state for 1, 4, 10, 20 and 40 % of the eigenstates considered in the simulation. For a 1% of eigenstates, the current peaks in-between -0.1 and 0 eV which corresponds to an extremely small transmission in Figure 3.16. That is to say the current peaks in the tunneling energy range. This peak is greatly reduced with increasing the eigenstates percentage, which indicates the CBR method overestimate the tunneling current when the eigenstates percentage is small. Tunneling current overestimation means the transmission function is overestimated in the tunneling regime.

n (m-3) Figure 3. 15: The electron density along the x-direction at midpoint along the y-direction for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation.

Figure 3. 16: The transmission function, on linear scale, versus energy for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation.

ID (µA/µm.eV)

78 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Figure 3. 17: The energy spectrum of the current for device 1 in the off-state for 1, 4, 10, 20 and 40 % of the eigenstates considered in the simulation.

Figure 3.18 depicts the transmission function, on log scale, versus energy for 1%, 4%, 20% and 40% of eigenstates. The transmission is incorrectly large between -0.1 and 0 eV for small percentage of eigenstates which results in the large error in current calculation in the off-states. In the on-state, however, the current is dominated by thermoionic emission over the barrier and, therefore accurately calculated using a small percentage of eigenstates.

Figure 3. 18: The transmission function, on log scale, versus energy for device 1 in the off-state for 1, 4, 20 and 40 % of the eigenstates considered in the simulation.

All of the above analysis have been repeated for device 2 (parameters and dimensions are given in Table 3.1). Table 3.5 gives the results obtained for device 2 in the on-state. Convergence couldn’t be reached 1% and 2% of eigenstates. The electron density, the transmission function and the energy spectrum of the current are given in Figure 3.19, 3.20 and 3.21 respectively for 3, 5, 6 and 40 % of the eigenstates. Table 3-5: CBR method results for device 2 in the on-state

N eigen

E eigen max

N iter

I (mA/um)

| I (%) |

1%

-0.3733

--

--

--

5858

2%

-0.0636

--

--

--

8787

3%

0.1883

6

3.6798

22.8%

11716

4%

0.3425

8

4.1394

13.1%

14645

5%

0.5192

7

4.3295

9.1%

17574

6%

0.6867

7

4.5587

4.3%

20503

7%

0.8293

7

4.6099

3.2%

23432

8%

0.9860

7

4.6432

2.5%

26361

9%

1.1542

7

4.6781

1.8%

29290

10 %

1.3038

7

4.6821

1.7%

117160

40 %

5.3862

7

4.7667

--

n (m-3)

2929

Figure 3. 19: The electron density along the x-direction at midpoint along the y-direction for device 2 in the on-state for 3, 5, 6 and 40 % of the eigenstates considered in the simulation.

80 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

ID (µA/µm.eV)

Figure 3. 20: The transmission function, on linear scale, versus energy for device 2 in the off-state for 3, 4, 6 and 40 % of the eigenstates considered in the simulation.

Figure 3. 21: The energy spectrum of the current for device 2 in the on-state for 3, 5, 6, and 40 % of the eigenstates considered in the simulation.

Due to the relatively large body thickness of device 2, the UMS results aren’t very accurate and the CBR method adds several iterations starting from the UMS results as an initial guess. Moreover, several modes belonging to several valleys share in the overall electron density, transmission function and current. This can be observed from the transmission function in Figure 3.20 where the transmission reaches a value of about 6 at 0.4 eV. It also can be observed from the energy spectrum of the current in Figure 3.21 where it has more than one peak. Table 3-6 gives the results obtained for device 2 in the off-state. The 50% eigenstates was taken as a reference because it was found that it gives almost the same results as the 100%. The electron density, the transmission function on linear scale, the transmission function on log scale and the energy spectrum of the current are given in Figure 3.22, 3.23, 3.24 and 3.25 respectively for 3, 20, 40 and 50 % of the eigenstates.

Table 3- 6: CBR method results for device 2 in the off-state.

N eigen

E eigen max

N iter

I (uA/um)

| I (%) |

2929

1%

-0.3260

--

--

--

5858

2%

0.0603

--

--

--

8787

3%

0.2940

3

0.5310

703.5

11716

4%

0.4909

4

0.4328

554.9

14645

5%

0.6252

4

0.3520

432.5

17574

6%

0.7692

3

0.2525

282

20503

7%

0.9304

3

0.1905

188.3

23432

8%

1.0887

3

0.1654

150.3

26361

9%

1.2111

3

0.1454

120

29290

10 %

1.3724

3

0.1306

97.6

58580

20 %

2.5247

2

0.0827

25.1

87870

30 %

3.7379

2

0.0724

9.53

117160

40 %

5.4590

2

0.0675

2.1

146450

50 %

7.4019

2

0.0661

--

n (m-3)

82 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Figure 3. 22: The electron density along the x-direction at midpoint along the y-direction for device 2 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation.

Figure 3. 23: The transmission function, on linear scale, versus energy for device 2 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation.

ID (µA/µm.eV)

Figure 3. 24: The transmission function, on log scale, versus energy for device 1 in the off-state for 3, 20, 40 and 50 % of the eigenstates considered in the simulation.

Figure 3. 25: The energy spectrum of the current for device 2 in the off-state for 3, 20, 40, and 50 % of the eigenstates considered in the simulation.

84 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Table 3-7 summarizes the results obtained for both the sample devices we studied with either on- or off-states. Let a small error in charge calculation using a small number of eigenstates means one self-consistent iteration difference with that using large number of eigenstates. Then, using Table 3-3, 3-4, 3-5 and 3-6, one can deduce the N eigench arg e given in Table 3-7. Thus, choosing a 6 % of eigenstates insures charge calculation with acceptable accuracy for either on- or off-states. Also let a small error in current calculation using a small number of eigenstates means less than 5% difference with the current of large enough number of eigenstates. Then, one can deduce the N eigencurrent given in Table 3-7. It is observed that current calculation in the on-state requires a small number of eigenstates in comparison with the off-state. A 6 % of eigenstates insures current calculation with acceptable accuracy for the onstate while a 40% is needed for the off-state. Therefore, it is suggested to increase the current calculation eigenstates with decreasing the gate voltage. Table 3-7: CBR method results summary. Device

On-state

Off-state

N eigench arg e

N eigencurrent

N eigench arg e

N eigencurrent

1

4%

2%

4%

20 %

2

5%

6%

6%

40 %

It is important to mention that our analysis was based on the 2 sample devices given in Table 3-2 This may not be accurate enough to deduce a general rule to be applied on any simulated dimensions and parameters. Our point of view was that the chosen devices were of typical dimensions; that are expected to be used in the next technology generations [ITRS]. A more general analysis covering wider range of dimensions and parameters was left as a future work.

3.7 SIMULATION TIME COMPARISON Now we have the different methods discussed throughout this chapter implemented and integrated into FETMOSS. It is the point now to compare their computational efficiency. For this purpose, the four methods (full matrix inversion, proposed method, RGF and the CBR) were used to simulate device 1. The drain voltage was kept constant at 0.7 V and the gate voltage was swept from 0.0 to 0.7 V with a step of 0.1 V. Thus, we have eight bias points. For the first bias point; i.e. VGS=0.0 V, the initial guess was taken to be the zero potential at various grid points in the device. The initial guess for any other bias point was taken from the solution of the preceding bias point, for example initial guess for VGS=0.1 V was taken from the solution of VGS=0.0 V. Figure 3.26, 3.27, 3.28 and 3.29 depicts the self-consistent error versus time for full matrix inversion, proposed method, RGF and the CBR respectively. A solution is found when the error drops below 1 mV. Once this criterion is met, the terminal current is calculated and a new bias point is initiated. This causes the error to jump to

Iteration error (V)

larger value, and the error starts decreasing again with iterations until the solution of the new bias point is found. The cycle was repeated until the eight bias points were completed.

Iteration error (V)

Figure 3. 26: The self-consistent error versus time using full matrix inversion of the Green’s function in the real-space.

Figure 3. 27: The self-consistent error versus time using proposed method for the Green’s function in the real-space.

Iteration error (V)

86 CHAPTER 3: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN REAL-SPACE

Iteration error (V)

Figure 3. 28: The self-consistent error versus time using the RGF in the real-space.

Figure 3. 29: The self-consistent error versus time using the CBR method.

It is important to observe that all the figures have almost the same values of error at various iterations. This reflects the fact that, all the methods have almost the same results due to the assumptions mentioned in section 3.1. The simulation time, however, differs considerably from one method to another. Full matrix inversion has the greatest simulation time while the CBR method has the smallest one. A summary of the total simulation time (ttotal), average simulation time per bias point (tbias) and the average simulation time per iteration (titration) is presented in Table 3.8. Table 3-8: Simulation time summary. These simulations were carried out on a home PC: Intel® Pentium 4 CPU 2.4GHz, 768 MB RAM. Method

ttotal (hr)

tbias (hr)

titeration (hr)

Full matrix inversion

75.7218

9.4652

1.2022

Proposed method

11.6525

1.4566

0.1850

RGF

5.0872

0.6359

0.0807

CBR

0.8661

0.1108

0.0135

Using the RGF or the proposed method introduces about one order of magnitude reduction in simulation time below that when using full matrix inversion. The CBR method yields the smallest simulation time with about two orders of magnitude reduction. In fact, the CBR method makes it practical to simulate and design a DG device using real-space simulation. The main disadvantage of the CBR method is that it trades off the accuracy with the simulation time. This doesn’t occur in either the RGF or the proposed method. As indicated earlier in chapter 2, the RFG is not capable of simulating devices with more than two metal contacts. This restriction doesn’t exist in either the CBR or the proposed method.

CHAPTER 4

DOUBLE-GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Computational efficiency is needed to make the self-consistent approach suitable for device design and characteristic prediction. The NEGF method has the advantage of being rigorous but the disadvantage of being heavy in computations [Ren03]. Real-space (RS) representation is the most accurate, yet complex, representation used in the NEGF [Ven02]. For nanoscale transistors, it is more convenient, to solve Schrödinger equation in the mode-space (MS) [Ren0, Wan05] where computational burden is affordable. The uncoupled mode-space (UMS) provides enormous saving in the computational burden but suffers from being valid only for ultra-thin body double-gate (DG) MOSFET [Ren03]. In contrast, the coupled mode-space (CMS) can be used for either thin or thick bodies but it is more computationally demanding than the UMS. In this chapter we review the MS approach and examine the validity of the UMS. Moreover, the possibility to reduce the CMS computations is examined. The original version of FETMOSS was capable of performing the UMS simulations of DG MOSFETs [Abd06]. The validity range wasn’t studied at that time because the RS wasn’t implemented yet. Now, after the RS addition mentioned in chapter three of this thesis, we are able to judge the accuracy of the UMS. In section 4.1, the MS approach as well as the assumptions to transform it to the UMS approach is reviewed. A study of the UMS approach validity for simulating DG MOSFETs is presented in section 4.2. The validity is studied by comparing with the RS results. The modes eigenfunctions and their dependence on the position along the channel and dependence on the gate voltage are also examined. A computationally efficient method to check the UMS validity (without comparing with RS) is then presented in section 4.3. The method was developed after the calculation of the coupling terms and examining them in details. Finally, some ideas reduce the CMS computational burden are suggested in section 4.4.

4.1

MODE-SPACE APPROACH

In the NEGF framework, a suitable set of basis functions is chosen in terms of which the operators, like the Hamiltonian operator and the Green’s function are represented. The MS approach is based on the assumption that the active device is

decoupled from the gate contacts [Ren03]. Decoupling is achieved by applying closed boundary condition at the insulator-metal interface. It may be also applied at the semiconductor-insulator interface if the conduction band offset between the semiconductor and the insulator is taken to be infinite. Applying the closed boundary condition at whether of the two interfaces gives rise to subbands or modes. The subbands are the eigenfunctions associated with the confinement in the gate confinement direction (y-direction in Figure 4.1). The model DG MOSFET is given in Figure 4.1 where it is divided into vertical slices with x spacing in the xdirection. The subbands of the structure are obtained by solving a 1D Schrödinger equation in the y-direction within each vertical slice, positioned at any x  x' , of the device along the x-direction: 

 2  2   n  ( x' , y )  EC ( x' , y )   n  ( x' , y )  E  n  ( x' )   n  ( x' , y ) 2m * y y 2

(4-1)

where  n  ( x' , y) is the eigenfunction in the y-direction at x  x' and E n  (x' ) represents the bottom of the subband at x  x' and n is the subband index.

Z Oxide

X n+

Y

n+

Source

Channel

Drain

Oxide

x Figure 4. 1: A model DG MOSFET divided into vertical slices. The horizontal spacing between the slices is Δx

At each vertical slice x' , the central difference approximation for the second derivative of  n  is applied: n 

2 j y 2



 jn1  2  jn    jn 

y 2

(4-2)

where  jn  =  n  ( x' , y j ) Substituting with eq. (4-2) into eq. (4-1), we obtain:  t y  jn1  2t y  jn   t y  jn1  EC j  jn   E n   jn 

(4-3)

90 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Writing eq. (4-3) at each node along the vertical slice x' and applying closed boundary condition at the semiconductor-insulator interface, a set of linear equations are obtained and cast in the matrix form:  ty 0 2t y  EC 1  t 2t y  EC 1  t y y   0     0  t y 2t y   0  0 

 0   1n      1n    n     n   0   1    1     n       E    (4-4)  0    n    EC N y 1  t y   Nny1   n    Nyn1   ty 2t y  EC N y    N    N y   y 

Eq. (4-4) is an eigenvalue problem that can be solved for N y different subbands. Each subband n has an eigenenrgy E n  and a corresponding eigenvector χ n  within each slice positioned at x' along the x-direction. Only the lowest N m subbands are considered others are discarded; higher subbands are ignored due to their high subband energies E n  and, consequently, low occupation probability. The essence of the MS representation is the expansion of the 2D wave function  ( x, y) in terms of the subbands, obtained from solving eq. (4-4), in the form:

 ( x, y) 

n N m

   ( x)    ( x, y) n

n

(4-5)

n1

where  n  (x) are the expansion coefficients. As indicated in the previous paragraph, only few subbands are taken into consideration and, as a result, the summation is truncated at n  N m . Substituting eq. (4-5) into (3-1), multiplying by  * gives [Fer97]:



m 

( x, y) and integrating over y

n N m  2  2 m  ( x) m  m   E ( x )  E  ( x )  Anm n  ( x)  l * 2 2m x x n 1





(4-6)

where Anm (x) is the operator given by: Anm ( x) 

2 2m * x

   n    2 n  * m  * m  2 dy  ( x , y )  ( x , y )  dy  ( x , y )  ( x, y)   2  x x x  

Eq. (4-6) is the CMS transformation of eq. (3-1). For each mode m , the right hand side involves a summation over all other modes and the m th mode itself. This summation gives rise to coupling between the modes.

The ultra-thin body DG MOSFETs channel is fully depleted and the shape of the potential EC ( x, y) along the y-direction varies slowly with the x- position, hence the variation of  with x can be neglected [Ren03] and eq. (4-6) becomes



 2  2 m  ( x)  E m  ( x) m  ( x)  El m  ( x) * 2 2m x x

(4-7)

Eq. (4-7) is called the uncoupled mode-space transformation of eq. (3-1). Eq. (4-7) represents the 1D transport model equation that will be solved using the NEGF for each subband m

4.1.1 Hamiltonian Matrix in Uncoupled Mode-Space Applying a central finite difference approximation for 2nd order derivative in eq. (4-7) gives: m  m  m   2im  i 1  2i   i 1  x 2 x 2

(4-8)

where im  =  m  ( xi ) Substituting with eq. (4-8) into eq. (4-7), we obtain

 t xim1  2t xim   t xim1  Eim im   Elim 

(4-9)

Writing eq. (4-9) at each node along the 1D grid in the x-direction applying closed boundary condition at the drain contact and source contact, a set of linear equations are obtained and cast in the matrix form: H m φm   E l Iφm 

where 2t x  E1m    tx 0  0   m  2t x  E2  t x 0    tx   Hm   0    0   m  0  t x 2t x  EN x 1  tx    m   0  0  tx 2t x  EN x   N x N x

φ m 

 1 m    m    2         m    N xm1    N x  

(4-10)

(4-11)

92 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

4.1.2 The NEGF in Uncoupled Mode-Space The closed boundary condition was applied at the source and drain contacts to obtain the mode Hamiltonian matrix in eq. (4-10). However, open boundary condition must be used to account for the transport along the x-direction. Within the NEGF framework, the correct (open) boundary condition is incorporated through the use of the drain and source self- energies Dm  and Sm  respectively for each mode m [Ren01]. The m th mode retarded Green’s function G m  is given by:



G m   El I  H m  Σ Sm   Σ Dm 



1

(4-12)

For our simple 1D problem, the self-energy can be obtained from fairly elementary arguments [Dat00]. The closed boundary condition was applied by setting 0m   0

and  Nmx  1  0 . The open boundary condition is applied by setting 0m   1m e ik1x and

 Nm 1   Nm e x

ik N x x

x





El  E1m   2t x 1  cos k1x   E Nmx   2t x 1  cos k N x x where is the Ek dispersion relation in the source and drain contacts. It can be obtained from the solution of the 1D Schrödinger equation inside the source and drain contacts where the subband energy is constant. Accordingly, the drain and source self-energy matrices for the m th mode are given by [Ren03]:

 Dm 

0     0

        0  ik N x x     txe  N x N x

 Sm 

 t x e ik1x 0   0 0           0

 

0

0     0 N  N x x

(4-13)

(4-14)

Once the self energy matrices  Dm  and  Sm  are computed from (4-13) and (4-14), the m th mode broadening functions, ΓSm  and ΓDm  , are obtained as:





(4-15)





(4-16)

Γ Sm   i Sm   S m 

Γ Dm   i Dm   D m 

The m th mode retarded Green’s function G m  is obtained using (4-12) and the m th mode spectral density functions due to the source/drain contacts can be afterwards obtained as:

ASm   G m  Γ Sm G m  ADm   G m  Γ Dm G m 



(4-17)



(4-18)

Then, the m th mode longitudinal-energy-resolved electron density at a grid point i is obtained by [Ren01]: n m  (i; El )   im 

2



1 ASm  i, i; El F ( El , E f S )  ADm  i, i; El F ( El , E f D ) 2x



(4-19)

The distribution function  m  ( x, y ) is used to convert the 2D electron density to the 2

1 comes from the process of 1D discretization x [Dat00]. In a discrete representation the wavefunction is normalized to the number of points N x . The real wavefunction in a continuous representation is normalized to the

usual 3D electron density. The factor

length of the domain N x x . The longitudinal energy resolved electron density

nm  (i, El ) is, further, summed over the N m subbands taken into consideration and the Si six conduction band valleys [Pie03] and, finally, the total electron density n(i) is obtained by integration over the longitudinal energy. The m th mode transmission coefficient from the source contact to the drain contact is defined in terms of the Green’s function and the broadening function as



m  TSD  Trace Γ Sm G m  Γ Dm G 

m 



(4-20)

The m th mode longitudinal energy resolved terminal current is, afterwards, obtained as I  m  ( El ) 



q m  TSD F ( El , E f S )  F ( El , E f D ) 2



(4-21)

over the N m subbands taken into consideration and the Si six conduction band valleys [Pie03] and, finally, the total terminal current is obtained by integration over the longitudinal energy.

4.1.3 UNCOPULED MODE-SPACE VALIDATION BY COMPARISON WITH REAL-SPACE At low drain voltage, the potential along the channel is almost constant and the assumption made in eq. (4-7) is accurate. However, at high drain voltage the assumption may not be valid except for ultra-thin body thickness [Ren03]. For that exception, there is no room for the potential to be changed along the y-direction. This

94 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

makes the shape of the potential along the y-direction is almost constant with the xposition. It was mentioned in [Ven03] that UMS results are accurate for Tsi  5nm as long as the device has a uniform body thickness. We call the body thickness at which the UMS results are no more accurate “the critical body thickness”. The error percentage between the UMS and RS results has not been studied versus the body thickness. Moreover, no information was given about the dependence of the critical body thickness on the gate voltage. The following study aims to explore these points. UMS approach simulation was carried out using FETMOSS [Abd06]. The RS approach has been integrated into FETMOSS and benchmarked as explained in chapter 3. Now we can simulate a given device with both UMS and RS and compare their results. RS results are taken as the reference to study the validity of the UMS. The simulated DG MOSFET device is shown in Figure 3.1 with the device doping, dimensions and material parameters listed in Table 4.1. The supply voltage VDD is taken to be 0.6 V. The body thickness is swept from 1.5 nm to 7 nm. Table 4-1: List of the device parameters used in the simulation together with the model device shown in Figure 3.1. Parameter Dimensions

Doping

Value

channel length

L

10 nm

oxide thickness

Tox

1.5 nm

silicon (body) thickness

TSi

1.5 : 7 nm

source and drain length

LS & LD

5 nm

channel doping

source and drain doping Material

0 cm-3

NA

1020 cm-3

ND

silicon relative permittivity oxide relative permittivity

 Si

11.7

 ox

top and bottom gate work function

3.9

m

0

0

4.25 eV

The error in terminal current between UMS and RS was studied at the two most important bias conditions, the on- and off-state. The on-state is at which both VGS and VDS = VDD while the off-state is at which VGS = 0 V and VDS = VDD . Figure 4.2 depicts the percentage error in the on- and off-current obtained by RS and MS approaches versus the body thickness Tsi . The error percentage is very small, and can be neglected, for ultra-thin body thickness and accordingly the MS approach is accurate for ultra-thin body thickness even at high drain voltage. The error percentage increases with increasing the body thickness and can reach high values for thick body especially at high gate voltage. The critical body thickness is no larger than 3.5nm at VGS  VDD and VDS  VDD . This result is inconsistent with the assumption that the critical body thickness is 5nm as given in [Ven03].

Figure 4. 2: Percentage error in the on-current I on and the off-current I off versus TSi. The percentage difference is calculated as I (%)  1  I UMS / I RS   100 .

In order to better understand the reasons behind the error percentage dependency on the gate voltage and the body thickness, the modes eigenfunctions in the y-direction  n  ( x' , y) are plotted versus y-position at different x-positions for the on- and off-states. Figure 4.3 depicts the 1st order mode eigenfunction  1 at x  0, 4, 8, 12, 16 and 20 nm , TSi  1.5nm in the off-state. The eigenfunction is constant along the x-direction. This is also the case for the on-state as shown in Figure 4.4. The same behavior is found for higher order modes. The 3rd order mode eigenfunction  3 , for example, is given in Figure 4.5 and 4.6 for the same conditions of Figure 4.3 and 4.4. Thus the assumption in eq. (4-7) is quite accurate and so no discrepancy is seen between the UMS and RS results at TSi  1.5nm .

96 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

 1

y (nm) Figure 4. 3: The off-state 1st order mode eigenfunction versus the y -position at different and at TSi  1.5nm

x -positions

 1

y (nm) Figure 4. 4: The on-state 1st order mode eigenfunction versus the y -position at different and at TSi  1.5nm

x -positions

 3 

y (nm) Figure 4. 5: The off-state 3rd order mode eigenfunction versus the y -position at different and at TSi  1.5nm

x -positions

 3 

y (nm) Figure 4. 6: The on-state 3rd order mode eigenfunction versus the y -position at different and at TSi  1.5nm

x -positions

98 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Now, we let’s check the situation at TSi  7nm . Figure 4.7 depicts the 1st order

mode eigenfunction  1 at x  0, 4, 8, 12, 16 and 20 nm , TSi  7nm and at off-state. The eigenfunction is no more constant along the x-direction. This is also the case at the on-state as shown in Figure 4.8. The same behavior is found for higher order modes. The 3rd order mode eigenfunction  3 , for example, is given in Figure 4.9 and 4.10 for the same conditions of Figure 4.7 and 4.8. Therefore, the assumption in eq. (4-7) is not correct leading to the discrepancy seen in Figure 4.2 between the UMS and RS results at TSi  7nm .

 1

y (nm) Figure 4. 7: The off-state 1st order mode eigenfunction versus the y -position at different and at TSi  7nm

x -positions

 1

y (nm) Figure 4. 8: The on-state 1st order mode eigenfunction versus the y -position at different and at TSi  7nm

x -positions

 3 

y (nm) rd

Figure 4. 9: The off-state 3 order mode eigenfunction versus the y -position at different and at TSi  7nm

x -positions

100 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

 3 

y (nm) rd

Figure 4. 10: The on-state 3 order mode eigenfunction versus the y -position at different and at TSi  7nm

x -positions

At TSi  7nm , the error percentage is about 20% and 60% in the off-and onstate respectively. It means larger error percentage VGS  VDD than at VGS  0 . However, Figure 4.7 to 4.10 show that the modes eigenfunctions varies with xposition slower at VGS  VDD than that at VGS  0 . This may look contradicting as slower variation with x-position means less coupling between the modes and, therefore, one may expect the UMS to be more accurate at VGS  VDD than at VGS  0 . Indeed, this is not the case as indicated by the error percentage in Figure 4.2. Seeking the reason for the larger error percentage at the on-state, the contribution of the individual modes to the total current, obtained from simulation results, is plotted in Figure 4.11. It depicts the percentage of the mth subband current I m to the total current I versus Tsi at both off- and on-states. Two deductions can be extracted from Figure 4.11, the first is that the contribution of the higher order modes increases with increasing the body thickness. This occurs because the larger the body thickness, the lower the subband energy and the higher probability for the mode to be occupied. The second deduction is that the contribution of the higher order modes also increases with increasing the gate voltage. This occurs because the higher the gate voltage, the lower the subband energy also. Figure 4.12 depicts the 3rd order mode subband energy, for example, versus the x - position at TSi  7nm and for the off- and onstates.

m=1 m=2 m=3 m=4

Figure 4. 11: Percentage of the mth subband current

I m to the total current I versus Tsi for the off-

and on-states.

x (nm) Figure 4. 12: 3rd order mode subband energy, for example, versus the x - position at for the off- and on-states.

TSi  7nm and

102 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

The dependence seen in Figure 4.12 is strongly expected; it is the concept upon which the transistor operation is built on. The gate voltage controls the conductivity of the channel by controlling the potential inside the channel. The higher the gate voltage, the higher the channel potential and the lower the potential energy and consequently, the lower the subband energies

4.2 COMPUTATIONALLY EFFICIENT METHOD TO CHECK THE UNCOUPLED-MODE SPACE VALIDITY The device design parameters, like the oxide thickness Tox and the metal gate work function  m , greatly affect the device electrostatics and, accordingly, the critical body thickness. The UMS approach is basically used to avoid the computational cost of the RS approach whenever possible. Therefore, it will be inconvenient to estimate the critical body thickness, for each time we tune the device design parameters, by comparing the UMS results against that of the RS. We have to find a faster method that can be used to estimate the critical body thickness, even if the method is not accurate as comparing against the RS approach. The right hand side of eq. (4-6), when ignored, resembles the source of error that appears between the UMS and the RS approach. It can be calculated after the UMS self-consistent loop converges; in order to estimate how large the neglected terms are. These calculations can be done for all the modes taken into account. The ignored terms depend on both the longitudinal energy El and the position along the xdirection. We define the following quantities: 1. coupling effect on the m th mode due to the n th mode Rnm El , x   Anm x n El , x 

(4-22)

2. coupling effect on the mth mode due to the mth mode itself Rmm El , x   Amm x m El , x 

(4-23)

3. total coupling effect on the m th mode:

Rm El , x    Rmn El , x 

(4-24)

n

Figure 4.13 depicts the coupling effect on the lowest order mode due to itself (R11) as a function of the longitudinal energy El and the position along the x direction at TSi  7nm in the off-state. The figure also shows the lowest order mode subband energy (indicated by the solid line) along the x -direction.

Oscillations

Figure 4. 13: R11 as a function of

El and x at TSi  7nm in the off-state. The scale on the right side

of the graph gives the value of the coupling term versus the intensity where brighter regions mean larger value of the coupling term. Moreover, the lowest order mode subband energy by the solid line versus x.

E 1 is indicated

There is no presence for the coupling term below the subband energy. This is identified by the dark region below the solid line. This happens simply because there are no states at this region. Above the subband energy, there is an observed oscillation (indicated by the arrow) in the coupling terms versus the longitudinal energy and the position along the x -direction. These oscillations are identified in the graph by the alternation of bright and dark regions. We relate this oscillation to the well known quantum interference effect. Interference occurs between the electron waves incident and reflected from the energy barrier; the reflected wave interferes constructively or destructively with the incident wave. At high temperature, the quantum interference effect may not be observed in the electron density because it is washed out by the statistics, however, it can be observed at low temperature [Ren03]. The quantum interference effect can be observed in the local density of states (LDOS) before applying the statistics. For example, the 1st mode 2D LDOS at Tsi  7nm in the off-state is shown in Figure 4.14. The interference pattern is clear especially for the drain region. In that region, the energy barrier is relatively large and most of the incident wave from the drain contact is reflected back.

104 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Figure 4. 14:

The 2D LDOS for the off-state 1st order mode versus

Tsi  7nm . E 1 is indicated by the solid line versus x .

El and x -position at

Figure 4.15 depicts the coupling effect on the lowest order mode due to the second order mode (R12) as a function of the longitudinal energy El and the position along the x -direction at TSi  7nm in the off-state. The figure also shows the lowest order as well as the second order modes subband energy (indicated by the solid line) along the x -direction. It is seen that the coupling effect between the two modes has larger values above the top of the energy barrier compared with the rest of the device. Also some oscillations can be identified where dark and bright regions alternate. These oscillations are quantum interference effect. However, it is difficult to explain the oscillation pattern because it is generated from the interaction of two modes, and not from a single mode with itself as in Figure 4.13. A similar behavior was obtained for the higher order coupling terms.

R12 as a function of El and x at Tsi  7nm in the off-state. E 1 and E 2  are indicated by the solid line versus x . Figure 4.15:

R11 as a function of El and x at Tsi  7nm in the on-state. E 1 is also indicated by the solid line versus x . Figure 4.16:

106 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Also Figure 4.16 and Figure 4.17 illustrate the previously shown coupling terms but for the on-states. By comparing the scale (on the right side) in Figure 4.14, 4.15, 4.16 and 4.17, it can be observed that the coupling terms in the off-state of higher values than those at the on-state. This is related to the variations of the subband energy along the x -direction; the subbands energies have faster variations in the offstate than in the on-state. Referring to Figure 4.14 and 4.16, the lowest order mode subband energy has a maximum value  0 eV for the off-state and  0.1 eV for the on-state while the minimum value  0.8 eV for both states.

R12 as a function of El and x at Tsi  7nm in the on-state. E 1 and E 2  are indicated by the solid line versus x . Figure 4.17:

Now, after acquiring a good insight about the coupling terms and their dependence on the gate voltage, the body thickness as well as the x-position, we are ready to develop a quick method to estimate the critical body thickness without comparing with the RS results. The estimation should take into account two important factors, namely, the coupling terms relating different modes and the modes contribution to total current. The coupling terms are energy as well as x-position dependent. We are going to consider the maximum value along the important part of the channel and over the energy range of interest. The terminal current at high VDS is mainly due to electrons incident from the source side and passing above (or through) the top of the energy barrier. Most of the electrons incident from the drain contact are reflected by the relatively large energy barrier. Therefore, the important part of the channel is the region around the energy barrier top. The energy upper limit is determined by the statistics; the condition that the probability of occupation isn’t negligible. The lower limit is taken as the mode subband energy in the drain region; below which there is no allowable electronic states. Figure 4.18 depicts the maximum absolute value of the total coupling effect on the m th mode Rm where

m  1, 2, 3, 4 and 5 versus TSi in the off-state. Also the on-state is given in Figure 4.19. It is observed that the coupling terms are smaller in the on-state than that in the off-state. In contrast, the UMS error is shown in Figure 4.2 to be larger in the on-state. That’s why accounting for the contribution of each mode to the total current is important to estimate the UMS validity.

TSi (nm) Figure 4.18: The maximum absolute value of the total coupling effect on the

m th mode Rm

where

m  1, 2, 3, 4 and 5 versus TSi in the off-state.

TSi (nm) Figure 4.19: The maximum absolute value of the total coupling effect on the

m  1, 2, 3, 4 and 5 versus TSi in the on-state.

m th mode Rm

where

108 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

To account for both the coupling terms and the modes contribution to the current, we define a new quantity R such that:

R

1 I2

I m

m

max |  I n Rnm | n

(4-24)

The new quantity weights the coupling effect the mth mode due to the nth mode by their corresponding percentage contribution in the total current. As shown in Figure 4.20, R has extremely small values for ultra-thin bodies. It starts to increase rapidly at Tsi about 5.1nm and 3.8nm in the off- and on-state respectively. These body thicknesses can be considered the critical values at the corresponding bias conditions. Figure 4.20 shows a percentage error of 5% between at Tsi about 5.2nm and 3.4nm in the off- and on-state respectively. Therefore, this method, based on R calculation, serves as a rough estimate for the UMS validity range at a given device parameters and bias conditions. The important point to be mentioned here is that this method is more computationally efficient than comparing with RS results. Indeed, calculating R almost needs no extra efforts from the simulation procedure of the UMS. Thus enormous saving in the computations is achieved. device designer can use the UMS simulation to tune the device parameters, in order to achieve a given specification. Then the validity range check method is used to know whether the body thickness is considerably smaller than the critical values or not. If yes, the results obtained from UMS are expected to be accurate and no more simulation is needed. If not, an accurate CMS or RS must be used to simulate the device electrical characteristics and makes further tuning if needed.

R(eV )

TSi (nm) Figure 4. 20: Total coupling effect of all the modes weighted by their percentage contribution (as defined in (4-24)) versus TSi in the off and on-state.

4.3

FAST MODE-SPACE SIMULATION

For ultra-thin body thickness, Figure 4.3, 4.4, 4.5 and 4.6 indicate that the modes eigenfunctions are independent on the x-position and the gate voltage. For this reason, we compared the eigenfunctions obtained from the UMS numerical simulation results with the well-known infinite potential well analytical eigenfunctions formulas, with constant potential inside the well, given by [Fer01]:

  n  y  for n odd  Bn cos 2 T n    Si   ( y)    B sin  n y  for n even  n  2T   Si  

(4-25)

where Bn is a normalization constant. In numerical simulation using FD method, the normalization constant Bn is obtained such that:



(n) 2 i

1

(4-26)

Ny

Figure 4.21 depicts the 1st order mode eigenfunction versus the vertical position at TSi  1.5nm . Both simulation result and analytical approximation are shown in the plot. Note that the simulation results are found to be independent on VGS . It is obvious that there is a perfect match between the simulation result and the analytical approximation. This match occurs because for ultra-thin body thickness, the potential along the gate confinement direction; i.e. the y-direction, is almost constant. The same case was found for higher order mode. For example, Figure 4.22 depicts the 3rd order mode eigenfunction for the same conditions. It is also obvious that a perfect match between the simulation result and the analytical approximation exist. Thus for ultrathin bodies, we can use the infinite potential well analytical formulas to reduce the simulation time needed to solve the Schrödinger equation in the gate confinement direction, eq. (4-4). This leads to fast UMS simulation.

110 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

 1

y (nm) Figure 4.21: The 1st order mode eigenfunction versus the vertical position at

TSi  1.5nm . Both

simulation result (SR) and analytical approximation (AA) are shown in the plot. SR obtained at VGS  0 and VGS  VDD are the identical and only one of them is shown.

 3 

y (nm) Figure 4.22: The 3rd order mode eigenfunction versus the vertical position at

TSi  1.5nm . Both

simulation result (SR) and analytical approximation (AA) are shown in the plot. SR obtained at VGS  0 and VGS  VDD are the identical and only one of them is shown

As indicated by Figure 4.2, UMS can’t be used to simulate thick bodies and CMS or RS has to be used instead. In the following discussion, we show how the CMS can be simulated.

Now, consider the maximum absolute values of the coupling terms Rnm at TSi  7nm given in Table 4-2 for the off-state and Table 4-3 for the on-state. There is a common observation among the whole modes; there is always one or two dominant term(s). For example, R31 is the dominate term for the 1st order mode, R13 and R53 are the dominant terms for the 3rd order mode and so on. The dominate terms are highlighted in the tables. Interestingly, this observation was found to be true at other body thicknesses than 7nm . For example, Table 4-4 and Table 4-5 give the maximum amplitude of the coupling terms at TSi  5nm . Table 4-2: The maximum absolute value of the coupling terms Rnm at TSi  7nm , VGS  0 and

VDS  VDD . The dominant terms are highlighted. n m

1 2 3 4 5

1

2

3

4

5

0.0035 0.0033 0.1435 0.0006 0.0199

0.0022 0.0029 0.0011 0.1523 0.0004

0.1445 0.0011 0.0066 0.0011 0.1368

0.0006 0.1527 0.0011 0.0062 0.0010

0.0192 0.0004 0.1233 0.0009 0.0049

Table 4-3: The maximum absolute value of the coupling terms

Rnm at TSi  7nm , VGS  VDD and

VDS  VDD . The dominant terms are highlighted. n m

1 2 3 4 5

1

2

3

4

5

0.0013 0.0030 0.0628 0.0005 0.0125

0.0030 0.0006 0.0009 0.0569 0.0003

0.0585 0.0009 0.0016 0.0008 0.0495

0.0005 0.0570 0.0008 0.0009 0.0009

0.0115 0.0003 0.0453 0.0009 0.0007

112 CHAPTER 4: DOUBLE GATE MOSFETS SIMULATION USING THE NEGF IN MODE-SPACE

Table 4-4: The maximum absolute value of the coupling terms

VDS

Rnm at TSi  5nm , VGS  0 and

 VDD . The dominant terms are highlighted. n

1

2

3

4

5

1 2 3 4 5

0.0013 0.0037 0.0994 0.0004 0.0100

0.0037 0.0011 0.0016 0.0905 0.0003

0.0974 0.0015 0.0021 0.0013 0.0734

0.0004 0.0927 0.0012 0.0017 0.0008

0.0090 0.0002 0.0661 0.0007 0.0012

m

Table 4-5: The maximum absolute value of the coupling terms

Rnm at TSi  5nm , VGS  VDD and

VDS  VDD . The dominant terms are highlighted. n

1

2

3

4

5

1 2 3 4 5

0.0003 0.0033 0.0484 0.0005 0.0072

0.0033 0.0002 0.0011 0.0400 0.0003

0.0436 0.0010 0.0005 0.0008 0.0318

0.0004 0.0382 0.0007 0.0004 0.0004

0.0062 0.0002 0.0273 0.0004 0.0002

m

The existence of few dominant terms among the coupling terms indicates that we may not consider full coupling among the modes during the CMS solution. This will reduce the system of equations we solve and will result in less computation than the full CMS solution. From Table 2, 3, 4 or 5 we deduce that the problem can be decoupled into two smaller problems. Modes 2 and 4 must be solved together due to their significant coupling. Also modes 1, 3 and 5 must be solved together for the same reason. Therefore, we have two uncoupled problems; one problem in two unknowns and the other one in three unknowns. The original problem is in five unknowns and 3

since the solution of a linear system of N unknowns involves N operations [Car95], hence the proposed decoupling is expected to make about 70% reduction in the simulation time for our case.

CHAPTER 5

DESIGN & SIMULATION STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

The International Technology Roadmap for Semiconductors (ITRS) predicts the future device performance by projection of the past technology generations [ITRS]. A major portion of semiconductor device production is devoted to digital logic. The ITRS divided logic ICs into three categories: 1) High Performance (HP) logic technology, 2) Low Operating Power (LOP) technology and 3) Low Stand-by Power technology (LSTP). The transistors for HP ICs have the highest performance but with the highest leakage current while the transistors for LSTP chips have the lowest leakage current but with the lowest performance. In between, the transistors for LOP chips have moderate performance and moderate leakage current. To avoid excessive power dissipation, a common approach for HP ICs is to fabricate more than one type of transistor on the chip, low threshold (HP) devices as well as others with higher threshold and sometimes larger Equivalent Oxide Thickness (EOT) to reduce the leakage current. The HP device is used just in critical paths, and the low leakage devices are used everywhere else. The ITRS is more concerned with HP transistor because this transistor is the technology driver. Table 5-1 summarizes most of the technology requirements for HP, LSTP and LOP DG MOSFETs for the years 2015 and 2022, which represent the end of near-term and long-term requirements respectively [ITRS].

Table 5-1: ITRS the technology requirements for HP, LSTP and LOP DG MOSFETs. HP

LSTP

LOP

Parameter Physical gate length L g (nm)

2015 10

2022 4.5

2015 16

2022 8

2015 13

2022 6

Silicon Thickness TSi (nm) Equivalent Oxide Thickness EOT (Å)

6

3

6

3

6

3

6

5

11

8

8

6

114

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END HP Extension lateral abruptness1 g (nm/dec). Metal gate work function  m

LSTP

LOP

1

0.45

1.6

0.8

1.3

0.6

midg ap

midga p

midgap  0.1eV

midgap  0.1eV

midgap

midgap

0.8

0.65

0.8

0.7

0.6

0.45

2000

4440

0.188

0.375

192

417

0.76

1.2

5.3 10 5

5.78  10 5

NMOS drive current in the onstate I ON ( A / m )

4590

5572

1404

1892

1786

1752

Total gate capacitance

1.05

0.684

1.36

0.84

1.218

0.794

0.180

0.08

0.69

0.31

0.41

0.2

Power supply voltage VDD (V) Maximum gate leakage current density2 J g,limit (A/cm2) Subthreshold offstate leakage current I sd ,leak (

2.2  10 2 4.8 10 2

A / m )

C g ,total

(fF/μm)

NMOSFET intrinsic delay3 τ (ps)

1. Extension

lateral abruptness is defined as the reciprocal of the drain/source doping gradient extending into the channel. 2. Maximum

gate leakage current density is the maximum allowed gate leakage current density at 25°C, and it is measured with the gate biased to VDD and the source, drain, and substrate all tied to ground. 3. NMOSFET

intrinsic delay is the intrinsic transistor delay for NMOS devices at 25°C that is given by: τ = (Cg,total × VDD) / ION.

Figure 5.1 (a) shows the FinFET device considered in this chapter. Within the source and drain regions, denoted by LS and LD respectively in Figure 5.1 (b), the doping is assumed to be constant. The doping profile starting just after LS or LD is assumed to be Gaussian due to the diffusion of the source and drain dopants into the channel. The doping gradient is characterized by the extension lateral abruptness parameter g (nm/dec) [ITRS]. The G-S/D underlap regions are denoted by LSU and LDU respectively. In this chapter, the device is assumed to be symmetric such that LS  LD , LSU  LDU  LU and g S  g D  g .

G G

S

D

TSi

G

Tox

h

G

BOX BOX (a) G

LS

LSU

TSi

(c)

Tox

LDU

LD

G

Lg (b)

Figure 5. 1: (a) 3-D view of a FinFET, (b) A-A′ cross- section and, (c) B-B′ cross-section.

A computational study of the device design issues for 9 nm gate length HP NMOS device was reported using 2D quantum simulation [Has04]. The on-current as well as the intrinsic switching speed were below the ITRS target with SiO2 as gate insulator. The series resistance and the limited scaling of SiO2 were identified as the limiting factors. It was shown that the series resistance is limited by the contact resistance, and the contribution of the source and drain regions to the resistance is negligible when the doping level is 2  10 20 cm 3 . It was suggested that the contact resistance has to be improved or a high-k gate insulator with EOT ≈ 5 Å is necessary to achieve the targets. It was also suggested that a hyper abrupt source/drain doping profile is necessary. However, in [Has04], the operating voltage of 0.4 V was used; whereas the ITRS projected operating voltage for HP 9-nm DG device is 0.8 V [ITRS].

A 2D simulation study for 13 nm gate length DG n-channel FinFET has demonstrated that, by optimizing the Gate-Source/Drain (G-S/D) underlap (Lu), the fin thickness (TSi) of a DG FinFET can be significantly increased without degradation in the speed compared to the conventional G-S/D overlap structure [Woo07]. When underlap exists, the effective channel length (Leff) is significantly longer than Lg in weak inversion, while Leff ~ Lg in strong inversion [Tri05]. Thus, Lg can be decreased

116

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

for higher ON-state current, or TSi can be increased for better manufacturability without degrading SCEs [Woo07]. It was also demonstrated that for optimal performance in the G-S/D underlap, the lateral doping gradient of the S/D does not need to be abrupt.

Recently, with the aid of 2D quantum simulation, a 10 nm gate length DG nchannel FinFET was designed [Has08]. The simulation showed that a full-volume inversion across the entire length of the channel can only be obtained, when the fin thickness TSi is smaller than 0.5 × Lg for the considered 10 nm DG FinFET device. Thus a fin thickness of 4 nm was chosen. To avoid excessive gate leakage, a gate oxide thickness of 1 nm was assumed for the FinFET device. The design was shown to meet the speed requirement of the ITRS but the on-current wasn’t met. The effect of parasitic series source/drain resistance wasn’t included; this will more reduce the on-current and the device switching speed. Finally, the gate material was assumed to be Poly Si. The design and results of this study are summarized in Table 5-2.

Table 5-2: Design parameters and simulation results reported in [Has08]. Parameter

Value

Lg

10 nm

TSi Lu ND NA EOT g VDD

4 nm 8 nm 1020 cm-3 1016 cm-3 10 Å 2 nm/dec 4.35 Ev 0.8 V

J g,limit

2706 A/cm2

I sd ,leak

0.75 μA/μm

I ON

3584 μA/μm 0.454 fF/ μm 0.101 ps

m

Cg,total



In this chapter the optimal design of a 10 nm FinFET will be studied. Some issues will be considered in this design that were ignored in the previously reported designs [Has04, Has08]. It will be shown that, both the on-current and the switching speed requirements can be satisfied simultaneously. Then, the use of a mid-gap work function will be considered. Moreover, it will be shown that device manufacturability can be improved by increasing the fin thickness (TSi) while ITRS targets can still be meet. In this chapter also, the use of high-k material gate dielectric for 2015 LSTP 16 nm gate length FinFET is investigated. For this purpose, a model for gate tunneling current was implemented [Hak08]. The implemented model calculates the gate tunneling current making use of the self-consistent solution results. It will be shown that SiO2 and Si3N2 can’t satisfy the gate leakage target of the ITRS while many other materials can. The least gate current density is found when Lanthanum Oxide (La2O3)

is used as a gate dielectric material. Finally in this chapter, a quantum-based simulation comparison between the 1D and the 2D gate capacitance of the DG MOSFET is carried out. The capacitance calculation was implemented in FETMOSS making use of the 2D self-consistent quantum simulation based on the NEGF. The 1D quantum simulation for the DG MOS capacitor was carried out using an online simulator “Shred” [Vas97, Vas06]. It will be shown that, the 1D capacitance calculation of nanosclae DG MOSFETs is no longer accurate. In section 5.1 the optimal design of a 10 nm FinFET will be studied. Then, in section 5.2 the use of a high-k material gate dielectric for 2015 LSTP 16 nm FinFET is investigated. Finally, a quantum-based simulation comparison between the 1D and the 2D gate capacitance of the DG MOSFET is carried out in section 5.3.

5.1

DESIGN OF A HIGH-POWER 10 NM FINFET

The 10 nm FinFET design reported in [Has08] wasn’t accompanied with an optimization procedure. In our study, the design is optimized to achieve the maximum possible on current or switching speed while keeping the leakage current within the acceptable ITRS limits. The optimization parameters are gate work function, the underlap regions length and the extension lateral abruptness parameter. Other parameters were kept constant and equal to that reported in [Has08]. The design parameters are summarized in Table 5-3. The top gate oxide thickness (see Figure 5.1 (c)) is assumed to be much thicker than the side gate oxide such that channels are formed under the side gate oxides only. Consequently, the simulation domain is assumed to be 2D as shown in Figure 5.1 (b). Quantum effects are taken into consideration using the 2D quantum ballistic NEGF in real-space as the transport model. Table 5-3: Design parameters for our study

Parameter Lg

Value

TSi Lu ND NA EOT g

m

4 nm variable 1020 cm-3 1016 cm-3 10 Å variable variable

VDD

0.8 V

10 nm

The optimization procedure is carried out as follows. The extension lateral abruptness parameter was varied from g  2 to g  5 nm/dec. For a given extension lateral abruptness, the underlap regions length is varied. Then, for a given combination of the extension lateral abruptness and the underlap regions length, the work function  m is found such that the source/drain subthreshold off-state leakage

118

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

current is just below the roadmap limit ( 0.76A / m ). The extension lateral abruptness values were chosen such that advanced fabrication process isn’t needed [Shi03, Kuz03]. The values used for the work function were chosen to be practically manufacturable [Li05].

Figure 5.2 depicts the needed work function versus the underlap length at different extension lateral abruptness. As the underlap region length is increased, the effective channel length is increased and the off-state leakage current is decreased. Thus the work function, needed for a given off-state leakage current, decreases. The same thing occurs when the extension lateral abruptness is decreased. Once the needed work function is found for each combination of the extension lateral abruptness and the underlap regions length, the on-current can be directly found. Figure 5.3 depicts the on-current versus the underlap length at different extension lateral abruptness. It is clear that for each extension lateral abruptness, there is an optimum value for the on-current. The existence of such an optimum-value behavior was explaneind in [Has04]. For a specific doping profile, if Lu is too small, it increases the SCE and thus the off-current. In order to meet the off-current specification, higher work function gate metal is necessary. This high value of gate work function reduces the gate overdrive at on-state, reducing the on-current. On the other hand, if Lu be too large, the gate cannot efficiently modulate the barrier of the whole channel, which also reduces the on-current. As a result an optimum value of LU can be obtained in between, this gives maximum on-current

Figure 5. 2: The metal gate work function needed to satisfy ITRS off-current (< 0.76A / m ) versus the underlap length at different extension lateral abruptness for the design in table 5.3.

Figure 5. 3: The on-current versus the underlap length at different extension lateral abruptness for the design in table 5.3. The work function was set to satisfy the ITRS off-current target (< 0.76A / m ).

It can be deduced from Figure 5.3 that, the optimum ballistic on-current value increases slightly with increasing the value of the extension lateral abruptness. It is expected that a similar behavior will be found even if the transport wasn’t ballistic. This expectation is based on the fact that, when the underlap region length is too long and/or extension lateral abruptness is too abrupt, the additional resistance degrades the on-current [Woo07]. It is important to note that the optimum on-current value occurs at a longer underlap for a less-abrupt (larger) extension lateral abruptness. Thus, the effect of the additional resistance may be cancelled out and both non-ballistic and ballistic current behave the same. The investigation of the optimum value of the oncurrent taking into account non-ballistic effects was left as future work. Table 5-4 summarizes the parameters that give the optimal on-current and the corresponding switching speed. In order to calculate the switching speed, the gate capacitance calculation was implemented in the FETMOSS simulator. The total charge (per unit width) in the device Q is calculated after the convergence of the selfconsistent solution and the capacitance is obtained from the relation:

C gg  Q / VG

(5-1)

where VG is taken to be 0.01 V around VGS  0.8 V and VDS is kept constant at 0.8 V.

120

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END Table 5-4: Design parameters for optimal on-current and the resulting switching speed. The rest of parameters are given in table 5.3. g (nm/dec) 2 3 4 5

Lu (nm) 4 5 7 8

 m (eV) 4.337 4.374 4.372 4.436

Ioff (μA/μm) 0.74 0.75 0.74 0.75

Ion (μA/μm) 6248 6469 6525 6734

Cgg (fF/μm) 0.383 0.375 0.356 0.341



(ps)

0.049 0.046 0.044 0.041

f T (THz) 20.4 21.7 22.7 24.4

It is interesting to note that, the obtained switching speed is about 4 times greater than that specified in the ITRS. There are some reasons for this gain: 1) we used an EOT of 1 nm in our design while is it specified to be 0.5 nm in the ITRS; this results in reducing the gate capacitance by a factor of two, 2) quantum mechanical effects further reduces the capacitance below the value estimated from parallel plate capacitance calculations, 3) the series resistance wasn’t included in our calculations, and finally 4) the on-current optimization carried out in our design. The ITRS requires the use of a mid-gap gate material for the FinFET devices. By mid-gap materials we mean those whose Fermi level is near the intrinsic energy level of Si that corresponds to a work function around 4.61 eV. The advantage of a mid-gap material is that it can be used as a gate for both n-channel and p-channel devices [Col08]. From the simulation study we carried out, we can deduce the design parameters that can be used in conjunction with a mid-gap gate material. This is done by drawing a horizontal line at 4.61 eV in Figure 5.2, and then finding the intersection with the drawn curves and finally rounding the overlap length (found from the intersection) to the nearest integer. Table 5-5 summarizes the parameters that can be used together with a near mid-gap material and satisfy the off-current target of the ITRS. It can be noticed that both the on-current and switching speed targets are satisfied for the different cases. Table 5-5: Design parameters for near mid-gap work function and the resulting on-current and switching speed. The rest of parameters are given in table 5.3. g (nm/dec) 2 3 4 5

Lu (nm) 1 2 4 6

 m (eV) 4.605 4.677 4.621 4.587

Ioff (μA/μm) 0.74 0.74 0.74 0.74

Ion (μA/μm) 5150 5384 5869 6162

Cgg (fF/μm) 0.420 0.4 0.32 0.3



(ps)

0.060 0.059 0.044 0.039

f T (THz) 16.7 16.9 22.7 25.64

The manufacturability of the FinFET is limited by the fin thickness and the extension lateral abruptness of the source/drain doping profile. A relaxed abruptness has been presented in the designs given in Table 5-4 and 5-5. It worth mentioning that the presented designs outperform the ITRS on-current and switching speed targets while keeping the off-current within the allowed limit. Extension lateral abruptness of 2, 3, 4 and 5 nm/dec have been used while the ITRS predicts the use of 1 nm/dec abruptness. Thus, novel dopant annealing processes [Shi03, Kuz03] may not be needed for the 10 nm FinFET device fabrication. In order to further increase the manufacturability of the FinFET, possibility of increasing the fin thickness is now investigated. In the designs presented in [Has04] and [Has08], the used fin thickness was 3 and 4 nm respectively while the ITRS target is 6 nm fin thickness. It has been demonstrated in [Woo07] that the fin thickness may

be increased without degradation in SCEs when underlap is used. Moreover, the use of low body doping makes the FinFET more immune to random dopant fluctuations and, thus, reduces its exposure to process variation [Fer06]. Exploiting these features, the procedure used to obtain the designs in Table 5-4 and 5-5 is repeated to design a FinFET with a 6 nm fin thickness and 1010 cm-3 body doping. The design parameters are given in Table 5-6. The results encountered during the optimization procedure are shown in Figure 5.4 and 5.5. Table 5-6: Design parameters for better manufacturability design.

Parameter Lg TSi Lu ND NA EOT G m VDD

Value 10 nm 6 nm variable 1020 cm-3 1010 cm-3 10 Å variable variable 0.8 V

Figure 5. 4.: The metal gate work function needed to satisfy ITRS off-current (< 0.76A / m ) versus the underlap length at different extension lateral abruptness for the design in table 5.6.

122

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

Figure 5. 5: The on-current versus the underlap length at different extension lateral abruptness for the design in table 5.6. The work function was set to satisfy the ITRS off-current target (< 0.76A / m ).

Table 5-7 summarizes the design parameters for optimal on-current and the resulting switching speed using 6 nm fin thickness. By comparing the results in Table 5-4 and Table 5-7, one can observe that the optimum FinFET on-current is reduced with increasing the fin thickness. This happens because the increased fin thickness causes a larger off-current (or smaller threshold voltage) due to the less quantization of the subbands energy. The increased off-current is decreased back to its ITRS target value by increasing the metal work function. The increased work function, by its turn, reduces the optimum value of the on-current. Nevertheless, both the on-current and switching speed values presented in Table 5-7 still satisfy the ITRS targets. This indicates that the manufacturability of the FinFET can be enhanced by increasing the fin thickness to 6 nm while keeping the device performance satisfactory from the ITRS point of view. Table 5-7: Design parameters for optimal on-current and the resulting switching speed using 6 nm fin thickness. The rest of parameters are given in table 5.6. g (nm/dec)

Lu (nm)

 m (eV)

Ioff (μA/μm)

Ion (μA/μm)

Cgg (fF/μm)

2 3 4 5

5 6 8 10

4.438 4.488 4.497 4.502

0.754 0.749 0.749 0.74

4788 4937 5130 5224

0.40 0.33 0.35 0.34



(ps)

0.067 0.053 0.055 0.052

f T (THz) 14.9 18.9 18.2 19.2

Finally, the use of a near mid-gap work function is also examined with a 6 nm fin thickness. The results are given in Table 5-8. The given values indicate that, it is possible to use a mid-gap work function with a 6 nm fin thickness while still

satisfying the ITRS targets. Again, comparing the results in Table 5-5 and Table 5-8 shows that, the optimum FinFET on-current is reduced with increasing the fin thickness. Table 5-8: Design parameters for near mid-gap work function and the resulting on-current and switching speed. The rest of parameters are given in table 5.6. g (nm/dec) 2 3 4 5

Lu (nm) 3 5 7 9

 m (eV)

Ioff (μA/μm) 0.75 0.756 0.748 0.744

4.591 4.567 4.571 4.572

Ion (μA/μm) 4273 4644 4905 5103

Cgg (fF/μm) 0.34 0.30 0.35 0.33



(ps)

0.064 0.052 0.057 0.052

f T (THz) 15.6 19.2 17.5 19.2

5.2 STUDY OF HIGH-K MATERIAL FOR LOWSTANDBY-POWER 16 NM FINFET The original version of FETMOSS doesn’t take into account the gate current in simulation. The most accurate way for the gate current calculation is to include it in the self-consistent solution using the NEGF, as indicated in chapter 3. This is, however, out of the scope of this study and was left as future work. In this section, the use of high-k material gate dielectric for 2015 LSTP 16 nm gate length FinFET is investigated using a simple gate current model [Hak08]. The model taken from the literature was originally for bulk MOSFETs. This model was adapted for FinFET (or DG MOSFET) and implemented in FETMOSS. The implemented model calculates the gate tunneling current from the results of the self-consistent solution (but not included in the self-consistent solution itself). The model is based on Tsu Esaki relation [Duk69, Tsu73] given by: 4qm * max TC ( E ) N S ( E )dE h 2 Emin E

Jg 

(5-2)

where E is the energy component normal to the Si-SiO2 interface, TC (E ) is the transmission coefficient calculated from the numerical solution of Schrödinger equation inside the dielectric material using the self-consistent potential and N S (E ) is called the supply function. The supply function describes the supply of carriers for tunneling and is given by: 

N S ( E )    f1 ( Etotal )  f 2 ( Etotal )dE parallel

(5-3)

0

where f1 is the energy distribution function in the channel while f 2 is that in the gate electrode. The supply function was assumed to be Fermi Dirac distribution with the Fermi level replaced by its quasi counterpart calculated from the self-consistent solution results. Consequently, the supply function is given by:

124

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

E  Ef1    1  exp( ) K BT   N S ( E )  K BT ln  E  Ef 2  )  1  exp(  K BT  

(5-4)

The minimum energy E min in eq. (5-2) is the minimum allowed energy for electrons tunneling from or to the conduction band. In bulk MOSFETs, the minimum energy should be the conduction band edge. In FinFETs, however, carrier confinement leads to energy quantization. Therefore, the minimum energy becomes the first subband energy. The maximum energy E max in eq. (5-2) is the tunneling barrier top energy after which thermoionic emission current flows. The total gate current is obtained by summing the currents of all subbands.

The model was implemented in the FETMOSS simulator and was verified by comparing with a published data on the gate current [Sai07]. The comparison results are shown in Figure 5.6. There is a good agreement between the presented model and the reference data within 9 % RMS difference. This is acceptable for our study where good estimation of the gate current is needed to compare different high-k materials.

Figure 5.6. Verification of our model results against that in [Sai07]. The device parameters and dimensions are: Lg= 20 nm, Tsi= 5 nm, Tox= 1 nm, NA= 1017 cm-3. The drain and source were connected together.

Now, the use of high-k material gate dielectric for 2015 LSTP 16 nm gate length FinFET is investigated. Table 5-1 lists parameters projected by ITRS for this device. The dielectric materials included in our investigation with their dielectric constant and conduction band offset from Si is listed in Table 5-9 [Nag05].

Table 5-9: Relative dielectric material constant and conduction band offset for dielectric materials used in our simulation taken from [Nag05]. Material Al2O3 HfO2 La2O3 SiO2 Si3N4 Ta2O5 TiO2 Y2O3 ZrO2

Dielectric constant 9 25 30 3.9 7 25 40 15 25

Conduction band offset (eV) 2.8 1.5 2.3 3.0 2.0 1.4 1.1 2.3 1.4

Figure 5.7 compares the gate current density using those dielectric materials against the Jg,limit of the ITRS indicated by the horizontal line. The materials below that line have less current density than the maximum allowed one. Unfortunately SiO 2 and Si3N2, that are most compatible with Si, don’t satisfy the gate leakage target of the ITRS. The least gate current density is found when Lanthanum Oxide (La2O3) is used as a dielectric material followed by Hafnium Oxide (HfO2). Due to technological reasons and compatibility issues, the material that Intel used in its recent 45 nm technology was HfO2 [Kuh08].

Figure 5. 7. High-k material study for LSTP 16nm technology node in the year 2015.

126

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

5.3 COMPARISON OF 1D AND 2D DG MOSFET CAPACITANCE Accurate calculation of the DG MOSFET capacitance is an important factor in estimating the device switching speed. 1D quantum simulation for the capacitance of a DG MOS capacitor can be carried out using an online simulator “Shred” [Vas97, Vas06]. Due to the 1D assumption, Shred doesn’t account for the S/D regions that reflect the 2D nature of the capacitance. Therefore, capacitance calculation was implemented in FETMOSS making use of the 2D self-consistent quantum simulation based on the NEGF. The total charge (per unit width) in the device Q is calculated after the convergence of the self-consistent solution and the capacitance is obtained from the relation C gg  Q / VG . A series of simulations were carried out sweeping the gate length from 5 nm to 50 nm where 2D effects are expected to be more pronounced for shorter gate lengths. Other simulation parameters are given in Table 5.10. The underlap region was assumed of zero length and S/D junctions were assumed to be abrupt. Moreover, it was assumed that the device is under equilibrium by connecting the S/D together and observing the effect of varying the gate voltage on the gate capacitance. All the mentioned assumptions were imposed in order to match the assumptions of Shred. It is important to note that, the FETMOSS is now capable of simulating the gate capacitance under non-equilibrium. This capability arises from the use of the NEGF in the capacitance calculations. The investigation of non-equilibrium gate capacitance and its bias and geometry dependent behavior was left as a future work.

Table 5-10: Simulation parameters for 2D capacitance study Parameter

Lg TSi Lu ND NA EOT g

m VDD

Value 5:50 nm 4 nm none 1020 cm-3 1016 cm-3 10 Å abrupt 4.35 0.8 V

Figure 5.8 depicts the gate capacitance per unit width versus VGS at different Lg obtained from FETMOSS. At high VGS, the capacitance increases as Lg increases as a result of the larger gate area. At low VGS the capacitance is small (but not zero) due to spillover of electrons from the S/D regions to the channel. The S/D doping was kept constant together with the channel doping that result in a constant capacitance in the spillover region. Consequently, the non-zero capacitance at zero VGS is constant capacitance for all the simulated gate lengths.

Cgg (fF/μm)

In order to compare the results obtained with FETMOSS to those of Shred, the capacitance is obtained per unit area (length and width) and the results are depicted in Figure 5.9. It is noticed that the capacitance of the 50 nm gate length device is relatively close to Shred capacitance due to the little impact of the 2D effects. For short gate lengths the 2D effects plays a major role in determining the device capacitance and 1D simulation becomes no more accurate. Since the zero VGS capacitance is almost constant (as seen in Figure 5.8) for different Lg, the capacitance per unit area of the shortest device is the greatest one at zero VGS.

Figure 5. 8: The gate capacitance per unit width versus VGS at different Lg obtained from FETMOSS

128

CHAPTER 5: DESIGN AND SIMULATION

STUDY OF DOUBLE-GATE MOSFETS AT THE ROADMAP END

Cgg (fF/μm2)

x 10

Figure 5. 9. The gate capacitance per unit area versus VGS at different Lg obtained from FETMOSS and compared to Shred.

From this study we conclude that 1D capacitance calculation of DG MOSFETs is no more accurate even for a 50 nm gate length. Interaction of the source/drain electric fields with the channel electric field is strong enough to make the 1D capacitance calculation invalid.

CONCLUSIONS AND FUTURE WORK A Two-Dimensional (2D) quantum transport based simulator for Double-Gate (DG) MOSFETs was developed. The simulator is based on the qausi-2D simulator “FETMOSS” originally developed by Dr. Tarik Abdolkader in his PhD thesis. True 2D quantum simulation was enabled using the Non-Equilibrium Green’s Function Formalism (NEGF) in the Real-Space (RS) representation. A computationallyefficient method based on Gauss elimination was proposed and implemented. The proposed method showed an exact match with the traditional NEGF while the simulation time was reduces by about an 85 %.

As one of the computationally efficient methods available in the literature, the Contact Block Reduction (CBR) method was reviewed and implanted in FETMOSS. Two sample device dimensions compatible with the ITRS were chosen and simulated using the CBR method. The method was shown to be computationally efficient only when the device is operating in the on-state. This is because only 6% of eigenstates is needed for acceptable accuracy. In the off-state, however, the method gives good accuracy with at least 40% of eigenstates, otherwise it gives wrong value for the tunneling current. Thus the CBR method is rendered in-efficient for the off-state simulations. The Uncoupled Mode-Space (UMS) validity for simulating DG MOSFETs was examined. This examination was carried out by comparing the UMS results to that of the RS added to the FETMOSS simulator. The UMS validity range was reported to depend on the device parameters and bias conditions. It was shown that the UMS may lose its validity for silicon thickness less than 5nm. A computationally efficient method to check the validity of the UMS, without comparing with RS, was developed. The method is based on calculating the modes coupling coefficients, the contribution of each mode to the total current and accordingly estimates the validity range of the UMS approach. Moreover, the possibility of reducing the Coupled Mode-Space (CMS) computational burden was investigated. Few coupling terms were found to be significant and , thus, other terms can be neglected to partially decouple the problem into simpler ones. This partial decoupling is expected to achieve about 70% reduction in the simulation time 2D quantum transport based simulation and design optimization of a 10 nm high power FinFET was presented. It was shown that both the on-current and the switching speed targets of the International Technology Roadmap for Semiconductor (ITRS) can be satisfied simultaneously while keeping the off-current within the specified limits. It was also shown that the fin thickness can be increased and the extension lateral abruptness can be relaxed while keeping the device performance within the limits. As a result, the FinFET manufacturability can be enhanced. Moreover, it was shown that it is possible to use mid-gap work function gate electrodes while still satisfying the ITRS targets.

130

CONCLUSIONS AND FUTURE WORK

A model for the tunneling gate current in the bulk MOSFET was successfully adapted to suits DG MOSFETs. It was implemented in FETMOSS and was used to study the effect of using high-k materials on the gate leakage in low standby power 16 nm FinFET. The gate current obtained for SiO2 and Si3N4 was shown to be exceeding the allowable limits of the ITRS. La2O3 and HfO2 were found to be the most promising materials showing the lowest gate current leakage. Finally, the effect of 2D electrostatics on the gate capacitace of the DG MOSFETs was studied. The study eas carried out by adding the capacitance calculations capability to the FETMOSS simulator. Then, 1D capacitace was calculated using the online available simulator “Shred”. It was shown that 1D capacitance isn’t valid even for a 50 nm DG MOSFET.

Future work based on the progress achieved in this thesis can be summarized in the following points: 1. Taking into account different scattering mechanisms in the NEGF simulation using the RS representation. 2. Gate current calculation self-consistently with the NEGF using the RS representation. 3. A broader study of the CBR accuracy for various device dimensions and operating conditions. 4. Implementing the CMS representation in FETMOSS and examine the validity of the partial decoupling scheme suggested. 5. Accounting for the technological and compatibility issues and restrictions when studying the impact of using high-k material on the gate current. 6. Optimizing the LSTP transistor for minimum gate-leakage current 7. Investigation of the non-equilibrium gate capacitance and its bias and geometry dependence. 8. Extending the capabilities of FETMOSS to simulate other nanoscale devices like multi-gate MOSFETs and carbon nanotubes.

REFERENCES [Abd04]

T. M. Abdolkader, W. Fikry, H.H. Hassan and O.A. Omar. “Solution of Schrödinger equation in double-gate MOSFETs using transfer matrix method,” IEE Electronics Letters, vol.40, pp. 1307-1308, 2004

[Abd06]

T. Abdolkader, W. Farouk, O. Omar and M. Hassan. “FETMOSS: software tool for 2D simulation of double-gate MOSFET,” Int. J. Numer. Model., vol.19, pp. 301-214, 2006

[Ahm07]

S. Ahmed, Gerhard Klimeck, Derrick Kearney, Michael McLennan, MP Anantram, "Quantum Simulations of Dual Gate MOSFET Devices: Building and Deploying Community Nanotechnology Software Tools on NanoHUB.org," J. High Speed Electron., in press (2007).

[Ale03]

Alexander L. Fetter, John Dirk Walecka, Quantum Theory of Many Particle Systems, Dover Publications, United States of America, 2003

[Ana06]

M. Anantram, M. Lundstrom, D. Nikonov ”Modeling of Nanoscale Devices,” Condensed Matter, Mesoscopic Systems and Quantum Hall Effect, 2006 [online]. Available: arXiv:cond-mat/0610247v2.

[Ant88]

P. Antognetti and G. Massobrio, Semiconductor Device Modeling with SPICE (McGraw- Hill, New York, 1988).

[Ban05]

I. Banerjee, M. Morse, J.-M. Verdiell and D. Docter, “Materials in semiconductor processing for high speed circuits: from electrical to optical circuits,” IEE Proc.-Circuits Devices Syst., Vol. 152, No. 5, October 2005 pp: 523-526.

[Bans05]

A. Bansal, B. C. Paul, and K. Roy, “Modeling and optimization of fringe capacitance of nanoscale DGMOS devices,” IEEE Trans. Electron Devices, vol. 52, no. 2, pp. 256 – 262, Feb 2005.

[Bog98]

Bogdan Majkusiak, Tomasz Janik, and Jakub Walczak, “Semiconductor Thickness Effects in the Double-Gate SOI MOSFET,” IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 45, NO. 5, MAY 1998 1127

[Boh07]

Mark Bohr, Kaizad Mistry and Steve Smith, Intel Demonstrates Highk + Metal Gate Transistor Breakthrough on 45 nm Microprocessors, Avialble: well.ee.nthu.edu.tw/~ycking/Course/Indus07/20070418_Intel_Highk.p df

132

REFERENCES

[Car95]

T. A. Carter, R. A. Tapia & A. Papakonstantinou: Linear Algebra, An Introduction to Linear Algebra for Pre-Calculus Students, Rice University, May 1995 [online]. Available: http://ceee.rice.edu/publications.html

[Che99]

B. Cheng, M. Cao, R. Rao, A. Inani, P. Vande , W. Greene, J. Stork, Z. Yu, P. Zeitzoff, and J. Woo, “The Impact of High- Gate Dielectrics and Metal Gate Electrodes on Sub-100 nm MOSFET’s,” IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 46, NO. 7, JULY 1999

[Col90]

J.P. Colinge, M.H. Gao, A. Romano-Rodriguez, H. Maes, and C. Claeys, “SILICON-ON-INSULATOR "GATE- ALL-AROUND DEVICE,” IEDM IEEE 90

[Dam03]

P. Damle “Nanoscale device modeling: From MOSFETs to molecules,” Ph.D. Dissertation, Purdue University, West Lafayette, May 2003

[Dan84]

Danielewics, P. “Quantum theory of non-equilibrium processes”, Ann. Phys. 152, 239, 1984

[Dat95]

Supriyo Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, Cambridge, UK, 1995

[Dat89]

Supriyo Datta, Quantum Phenomena, Modular Series on Solid-state Devices, vol. VIII, New York, Addison Wesley, 1989

[Dat00]

S. Datta, “Nanoscale device modeling: the Green’s function method,” Superlatt. Microstruc., vol. 28, pp. 253-278, 2000

[Dat05]

S. Datta, Quantum Transport University Press, 2005

[Dig90]

Digh Hisamoto, Toru Kaga, Yoshifumi Kawamoto and Eiji Takeda, “A Fully Depleted Lean-Channel Transistor (DELTA)-A Novel Vertical Ultrathin SO1 MOSFET,” IEEE ELECTRON DEVICE LETTERS, VOL. 11, NO. 1, JANUARY 1990

[Drag98]

Dragica Vasileska, William J. Gross and David K. Ferry, “Modeling of deep-submicrometer MOSFETs: random impurity effects,threshold voltage shifts and gate capacitance attenuation” ,IWCE-6. Extended Abstracts of 1998 Sixth International Workshop on Computational Electronics, Oct 1998, pp: 259-262 Dragica Vasileska and Stephen M. Goodnick, Computational Electronics, 2006 by Morgan & Claypool

[Dra06]

[Duk69]

Atom to Transistor, Cambridge

C. B. Duke, “Tunneling in Solids”. Academic Press, 1969.

[Fer97]

D. Ferry and S. Goodnick, Transport in Nanostructures, Cambridge University Press, Cambridge, UK, 1997

[Fer01]

D.Ferry, Quantum Mechanics: An Introduction For Device Physicists And Electrical Engineers, 2nd edition, Institute of Physics Publishing, 2001

[Fer06]

D. K. Ferry, R. Akis, A. Cummings, M. J. Gilbert, and S. M. Ramey, “Semiconductor Device Scaling: Physics, Transport, and the Role of Nanowires” IEEE-NANO June 2006, Vol.2, pp: 415- 418

[Fra98]

David J. Frank, Yuan Taur and Hon-Sum P. Wong, “Generalized Scale Length for Two-Dimensional Effects in MOSFET’s”, IEEE ELECTRON DEVICE LETTERS, VOL. 19, NO. 10, OCTOBER 1998 385, pp: 385-387

[Fra01]

DAVID J. FRANK, ROBERT H. DENNARD, EDWARD NOWAK and PAUL M. SO ,”Device Scaling Limits of Si MOSFETs and Their Application Dependencies” PROCEEDINGS OF THE IEEE, VOL. 89, NO. 3, MARCH 2001

[Fra02]

D. J. Frank, “Power-constrained CMOS scaling limits” IBM J. RES. & DEV. VOL. 46 NO. 2/3 MARCH/MAY 2002. Online: http://www.research.ibm.com/journal/rd/462/frank.txt

[Gor97]

C. M. Goringe, D. R. Bowler and E. Hernández, “Tight-Binding Modelling of Materials” Rep. Prog. Phys. 60 1447 (1997).

[Hac02]

Hackenbuchner S, Sabathil M, Majewski J A, Zandler G, Vogl P, Beham E, Zrenner A and Lugli P 2002 Physica B 314 145

[Hak08]

M. A. Hakim “Simulation & Parameter Extraction of double gate MOSFET: Modeling Of gate Leakage Current,” B.Sc. Dissertation, Ain Shams University, Cairo, July 2008

[Has04]

S. Hasan, J.Wang, andM. Lundstrom, “Device design and manufacturingissues for 10 nm-scale MOSFETs: A computational study,” Solid State Electron., vol. 48, pp. 867–75, 2004.

[Has08]

Hasanur R. Khan, Denis Mamaluy and Dragica Vasileska, “Approaching Optimal Characteristics of 10-nm High-Performance Devices: A Quantum Transport Simulation Study of Si FinFET, “IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 55, NO. 3, MARCH 2008.

[Han89]

W. Hansch, T. Vogelsang, R. Kircher, and M.Orlowski, ”Carrier Transport near the Si/SiO2 interface of a MOSFET,” Solid State Electron., vol. 32, no. 10, pp. 839–849, 1989.

134

REFERENCES

[Hau96]

H, Haug and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, Springer-Verlag Berlin Heidelberg 1996

[Hay00]

William H. Hayt and John A. Buck, Engineering Electromagnetics, sixth edition, McGraw-Hill, 2000

[Hu96]

C. Hu, “Gate oxide scaling limits and projection,” in IEDM Tech. Dig., 1996, p. 319.

[ITRS]

International Technology Roadmap for Semiconductors. Availabe: http://www.itrs.net/about.html

[Jan98]

T. Janik, and B. Majkusiak, “Analysis of the MOS transistor based on the self-consistent solution to the Schrodinger and Poisson equations and on the local mobility model,” IEEE Trans. Electron Devices, vol. 45, no. 6, pp. 1263–1271, June 1998.

[Jau94]

Jauho A.P., Wingereen N. S., and Meir Y. “Time-dependent transport in interacting and non-interacting resonant tunneling systems” Phys. Rev. B, 50, 5528, 1994

[Jon01]

Jong-Tae Park, Jean-Pierre Colinge, and Carlos H. Diaz,” Pi-Gate SOI MOSFET” IEEE ELECTRON DEVICE LETTERS, VOL. 22, NO. 8, AUGUST 2001

[Kad62]

Kadanoff, L.P. and Baym, G., Quantum Statistical Mechanics, Frontiers in Physics Lecture Note Series, Benjamin/ Cummings, 1962

[Kah60]

D. Kahng and M. M. Atalla, “Silicon–silicon dioxide field induced surface devices,” presented at the IRE Solid-State Device Res. Conf., Pittsburgh, PA, June 1960.

[Kaw99]

H. Kawaura, T. Sakamoto, and T. Baba, “Direct source–drain tunneling current in subthreshold region of sub-10-gate EJMOSFETs,” Si Nanoelectronics Workshop Abstracts, June 1999, pp. 26–27.

[Key75]

Keyes, R. W., "Effect of Randomness in the Distribution of Impurity Ions on FET Thresholds in Integrated Electronics", Journal of SolidState Circuits, August, 1975, pp. 245- 247.

[Kli98]

G. Klimeck, D. Blanks, and R. Lake, Nanoelectronic Modeling Tool (NEMO), User’s Manual and Reference Guide, Raytheon, 1998.

[Kli02]

Gerhard Klimeck, Fabiano Oyafuso, Timothy B. Boykin, R. Chris Bowen, and Paul von Allmen , "Development of a Nanoelectronic 3-D (NEMO 3-D) Simulator for Multimillion Atom Simulations and Its Application to Alloyed Quantum Dots" (INVITED), Computer

Modeling in Engineering and Science (CMES) Volume 3, No. 5 pp 601-642 (2002). [Kha06]

H. Khan, D. Mamaluy, D. Vasileska , “Self-consistent Treatment of Quantum Transport in 10 nm FinFET Using Contact Block Reduction (CBR) Method,” Journal of Computational Electronics , 2006

[Kha07]

Hasanur Rahman Khan, “QUANTUM TRANPORT SIMULATION OF NANO-SCALE FINFET DEVICES,”Ph.D dissertation , ARIZONA STATE UNIVERSITY, 2007.

[Kha87]

Khan, F. S., Davies, J. H. and Wilkins, J. W. “Quantum transport equations for high electric fields” Phys. Rev. B, 36, 2578, 1987

[Kik05]

Kikuji Hirose , Tomoya Ono, Yoshitaka Fujimoto and Shigeru Tsukamoto, First-Principles Calculations In Real-Space Formalism: Electronic Configurations And Transport Properties Of Nanostructures, Imperial College Press, Jan 2005

[Kim04]

Hyun-Woo Kim, Ji-Young Lee, Jangho Shin, Sang-Gyun Woo, HanKu Cho, and Joo-Tae Moon, " Experimental Investigation of the Impact of LWR on Sub-100-nm Device Performance ", IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 51, NO. 12, DECEMBER 2004

[Kim05]

K. Kim, K. K. Das, R. V. Joshi, and C.-T. Chuang, “Leakage Power Analysis of 25-nm Double-Gate CMOS Devices and Circuits,” IEEE Trans. Electron Devices, vol. 52, no. 5, pp. 980–986, May 2005.

[Kyl65]

Kyldysh, L.V. “Diagram technique for non-equilibrium processes”, Sov. Phys. JETP20, 1018, 1965.

[Kuh08]

K. Kuhn, C. Kenyon, A. Kornfeld, M. Liu, A. Maheshwari, W. Shih, S. Sivakumar, G. Taylor, P.Voorn and K. Zawadzki “Managing Process Variation in Intel’s 45nm CMOS Technology,” Intel Technology Journal Volume 12 pp. 93:109, June 2008

[Kuz03]

V. I. Kuznetsov, A. J. M. M. van Zutphen, H. R. Kerp, P. G. Vermont, X. Pages, J. J. van Hapert, K. van der Jeugd, and E. H. A. Granneman, “Continuity in the development of ultra shallow junctions for 130-45 nm CMOS: The tool and annealing methods,” in Proc. IEEE Int. Conf. RTP, 2003, pp. 63–74.

[Lan76]

Langreth , D. C. , In Linear and Non-liner Electron Transport in Solids, eds. J. T. Devereese and E. van Doran, NATO Advanced Study Institute Series B, vol.17, p.3, Plenum, New York, 1976

136

REFERENCES

[Lee05]

K. Lee, C.-H. Lee, D. Park, and Y.-K. Choi, “A Study of NegativeBias Temperature Instability of SOI and Body-Tied FinFETs,” IEEE Electron Device Letters, vol. 26, no. 5, pp. 326–328, May 2005.

[Lel00]

Leland Chang, Peter Yu, “Ballistic Transport in Silicon MOSFET’s,” Physics 250 Term Paper, May 2000 [online]. Available: http://hkn.eecs.berkeley.edu/~leland/coursework/phys250.pdf

[Lem04]

M.C. Lemme, T. Mollenhauer, W. Henschel, T. Wahlbrink, M. Baus, O.Winkler, R. Granzner, F. Schwierz, B. Spangenberg and H. Kurz, Subthreshold behavior of triple-gate MOSFETs on SOI Material,” Solid State Electronics 48-4, 529 (2004).

[Li05]

Li T-Z, Hu C-H, Ho W-L, Wang C-H and Chang C-Y, “Continuous and precise work function adjustment for integrable dual metal gate CMOS technology using Hf-Mo binary alloys,” IEEE Trans. Electron Devices vol. 52, 2005, pp. 1172–1179.

[Lib98]

Richard Liboff, Introductory Quantum mechanics, Addison Wesley, 1998.

[Lil30]

J. E. Lilienfeld, “Method and apparatus for controlling electric currents,” U.S. Patent 1 745 175, 1930.

[Lo97]

S.-H. Lo, D. A. Buchanan, Y. Taur, and W. Wang, “Quantummechanical modeling of electron tunneling current from the inversion layer of ultra-thin-oxide nMOSFETs,” IEEE Electron Device Lett., vol. 18, p. 209, May 1997.

[Loj07]

Bo Lojek, History of Semiconductor Engineering, Springer-Verlag Berlin Heidelberg 2007

[Lun00]

Mark Lundstrom, Fundamental of carrier transport, Cambridge University Press, Cambridge, UK, 2000

[Lun06]

Mark Lundstrom and Jing Guo, Nanoscale Transistors: Device Physics, Modeling and Simulation, Springer Science Business Media Inc. 2006

[Maj01]

H. Majima, Y. Saito and T. Hramoto, “Impact of Quantum Mechnaical Effects on Design of Nano-scale Narrow Channel n- and p-type MOSFETs”, IEEE IEDM March, 2001 pp: 733:736

[Mah87]

Mahan, G. D. “Quantum transport equation for electric and magnetic fields”, Ohys. Rep. 145, 251, 1987

[Mam03]

D. Mamaluy, M. Sabathil, and P. Vogl, “Efficient method for the calculation of ballistic quantum transport,” JOURNAL OF APPLIED PHYSICS VOLUME 93, NUMBER 8 15 APRIL 2003

[Mam04]

DENIS MAMALUY, ANAND MANNARGUDI AND DRAGICA VASILESKA , “Electron Density Calculation Using the Contact Block Reduction Method, “Journal of Computational Electronics 3: 45–50, 2004 Kluwer Academic Publishers. Manufactured in The Netherlands.

[Mam05]

D. Mamaluy , D. Vasileska, M. Sabathil, T. Zibold, and P. Vogl,” Contact block reduction method for ballistic transport and carrier densities of open nanostructures,” PHYSICAL REVIEW B 71, 245321, 2005

[Mam07]

Hasanur R. Khan, Denis Mamaluy and Dragica Vasileska “Quantum Transport Simulation of Experimentally Fabricated Nano-FinFET, “IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 54, NO. 4, APRIL 2007

[Mam08]

D. Mamaluy (private communication), 2008

[Mat02]

MATLAB, The language of Technical computing, Release 13, The Mathworks, Inc., 2002.

[Mat04]

L. Mathew, Yang Du ,A.V-Y Thean, M. Sadd, A. Vandooren, C. Parker, T. Stephens, R. Mora, R. Rai, M. Zavala, D. Sing, S. Kalpat, J. Hughes, R. Shimer, S. Jallepalli, G. Workman, W. Zhang2,J.G. Fossum2, B. E. White,B. Y. Nguyen, J. Mogab., “CMOS Vertical Multiple Independent Gate Field Effect Transistor (MIGFET),”Proceeding IEEE International DOI Conference, 187 (2004).

[Mart59]

Martin, P. C. and Schwinger, J. “Theory of many-particle systems. I”, Ohys. Rev. 115, 1342, 1959

[Muk06]

Mukunda P.Das, Frederick Green “Ballistic transport is dissipative: the why and how,” J. Phys.: Condens. Matter 17 V13-V16, 19 Jan 2006 [online]. Available: http://arxiv.org/pdf/cond-mat/0601459

[Nag05]

P.V. Nagaraju, Amitava DasGupta, “Study of gate leakage current in symmetric double gate MOSFETs with high-k/stacked dielectrics,” Science Direct, Oct.2005.

[NEM]

NEMO 1-D: the first NEGF-based TCAD cobweb.ecn.purdue.edu /~gekco /pubs_src/P50.pdf

[Pas06]

V. Passi, B. Olbrechts, J.P. Raskin, “Fabrication of a Quadruple Gate MOSFET in Silicon-on-Insulator technology,” Abstarcts of the NATO Advanced Research Workshop on Nanosclaed Semiconductir-onInsulator Structures and Devices, 11 (2006).

[Pie03]

R. Pierret, Advanced Semiconductor Fundamental, Prentice Hall, 2003

tool.

Available:

138

REFERENCES

[Ram86]

Rammer, J. and Smith, H. “Quantum field-theoretical methods in transport theory of metals”, Rev. Mod. Phys. 58, 323, 1986

[Rah05]

Anisur Rahman, Mark S. Lundstrom, and Avik W. Ghosh, “Generalized effective-mass approach for n-type metal-oxidesemiconductor field-effect transistors on arbitrarily oriented wafers,” J. Appl. Phys. 97, 053702, 2005

[Ren01]

Z. Ren “Nanoscale MOSFETs : physics, simulation and design” Ph.D. Dissertation, Purdue University, West Lafayette, October 2001

[Ren03]

Z. Ren, R. Venugopal, S. Goasguen, S. Datta, and M. S. Lundstrom, “nanoMOS 2.5: A Two-dimensional simulator for quantum transport in Double-gate MOSFETs,” IEEE Trans. Electron Devices, vol. 50, pp. 1914–1925, Sep. 2003.

[Ric01]

C. Richter, A. Hefner, and E. Vogel , “A Comparison of QuantumMechanical Capacitance–Voltage Simulators”, IEEE ELECTRON DEVICE LETTERS, VOL. 22, NO. 1, JANUARY 2001 35.

[Sab03]

M. SABATHIL AND S. BIRNER, D. MAMALUY and P. VOGL, “Efficient Computational Method for Ballistic Currents and Application to Single Quantum Dots,” Journal of Computational Electronics 2: 269–273, 2003

[Sab04]

M. Sabathil, D. Mamaluy and P. Vogl , “Prediction of a realistic quantum logic gate using the contact block reduction method,” INSTITUTE OF PHYSICS PUBLISHING SEMICONDUCTOR SCIENCE AND TECHNOLOGY, Semicond. Sci. Technol. 19 (2004) S137–S138 PII: S0268-1242(04)72179-3

[Sai07]

Saibal Mukhopadhyay, Keunwoo Kim, Jae-Joon Kim, Shih-Hsien Lo, Rajiv V. Joshi, Ching-Te Chuang, Kaushik Roy, “Estimation of gateto-channel tunneling current in ultra-thin oxide sub-50 nm double gate devices”, Microelectronics Journal, Vol. 38, pp. 931-941 August-Sep. 2007.

[Sek84]

T. Sekigawa and Y. Hayashi, “Calculated Threshold-Voltage Characteristics of an XMOS Transistor Having an Additional Bottom Gate,” Solid State Electron. 27, 827 (1984).

[Sha01]

Sharma, Rajnish K, Kumar, Ashok, Anthony, John M, “Advances in high-k dielectric gate materials for future ULSI devices,” Jome, June01 [online].Available:http://findarticles.com/p/articles/mi_qa5348/is_2001 06/ai_n21473896/pg_1

[Shi03]

A. Shima, H. Ashihara, T. Mine, Y. Goto, M. Horiuchi, Y. Wang, S. Talwar, and A. Hiraiwa, “Self-limiting laser thermal process for ultrashallow junction formation of 50-nm gate CMOS,” in IEDM Tech. Dig., Dec. 2003, pp. 493–496.

[Sin06]

N. Singh, A. Agarwal, L.K. Bera, T.Y. Liow, R. Yang, S.C. Rustagi, C.H. Tung, R. Kumar, G.Q. Lo, N. Balasubramanian, D.L. Kwong, “ High performance fully depleted silicon nanowire (diameter