3D Simulation of Nanowire FETs using Quantum Models
Vijay Sai Patnaik, Ankit Gheedia and M. Jagadesh Kumar
The authors are with the Department of Electrical Engineering, Indian Institute
of Technology, Huaz Khas, New Delhi 110 016, India (Email: vijaysai.patnaik@gmail.com)
Abstract— After more than 30 years of validation of Moore’s law, the CMOS technology has already entered the nanoscale (sub100nm) regime and faces strong limitations. The nanowire transistor is one candidate which has the potential to overcome the problems caused by short channel effects in SOI MOSFETs and has gained significant attention from both device and circuit developers. In addition to the effective suppression of short channel effects due to the improved gate strength, the multigate NWFETs show excellent current drive and have the merit that they are compatible with conventional CMOS processes. To simulate these devices, accurate modeling and calculations based on quantum mechanics are necessary to assess their performance limits, since crosssections of the multigate NWFETs are expected to be a few nanometers wide in their ultimate scaling. In this paper we have explored the use of ATLAS including the Bohm Quantum Potential (BQP) for simulating and studying the shortchannel behaviour of nanowire FETs.
Introduction
Some of the fundamental problems of MOSFETinspired devices for sub10nm channel length are expected to be electrostatic limits, sourcetodrain tunneling, carrier mobility, process variations, and static leakage. A simultaneous concern is that of power scaling. Nevertheless. it appears that emerging device architectures can extend the CMOS lifetime and provide solutions to continue scaling into the nanometer range, or at least until the 10 nm wall is reached [13]. But what comes after this limit? The nanowire transistor is one candidate which has gained significant attention from both device and circuit developers because of its potential for building highly dense and high performance electronic products. Nanowire transistors can be made using different materials on low cost substrates such as glass or plastics. Si and Ge nanowire transistors are of particularly more importance because of their compatibility with CMOS technology[47].
The objective of this paper is to explain how ATLAS3D can be used to simulate silicon nanowire field effect transistors. We demonstrate that by using the Bohm Quantum Potential to include the confinement effects and by calibrating the BQP model against the 2D SchrodingerPoisson simulation, the scaling behaviour of silicon nanowire FETs can be successfully predicted.
Quantum Models in ATLAS3D
Quantum3D provides a set of models for simulation of various effects of quantum confinement and quantum transport of carriers in semiconductor devices such as the silicon nanowire FETs shown in Fig. 1. A selfconsistent Schrödinger – Poisson solver allows calculation of bound state energies and associated carrier wave function self consistently with electrostatic potential.
Fig.1 3D Nanowire FET generated using ATLAS. 
A. Selfconsistent Coupled SchrodingerPoisson Model [8]
To model the effects of quantum confinement, ATLAS also allows the solution of Schrödinger’s Equation along with the fundamental device equations. The solution of Schrödinger’s Equation gives a quantized description of the density of states in the presence of quantum mechanical confining potential variations.
The SCHRO parameter of the MODEL statement enables the selfconsistent coupled SchrödingerPoisson model for electrons. The P.SCHRO parameter in the MODEL statement enables the SchrödingerPoisson model for holes. ATLAS2D can solve Schrödinger equation in 1D slices or 2D plane. In the case of cylindrical coordinates, Schrödinger equation is solved in radial direction for different orbital quantum numbers and for all slices perpendicular to the axis. By default, a 1D equation is solved. To enable the 2D model, we need to specify the option 2DXY.SCHRO on the MODEL statement. ATLAS3D solves a 2D Schrödinger equation in yz slices, perpendicular to x axis and the option 2DXY.SCHRO is optional.
When the quantum confinement is in one dimension (along yaxis), the calculation of the quantum electron density relies upon a solution of a 1D Schrödinger equation solved for eigen state energies E_{iv}(x) and wavefunctions _{iv}(x, y) at each slice perpendicular to xaxis and for each electron valley (or hole band) v as described below:
(1)
Here, (x, y) is a spatially dependent effective mass in ydirection for the vth valley and EC(x, y) is a conduction band edge. The equation for holes is obtained by substituting hole effective masses instead of electron ones and valence band edge E_{V}(x,y) instead of E_{C}(x, y). Analogously, in cylindrical coordinates, the equation for the radial part of the wave function R_{im}v(r) for each orbital quantum number m reads as
(2)
In the case where the mass is isotropic (i.e., mass is the same in all directions), only one solution to Schrödinger’s equation is obtained with the appropriate mass. ATLAS will automatically determine the appropriate number of valleys based on the material in question. We can, however, limit the number of directions in k space by using the NUM.DIRECT. To choose the number of electron valleys, we need to specify the NUM.DIRECT parameter on the MODELS statement. If you set NUM.DIRECT to 1, we will obtain the solution for isotropic effective mass (MC on the MATERIAL statement). If we set NUM.DIRECT to 3, we obtain a solution for a single lateral effective mass and two transverse masses (ML, MT1 and MT2 on the MATERIAL statement). In a special case of a 1D confinement and equivalent transverse masses, we can set NUM.DIRECT to 2 and only two solutions will be obtained (ML and MT1) with appropriate degeneracy factors. To specify how many valence bands to consider, we need to specify the NUM.BAND parameter. If we set NUM.BAND to 1, you will obtain a Schrödinger solution for only one valence band (MV on the MATERIAL statement). If you set NUM.BAND to 3, you will cause solutions for heavy holes, light holes and holes in the split off band (MHH, MLH, and MSO on the MATERIAL statement).
Using FermiDirac statistics, the discrete nature of the quantized density of states reduces the integral over energy to a sum over bound state energies. The expression for the electron concentration then becomes:
(3)
for 1D confinement and
(4)
for a 2D confinement case and cylindrical case, where F1/2 is the FermiDirac integral of order 1/2.
ATLAS solves the onedimensional Schrödinger’s Equation along a series of slices in the y direction relative to the device. The location of the slices in the y direction is developed in two ways. For rectangular ATLASdefined meshes, the slices will automatically be taken along the existing mesh lines in the ATLAS mesh. If the mesh is nonrectangular or not an ATLAS defined mesh or both, we must specify a rectangular mesh. To do this, we need to specify the locations of individual mesh lines and their local spacing using the SX.MESH and SY.MESH statements like that to the specification of a device mesh using the X.MESH and Y.MESH or a laser mesh using the LX.MESH and LY.MESH statements. The 2D Schrödinger equation in xy plane can be solved on a general nonrectangular mesh. Specifying a rectangular mesh in addition to imported mesh is unnecessary.
Once the solution of Schrödinger’s Equation is taken, carrier concentrations calculated from Equation (3) and Equation (4) are substituted into the charge part of Poisson’s Equation. The potential derived from the solution of Poisson’s Equation is substituted back into Schrödinger’s Equation. This solution process (alternating between Schrödinger’s and Poisson’s equations) continues until convergence and a selfconsistent solution of Schrödinger’s and Poisson’s equations is reached.
Since the wave functions diminish rapidly from the confining potential barriers in the Schrödinger solutions, the carrier concentrations become small and noisy. We can refine these carrier concentrations by setting a minimum carrier concentration using the QMINCONC parameter on the MODELS statement. This parameter sets the minimum carrier concentration passed along to the Poisson solver and the output to the structure files. The transition between the Schrödinger solution and the minimum concentration is refined between 10¥QMINCONC and QMINCONC so that it is continuous in the first derivative.
We can use the SAVE statement or the OUTFILE parameter on the SOLVE statement to write the solutions of the selfconsistent system into a structure file. These structure files will then contain the selfconsistent potential and electron or hole concentrations. The Eigen energies and functions can also be written to the structure file by specifying the EIGENS parameter on the OUTPUT statement. This parameter specifies the number of Eigen energies/wave functions to be written.
The number of Eigen values solved is limited to a number of 2 less than the total number of grid points in the Y direction. Note that the selfconsistent solution of Schrödinger’s Equation with the Poisson’s Equation doesn’t allow solutions for the electron and hole continuity equations in the current ATLAS version. Nonselfconsistent solutions, however, can be obtained by setting the POST.SCHRO parameter in the MODELS statement. These nonselfconsistent solutions are obtained by solving Schrödinger’s Equation only after convergence is obtained. That way, you can obtain Schrödinger solutions with the electron and hole continuity equations.
B. Bohm Quantum Potential (BQP) [9]
This model was developed for SILVACO by the University of Pisa and has been implemented into ATLAS with the collaboration of the University of Pisa. This is an alternative to the Density Gradient method and can be applied to a similar range of problems. There are two advantages to using Bohm Quantum Potential (BQP) over the density gradient method. First, it has better convergence properties in many situations. Second, you can calibrate it against results from the SchrödingerPoisson equation under conditions of negligible current flow.
The model introduces a position dependent Quantum Potential, Q, which is added to the Potential Energy of a given carrier type. This quantum potential is derived using the Bohm interpretation of quantum mechanics and takes the following form
(5)
where and are two adjustable parameters, M_{1} is the inverse effective mass tensor and n is the electron (or hole) density. This result is similar to the expression for the quantum potential in the density gradient model with = 0.5, but there are some differences about how they are implemented.
The Bohm Quantum Potential (BQP) method can also be used for the Energy balance and hydrodynamic models, where the semiclassical potential is modified by the quantum potential the same way as for the continuity equations.
The iterative scheme used to solve the nonlinear BQP equation along with a set of semiclassical equations is as follows. After an initial semiclassical solution has been obtained, the BQP equation is solved on its own Gummel iteration to give Q at every node in the device. The semiclassical potential is modified by the value of Q at every node and the set of semiclassical equations is then solved to convergence as usual (using a Newton or Block iterative scheme). Then, the BQP equation is solved to convergence again and the process is repeated until selfconsistency is achieved between the solution of the BQP equation and the set of semiclassical equations. The set of semiclassical equations solved can be any of the combinations usually permitted by ATLAS.
To use the BQP model for electrons (or holes), we need to specify BQP.N (BQP.P) in the MODELS statement. We can also set the parameter values (a and g) and the direction of the quantization (confinement).
C. Calibration against SchrödingerPoisson Model
We can obtain close agreement between BQP and the results of SchrödingerPoisson (SP) calculations for any given class of device. ATLAS has a SchrödingerPoisson model that can model spatial confinement in only one direction. Therefore, calibration is currently restricted to this case. To obtain comparisons with SP results, ATLAS recommends to use either the new quasistatic capacitancevoltage profile feature or compare chargevoltage curves. This will ensure similar charge control properties between the two models.
The first part of the calibration is to choose a suitable biasing for the device. There should be negligible current flow and quantum confinement effects that manifest at the chosen biases. The second part of the calibration is to set the appropriate BQP parameters in the MATERIAL or MODELS statements, and to set CARRIERS=0 in the METHOD statement.
This will cause the BQP equation to be coupled with Poisson’s equation using the charge density terms.
(6)
This gives the same results as solving the current continuity equations with the constraint of zero current density.
The third part of calibration is to choose the quantity to compare with SP results. For example, for a MOSFET holding the drain and source voltages at the same bias and ramping the gate bias will give us a bias dependent capacitance with negligible current flow. So for an NMOS, we may have the statement
SOLVE VGATE=0.0 NAME=GATE VSTEP=0.01 VFINAL=2.0 QSCV
to give us the quasistatic CV curve as it is biased into inversion. It is best to use a fine voltage step with QSCV to give good resolution. The process can be repeated by setting the SP model in the MODELS statement instead of BQP to obtain the same set of curves for the SP model. The BQP model is then rerun with different sets of parameters until an acceptable agreement with the curves produced by the SP model is achieved.
In our simulations, we have used the Quasistatic capacitance voltage curves for comparison with the Schrödinger Poisson equations. Later on, the electron concentration variation with depth was also compared. It was found that calibration against either one of these quantities also gave the best possible calibration when the other quantity was compared. Hence, comparing both was unnecessary.
D. Post Calibration runs [8]
After obtaining the parameters for BQP, either the DriftDiffusion or energy balance (hydrodynamic) equations can be solved as usual. For cases where Lattice Heating is important then LAT.TEMP can be enabled at the same time.
The iteration scheme uses a modified version of BLOCK. Set BLOCK in the METHOD statement (although NEWTON and GUMMEL are ignored, BLOCK is always used if the BQP model is set). If an Energy Balance model is chosen (HCTE.EL or HCTE.HO on the MODELS statement), then an alternative iteration scheme will become available by specifying BQP.ALTEB in the METHOD statement. This method is slower and is only available in ATLAS2D but may have better convergence properties.
By using BQP.NOFERMI, the BQP equation will only use its Boltzmann statistics form. Without this parameter, the statistics used are those specified for the other equations. With fermi statistics, the convergence can be poor for very high carrier densities, and this parameter can circumvent some convergence properties. Recalibrate the BQP parameters if you set BQP.NOFERMI.
To speed up convergence, we specified the NOCURRENT parameter on the first SOLVE statement after SOLVE INIT. It prevented the need to use the QFACTOR parameter as was necessary for the Density Gradient method. QSCV enables the quasistatic capacitance calculation and output.
3D Quantum Simulation Results
In this section, we present the results of the 3D quantum simulation of a silicon nanowire structure.
A. Device structure
The simulated device is a threedimensional structure with gates all around the silicon channel as shown in Fig. 1 and its cutplane is shown in Fig.2. Device simulations are performed using Silvaco’s 3D ATLAS device simulation environment. Half of the device is constructed in a 2D platform and then extended in the z plane to create a 3D nanowire structure for simulations. The device cross section is 5 nm X 5 nm (rectangular) while its effective channel length is varied between 5 and 100 nm. The channel region is assumed to be p doped (NA = 10^{16} cm^{3}) and the source and drain extension regions are heavily ndoped (N_{D} = 10^{20} cm^{3}) with abrupt doping pro€les. A metallic gate with the midgap work function of 4.7 eV is assumed and the gate oxide is 2 nm thick.
Fig.2 YZ cutplane view of the Nanowire.

B. Quantum Simulation using BQP Model
As discussed earlier, the Bohm Quantum Potential (BQP) model [9] calculates a position dependent potential energy term using an auxiliary equation derived from the Bohm interpretation of quantum mechanics. This extra potential energy modifies the electron and/or hole distribution. Fig. 3 shows an isosurface of electron Bohm Quantum Potential value of 0.1V and and Fig. 4 shows its cutplane across the center of the channel showing the electron BQP variation. It is observed that there is localization near the perimeters of the device. The effect of the EBQP is to reduce the electron density around the edges of the channel, where EBQP is negative and increase electron density in the center, where EBQP is positive. The 2D electron concentration is plotted in Fig. 5 for 0.5 V bias applied to the gate. It clearly shows electrons are repelled from the interface Si/Oxide in all four directions equally.
Fig.3 An isosurface of constant electron Bohm Quantum
Potential with a value 0.1 V.

Fig.4 Cutplane across the center of the Silicon Nanowire
FET channel showing the 2D electron BQP variation.

Fig. 5 Cutplane across the center of the Silicon
Nanowire FET channel showing the 2D electron concentration variation.

C. Calibration of BQP model
For calibration of the BQP model with Schrodinger Poisson results, we chose to compare the quasistatic capacitancevoltage curves. In ATLAS, the low frequency CV curve can be computed in static operation: the charge concentration is integrated in the whole structure and then this quantity is derived as C=dQ/dV. In Fig. 6 we show the calibration of Bohm Quantum Potential model against the SchrodingerPoisson model for a NWFET of channel length 50 nm being biased into inversion. The calibration takes place under conditions of negligible current flow. The voltage range was 01V volts which includes the depletion to inversion transition region.
Fig.6 Calibration of BQP parameters for the 50 nm
channel length Nanowire FET.

First, the Quasistatic CV curves using the SchrodingerPoisson model were evaluated. Quantum con€nement is strong in both the x and y directions, so a complete solution of the 2D Schrodinger equation in two dimensions is required. QSCV is specified on the solve statement and a fine voltage increment is used in order to give accurate calculation of the quasistatic capacitance as a function of voltage. Next, the CV curves over the same voltage range and with the same voltage increment as for the SP model are simulated. Multiple different sets of BQP parameters are used for each of these simulations, in order to get a close calibration as shown in Fig. 7 where the closest overlap is obtained for gamma = 0.01 and alpha = 0.3 for the 50 nm channel length nanowire FET. The electron concentration profiles were also compared using the calibrated BQP model values to confirm the calibration as shown in Fig. 8 obtained at V_{DS} = 0.5 V.
Fig.7 Total quasistatic CV showing the closest match
between BQP and SP.

As can be observed in Fig. 8, the effective quantum potential method is capable to generate accurate shape of the carrier density pro€le as the applied voltage is varied. The method is also capable to accurately reproduce the carrier density per unit area obtained with the direct solution of the Schrödinger equation. For the particular device whose calibration is performed (in this case, the 50 nm length NWFET), the calibrated set of BQP parameters are used for all the further calculations.
Fig. 8 Plot of electron concentration variation across
the channel using the calibrated BQP model.

D. Device characteristics using calibrated BQP model
The transfer and output characteristics for the device were obtained using the calibrated BQP model. The effect of using the BQP model is shown in Fig. 9 by comparing against the results obtained from semi classical simulation, which does not include quantum confinement effects. The effect of BQP model is to increase the drain current compared to the current obtained when quantum confinement is not included. The output characteristics computed for different gate voltages using the BQP model are shown in Fig. 10.
Fig. 9 Transfer characteristics computed for a 50
nm channel lengh nanowire FET using semiclassical and quantum methods
(V_{DS} = 0.5 V).

Fig. 10 Output characteristics computed VGS = 0.6,
.7 and 0.8 V for L = 50 nm using BQP model.

E. Channel Length scaling
We studied the scaling behavior of SNWFETs by varying the length of the channel in the device. Fig. 11 shows the variation in transfer characteristics when the channel length is reduced from 100 nm to 10 nm. The cross section of the channel was 5 nm x 5nm in each case. It can be seen that the device characteristics degrade as the channel length is reduced. As the channel length becomes smaller, the drain current does not saturate, as is seen for the case of 10 nm channel length device shown in Fig. 12 for V_{GS} = 0.6 V. This is due to the Drain Induced Barrier Lowering (DIBL).
Fig. 11 Transfer characteristics of a Silicon Nanowire
FET with a cross section 5 nm x 5 nm and for channel lengths L =
100, 50, 40, 30, 20, 15, and 10 nm.

Fig. 12 Output characteristics of a Silicon Nanowire
FET with a cross section 5 nm x 5 nm and for channel lengths L = 100,
50, 40, 30, 20, 15, and 10 nm.

To understand the DIBL effects, the conduction band energy along the channel from source to drain is shown in Fig. 13 for different drain voltages for a 20 nm channel length nanowire FET. This figure clearly demonstrates the lowering of the barrier potential on the source side as the drain voltage increases. Threshold voltage variation as a function of channel length shown in Fig. 14 also demonstrate the lack of gate control due to short channel effects.
Fig. 13 Conduction band energy at Y=0 as a function
of X at different drain voltages V_{DS} = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 V
and V_{GS} = 0.6 V.

Fig. 14 Threshold voltage variation as a function
of channel length of a Silicon Nanowire FET with a cross section 5 nm
x 5

Conclusions
ATLAS Quantum provides a set of models for simulation of various effects of quantum confinement and quantum transport of carriers in semiconductor devices. We have performed a quantum simulation of Silicon Nanowire FieldEffect Transistors by using the Bohm Quantum Potential model to include quantum confinement effects in the simulations. The BQP model was successfully calibrated against 2D Schrodinger Poisson simulation results under conditions of negligible current flow by matching the Quasistatic Capacitance Voltage curves. The scaling behavior of the silicon nanowire €eld effect transistors has been investigated, and it has been found that their device characteristics sensitively depend on the channel length. The device performance has been found to degrade sharply for channel lengths in the sub 20 nm regime.
References
 H.S. P. Wong, “Beyond the conventional transistor,” Solid State Electronics, vol. 49, pp. 755762, May 2005.
 J. T. Park, J. P. Colinge, “Multiplegate SOI MOSFETs: device design guidelines,” IEEE Trans. Electron Devices, vol. 49, no. 12, pp. 2222 2229, Dec. 2002.
 J. P. Colinge, J.W. Park, and W. Xiong, “Threshold voltage and subthreshold slope of multiplegate SOI MOSFETs,” IEEE Electron Device Lett., vol 24, no 8, pp. 515517. Aug. 2003.
 Y. Cui, Z. Zhong, D. Wang, W. Wang, and M. Lieber, “High performance silicon nanowire €eld effect transistors,” Nano Lett., vol. 3, pp. 149152, 2003.
 D. K. Ferry, “Quo Vadis Nanoelectronics,” Phys. Stat. Sol. (c), vol. 5, no. 1, pp. 1722, 2008.
 D. R. Bowler, “Atomicscale nanowires: Physical and electronic structure,” J. Phys. Cond. Matt., vol. 16, pp. R721R754, 2004.
 S.M. Koo, M. D. Edelstein, L. Qiliang, C. A. Richter, and E. M. Vogel, “Silicon nanowires as enhancementmode Schottky barrier fieldeffect transistors,” Nanotechnology, vol. 16, pp. 14821485, 2005.
 ATLAS User’s Manual Volume I, Silvaco International, March 2007.
 G. Iannaccone, G. Curatola, G. Fiori, “Effective Bohm Quantum Potential for device simulators based on drift diffusion and energy transport”, Simulation of Semiconductor Processes and Devices 2004, vol. 2004, pp. 275278, Sept. 2004.