A 2-D/3-D Schrödinger-Poisson Drift-Diffusion Numerical Simulation of Radially-Symmetric Nanowire MOSFETs

1. Introduction

The phenomenal success of CMOS technology, and, then the progress of the information technology, can be attributed without any doubt to the scaling of the MOS transistor, which has been pushed during more than thirty years to increasingly levels of integration and performances. Then, MOSFETs have been fabricated always smaller, denser, faster and cheaper in order to provide ever more powerful products for digital electronics. Recently, the scaling rate has accelerated, and the MOSFET gate length is now less than 40 nm, with devices entering into the nanometer world [1]-[2]. The so-called “bulk” MOSFET is the basic and historical key-device of microelectronics: its dimensions have been reduced more than ~10³ times during the three past decades. However, the bulk MOSFET scaling has recently encountered significant limitations, mainly related to the gate oxide (SiO₂) leakage currents [3]-[4], the large increase of parasitic short channel effects and the dramatic mobility reduction [5]-[6] due to highly doped Silicon substrates precisely used to reduce these short channel effects. Technological solutions have been proposed in order to continue to use the “bulk solution” until the 45 nm ITRS node. Most of these solutions envisage the introduction of high-permittivity gate dielectric stacks (to reduce the gate leakage, [4], [7]-[8]), midgap metal gate (to suppress the Silicon gate polydepletion-induced parasitic capacitances) and strained Silicon channel (to increase carrier mobility, [9]). However, in parallel to these efforts, alternative solutions to replace the conventional bulk MOSFET architecture have been proposed and studied in the recent literature. These options are numerous and can be classified in general according to three main directions: (i) the use of new materials in the continuity of the “bulk solution”, allowing increasingly MOSFET performances due to their dielectric properties (permittivity), electrostatic immunity (SOI materials), mechanical (strain), or transport (mobility) properties; (ii) the complete change of the device architecture (e.g. multiple-gate devices, Silicon nanowires MOSFET) allowing better electrostatic control, and, as a result, intrinsic channels with higher mobilities and currents; (iii) the exploitation of certain new physical phenomena that appear at the nanometer scale, such as quantum transport, substrate orientation or modification of the material band structure in devices/wires with nanometer dimensions [2], [10].

Nanowire MOSFETs with completely surrounding gate are presently considered as one of the possible solutions for replacing bulk (single-gate) devices and continuing MOSFET scaling in the nanometer scale [11]-[12]. The main advantage of this architecture is to offer a reinforced electrostatic coupling between the conduction channel and the gate electrode, which considerably reduces short-channel effects (SCE) compared to conventional bulk devices [11]-[13]. Then, the constraints on channel doping levels can be relaxed and nanowire MOSFET devices can be designed with intrinsic (pure Si) channels. This offers considerable advantages, especially in terms of mobility and elimination of doping fluctuations. The intrinsic channel can also be beneficial with regard to source-to-drain transport due to the high probability of ballistic transport.

Nanowire MOSFETs are generally designed with very thin silicon films in order to reinforce the electrostatic control of the gate over the channel. The ultra-thin silicon film creates a sufficiently narrow rectangular potential well for inducing the quantization of carrier energy. Carriers are then confined in a rectangular quantum well having feature size close to the electron wavelength. This gives rise to a splitting of the energy levels into subbands (two-dimensional (2-D) density-of-states) [14]-[15], such that the lowest of the allowed energy levels for electrons (resp. for holes) in the well does not coincide with the bottom of the conduction band (resp. the top of the valence band). In nanowire devices the carrier energy is quantified in the two directions perpendicular to the source-to-drain axis (leading to one-dimensional (1-D) density-of-states). The total density of states in a 1-D system is less than that in a three-dimensional (3-D) (or classical) system, especially for low energies. Carriers occupying the lowest energy levels behave like quantized carriers while those lying at higher energies, which are not as tightly confined in the potential well, can behave like classical (3-D) particles with three degrees of freedom. As the surface electric field increases, the system becomes more quantized and more and more carriers become confined in the potential well. The quantum mechanical confinement considerably modifies the carrier distribution in the channel: the maximum of the inversion charge is shifted away from the interface into the Silicon film. Because of the smaller density of states in the 1-D system, the total population of carriers is smaller for the same Fermi level than in the corresponding 3-D (or classical) case. This phenomenon affect the net sheet charge of carriers in the inversion layer, thus requiring a larger gate voltage in order to populate a 1-D inversion layer to have the same number of carriers as the corresponding 3-D system. This leads to an increase of the threshold voltage of the MOSFET, which is an important issue, especially as the power supply voltages drop to lower levels. The gate capacitance and carrier mobility are also modified by quantum effects. These considerations indicate that the wave nature of carriers can no longer be neglected in very thin nanowire MOSFETs and has to be considered in simulation studies.

Modeling and simulation of nanowire MOSFETs devices are currently experiencing a growing interest due to unique capabilities: (i) simulation provides useful insights into device operation since all internal physical quantities that cannot be measured on real devices are available as outputs in simulation; (ii) the predictive capability of simulation studies makes possible the reduction of systematical experimental investigation of these new ultra-scaled devices; (iii) simulation offers the possibility to test hypothetical devices which have not yet been manufactured. Since computers are today considerably cheaper resources, simulation is becoming an indispensable tool for the device engineer, not only for the device optimization, but also for specific studies such as the investigation of new physical phenomena specific to the ultra-short channels (quantum confinement of carriers or short-channel electrostatic effects).

This chapter presents a simulation study of electrostatics and electronic transport in radially-symmetric nanowire MOSFETs by quantum drift-diffusion numerical simulation. In this chapter, we will investigate the operation of radially-symmetric nanowires using a 2-D/3-D Schrödinger/Poisson solver, called Cylmos. Various methods have been suggested to model quantum confinement effects. Among the approaches that are compatible with classical device simulators based on the drift-diffusion approach, the physically most accurate method is to include the Schrödinger equation into the self-consistent computation of the device characteristics [15]. Then, the solver described in this work provides self-consistent solutions of the 2-D Schrödinger equation and the 3-D Poisson equation in cylindrical coordinates [16]-[18], coupled with the drift-diffusion transport equation; this approach is commonly called in literature “quantum drift-diffusion”. Simulated drain current versus gate voltage characteristics have been compared to data obtained from simulation with commercial code with an excellent agreement. Cylmos also provides a lot of additional information and valuable physical insights (such as the 3-D profile of electrostatic potential, classical and quantum carrier densities in the channel, the energy levels and total inversion charge) used to investigate the influence of short-channel and quantum-mechanical effects. Finally, the quantum drift-diffusion code will be used to analyze electrical parameters such as threshold voltage, drain-induced barrier lowering (DIBL) effect and subthreshold swing in radially-symmetric nanowire MOSFETs. The difference between classical (i.e., without quantum confinement) and quantum threshold voltage, as function of different film silicon thicknesses in nanowire MOSFETs will be also discussed.

The chapter is organized as follows. In section 2 we describe in detail the theoretical background of our approach and the simulation code. After the description of the simulated devices, the Schrödinger, Poisson and current continuity solvers are presented. The complete discretization of these equations and the calculation algorithm are explicitly presented. In Section 3 we compare Cylmos to results obtained using a commercial simulator, particularly in terms of output parameters of the Schrödinger module and drain current. Finally, after this validation step, we use Cylmos to analyze and discuss the operation and performance of circular nanowire-based MOSFETs, including the off-state current, threshold voltage and short channel effects.

2. Theoretical background and description of the simulation code

Cylmos is based on the numerical solving of the Poisson-Schrödinger system coupled with the drift-diffusion equation. The Poisson equation is solved on the entire 3-D structure. To solve the Schrödinger equation, the device is divided into parallel vertical slices (y-z plane) (one slice per mesh point in the x direction). The 2-D Schrödinger equation is solved in each slice to obtain the wave functions, the quantum energy levels and the charge density. The solution of the Schrödinger equation in 2-D Cartesian coordinates for the nanowires with square or circular section is straightforward to implement numerically, but it requires very extensive computation time [19]-[22]. To simplify the calculations, in radially-symmetric nanowires it is possible to reduce the size of the two Poisson and Schrödinger equations. By expressing these equations in cylindrical coordinates and using the property of cylindrical symmetry of the structure (which implies a symmetry of the potential and wave function), the Poisson equation becomes a 2 dimensional equation and the Schrödinger equation becomes one dimensional. The Cartesian coordinates (x, y, z) are converted in cylindrical coordinates (x, r, θ), as shown in Fig.1, and considering these later, the structure is symmetrical with respect to the coordinate θ. Therefore, the coordinate θ will be ignored and the Poisson equation is solved in two dimensions, r and x, on the 2-D mesh shown in Fig. 2(a). In the same way, the circular symmetry allows us to simplify the Schrödinger equation, which will be solved on a 1-D mesh along the radial direction r, as shown in Fig. 2(b). The Schrödinger equation will be solved along vertical cut-lines, perpendicular to the x axis, in each mesh-point i of the x-axis. The advantage of this transformation is that, contrary to the 3-D-Poisson/2-D-Schrödinger equations, the 2-D-Poisson/1-D-Schrödinger system is numerically less CPU consuming. Nevertheless, this system of equations is more difficult to implement due to the form of the Laplacian operator [16]. In the following, we detail the equations for a n-channel fully depleted nanowire MOSFET, but similar equations can be derived for p-channel structures.

Figure 1.

Schematic representation of the radially-symmetric nanowire MOS transistor simulated in this work. The main geometrical parameters are also defined.

Figure 2.

Schematics of the mesh used in Cylmos code and definition of mesh points in x-direction and r-direction.

The physical description of electrostatics of the nanowire device requires solving the coupled Poisson equation with the Schrödinger equation. The general Poisson equation is given by:

$$\nabla \left(\epsilon \nabla V\right)=-\rho$$E1

where ρ is the charge density, Vis the electrostatic potential and ε is the material-dependent permittivity. This is a 3-D equation which can be expressed either in Cartesian coordinates or in cylindrical coordinates. In this work we will use the cylindrical coordinates, as it will be shown in the following. The general independent of time Schrödinger equation, expressed in terms of effective mass, is given by:

$$\left\{-\frac{{\hslash }^2}{2}\left({\nabla }_u\cdot \frac{{\nabla }_u}{{\widehat{m}}^{*}}\right)+U(u)\right\}\Psi =E\Psi$$E2

where ℏ is the reduced Planck constant, Ψ is the wave function and E the energy of the particle, $$1/{\widehat{m}}^{*}$$ is the inverse effective mass tensor and U(u) is the potential energy. The potential energy U(u) is related to the electrostatic potential V(u) via the relation:

$$U(u)=-qV(u)+\Delta E_C(u)$$E3

where ΔE_C(u) represents the band offset at the interface between two materials (i.e. the band offset between the conduction band of the semiconductor and the conduction band of the oxide). Similarly to the Poisson equation, the general form of the Schrödinger equation (2) will be expressed in paragraph 2.3 in cylindrical coordinates. These two equations will be solved self-consistently and the solution of this system of equations will be coupled with the courant continuity equation which is written:

$$\nabla J_n=qR+q\frac{\partial n}{\partial t}$$E4

where R is the net rate of generation-recombination and J_nis the electrons current density given by the following equation in the case of the drift-diffusion model:

$$J_n=-q{\mu }_{n}n\nabla \varphi +q{D}_n\nabla n$$E5

where µ_n is where the electron mobility, D_n is the coefficient of thermal diffusion, and n is the electron density. In equation (5) the first term on the right corresponds to the drift component of the total current density and the second term to the diffusion component.

Poisson, Schrödinger and drift-diffusion equations are solved using a finite difference scheme with a uniform mesh on a domain including the channel, the source and drain regions, the gate-oxide layer and the gate electrodes (Fig. 2(a)). The electric field penetration in the source/drain and electron wave function penetration in the gate-oxide can be thus taken into account.

The general flowchart of our code is presented in Fig. 3. In a first step, the quantized energy levels as well as their associated electron wave functions and populations can be obtained from the self-consistent solving of the Schrödinger and Poisson equations. The next step concerns the continuity equation which is numerically solved to ensure a constant drain current along the device channel. The continuity equation gives the numerical solution of the quasi-Fermi level, used further in the calculation of the charge density in the Schrödinger-Poisson system. After reaching the convergence, the drain current density is finally evaluated from the drift-diffusion formalism. In the following, we will explain in detail the three modules corresponding to the solving procedures for the Poisson, Schrödinger and current continuity equations. The discretization of theses equations and the simplifying assumptions used in the calculations, assumptions mainly related to the cylindrical symmetry of the structure, will be systematically explained and justified.

2.1. Description of the simulated devices

The description of the 3-D architecture of the radially-symmetric nanowire MOSFETs considered in the simulation is presented in Fig. 1. The definition of the main geometrical parameters of this structure is also shown in Fig. 1. The structure is symmetric with intrinsic thin silicon film and highly doped source and drain (N_SD = 3×10²⁰ cm^-3). A midgap metal gates (Φ_m=4.61 eV) and a 1.0-nm-thick gate-oxide have been also considered. The source is grounded, and the gate and the drain biased at V_G and V_D, respectively. In order to investigate the influence of short-channel effects and quantum mechanical confinement, different gate lengths and nanowire diameters will be simulated in section 4.

Figure 3.

Flowchart of the numerical code illustrating the different problems to solve in the quantum drift-diffusion calculation of the output drain current.

2.2. Poisson solver

For a given channel geometry and bias conditions, the 2-D Poisson equation is solved (see flowchart in Fig. 3) considering a net charge density given by [23]:

$$\rho =q\left[N_{DOPING}^{}-n\right]$$E6

where q is the absolute value of the electron charge, N_DOPING is the doping atom concentrations in the Si film (N_DOPING = +N_D in source and drain regions and N_DOPING = -N_A in the channel) and n(x, r) is the electron density. This electron density is given by the Schrödinger module (see paragraph 2.3).

The Poisson equations expressed in cylindrical coordinates (x, r, θ) and ignoring the dependence in θ due to the circular symmetry of the electrostatic potential can be written under the form:

$$\frac{{\partial }^{2}V}{\partial x^2}+\frac{1}{r}\frac{\partial V}{\partial r}+\frac{{\partial }^{2}V}{\partial r^2}=-\frac{\rho (x,r)}{\epsilon }$$E7

This equation is discretized using a three-point finite difference scheme applied to the grid shown in Fig. 2, as follows:

$$V_{i,j}=\left[\frac{{\rho }_{i,j}}{{\epsilon }_j}+A_{i,j}+B_{i,j}\right]/{\Sigma }_{i,j}$$E8

where A_i,j, B_i,j and Σ_i,j are given by:

$$A_{i,j}=\frac{V_{i-1,j}+V_{i+1,j}}{\Delta x^2}+V_{i,j-1}\left(\frac{1}{\Delta r^2}-\frac{1}{r_j}\frac{1}{2\Delta r}\right)+V_{i,j+1}\left(\frac{1}{\Delta r^2}+\frac{1}{r_j}\frac{1}{2\Delta r}\right)$$E9

$$B_{i,j}=\frac{{\epsilon }_{j+1}-{\epsilon }_{j-1}}{2\Delta r}\times \frac{V_{i,j+1}-V_{i,j-1}}{2\Delta r}\times \frac{1}{{\epsilon }_j}$$E10

$${\Sigma }_{i,j}=2\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta r^2}\right)$$E11

where the index i identifies the mesh points in x-direction (i is an integer, i=1,..,M), j identifies the mesh points in r-direction (j is an integer, j=1,..,N), Δx and Δr are the mesh sizes between two mesh points in the x-direction and r-direction, respectively (Fig. 2) and r=j‧Δr.Δx and Δr are considered constant in this work. The electrostatic potential calculated by eq. (8) is then used in the solving of the Schrödinger equation for the calculation of the quantum charge density (see below).

2.3. Schrödinger solver

The Schrödinger equation in cylindrical coordinates (x, r, θ) is given by:

$$-\frac{{\hslash }^2}{2}\left[\frac{\partial }{\partial x}\left(\frac{1}{m^{*}}\frac{\partial \Psi }{\partial x}\right)+\frac{\partial }{\partial r}\left(\frac{1}{m^{*}}\frac{\partial \Psi }{\partial r}\right)+\frac{1}{r}\frac{1}{m^{*}}\frac{\partial \Psi }{\partial r}+\frac{1}{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{m^{*}}\frac{\partial \Psi }{\partial \theta }\right)\right]+U\Psi =E\Psi$$E12

where m*(x, r, θ) is a position dependent effective mass, U(x, r, θ) is the potential energy of the conduction-band edge. It is assumed that m* is constant within the silicon nanowire so that it is not dependent on θ. The effective mass is also considered energy-independent (non-parabolicity effects are thus neglected).

Due to the cylindrical symmetry, we may consider wave functions of the form:

$$\Psi (x,r,\theta )=\Psi (x,r)e^{im\theta }$$E13

where m=0, ± 1, ± 2, ± 3,… is the angular momentum quantum number. This wave function is normalized according to:

$$2\pi {\displaystyle \underset{0}{\overset{R}{\int }}r{\left|\Psi (x,r)\right|}^{2}dr}=1$$E14

where the factor 2π is due to the integration with respect to θ.

Using relation (13), equation (12) becomes:

$$-\frac{{\hslash }^2}{2}\left[\frac{\partial }{\partial x}\left(\frac{1}{m^{*}}\frac{\partial \Psi }{\partial x}\right)+\frac{\partial }{\partial r}\left(\frac{1}{m^{*}}\frac{\partial \Psi }{\partial r}\right)+\frac{1}{r}\frac{1}{m^{*}}\frac{\partial \Psi }{\partial r}\right]+\left(U+\frac{{\hslash }^2}{2}\frac{1}{m^{*}}\frac{m^2}{r^2}\right)\Psi =E\Psi$$E15

We recall here that the Schrödinger equation is solved along vertical cut-lines (1-D mesh along the radial direction), perpendicular to the x axis, in each mesh-point i of the x-axis (Fig. 2(b)). The procedure for solving the Schrödinger equation will be the same for each mesh-point i. In the following, we set a fixed mesh point i (then the dependence on x in the equation (14) vanishes); but this procedure will be repeated M times for each mesh-point of the x-axis (this will be explained in detail at the end of this paragraph).

As the Poisson equation, Schrödinger equation (14) is discretized using finite differences scheme to obtain a standard matrix eigen-value problem:

$$H\Psi =E\Psi$$E16

where H is the general Hamiltonian matrix. In practice, equation (16) is solved for each value of the angular momentum quantum number m separately. In the following, we will describe in detail the Hamiltonian H for a given m value.

To simplify the equations, we consider firstly m=0, and we will address the other cases later in this paragraph. With this condition, the discrete form of the Schrödinger equation for the ladder of states with m=0 on a mesh point j along the radial direction is written as:

$$-\frac{{\hslash }^2}{2}\left[\frac{{\Psi }_{j+1}-{\Psi }_j}{m*\Delta r^2}\right]+\frac{{\hslash }^2}{2}\left[\frac{{\Psi }_j-{\Psi }_{j-1}}{m*\Delta r^2}\right]-\frac{{\hslash }^2}{2}\frac{1}{(j-1)\Delta r}\left[\frac{{\Psi }_{j+1}-{\Psi }_{j-1}}{2m*\Delta r}\right]+U_{i,j}{\Psi }_j=E{\Psi }_j$$E17

where the effective mass has been considered independent on r (since we simulate here a fully silicon nanowire). In U_ij, the index i corresponds to the cut-line on the x axis for which the Schrödinger equation is solved. The potential energy term U_ij is obtained using eq. (3) which is discretized as:

$$U_{i,j}=-q{V}_{i,j}+\Delta E_C{}_{i,j}$$E18

where V_i,j is the electrostatic potential obtained from the resolution of the Poisson equation as described in the previous paragraph (Poisson module).

We explain below how the effective masses are calculated in the cylindrical silicon nanowire considered here.The conduction band of silicon has 6 equivalent minima. The variation of the energy of the conduction band near the minimum is not isotropic and constant energy surfaces are ellipsoids around each of the axes (Fig. 4). In each minimum, the electrons propagating along the axis of the ellipsoid are characterized by a longitudinal effective mass, m_l, and those that propagate in a direction perpendicular to this axis are characterized by a transverse effective mass, m_t. In the circular nanowire simulated here, the two valleys along the x direction (transport direction) have a confinement mass equal to m_t and a transport mass equal to m_l. The other four valleys following the directions y and z are anisotropic.

Figure 4.

Schematic illustration of electron confinement in a Silicon nanowire with circular cross-section.

In order to exploit the cylindrical symmetry, we use an approximate isotropic effective mass in the radial direction. Because of this assumption, now there are only two different types of conduction band valleys: the valley along the x direction which is doubly degenerate and the valley along the radial (r) direction whose degeneracy is equal to 4. The effective masses (m_r, m_x) along radial and axial directions for the two types of valleys are given in Table 1.

Valley	Number of equivalent valleys	Radial effective mass, mr (confinement)	Axial effective mass, mx(transport)
1-2	2	mr1 = mt	ml
3-6	4	$$m_{r2}=\frac{2\text{}{m}_t{m}_l}{m_t+m_l}$$	mt

Table 1.

Definition of the effective masses along the radial and axial directions in cylindrical nanowire (m_t=0.19×m₀, m_l=0.98×m₀, m_r2=0.35×m₀, m₀= 9.11×10^-31Kg).

It is important to note that in addition to the solving in one dimension of the Schrödinger equation, which drastically reduces the computation time compared to a 2-D solving, the computation time is further reduced because of the smaller number of valleys taken into account in the calculations compared to the case where an anisotropic effective mass is considered in the radial direction.The choice of an approximate isotropic effective mass for r-direction is justified by the fact that the eigenvectors and eigenenergies resulting from the cylindrical solver are very similar to those of an exact 2-D solver taking into account the anisotropy of the effective mass [24]. In the following, the effective mass m^* used in the Schrödinger equation is equal to the radial effective mass m_r, which may have two distinct values, m_r1 for the valleys 1 and 2 and m_r2 for the valleys 3, 4, 5 and 6, as shown in Table 1.

Let’s go back now to the solving the Schrödinger equation. Using the discrete form given by eq. (17), the Schrödinger equation for m=0 may be set in the form of a matrix equation:

$$H^0{\Psi }^0=E^0{\Psi }^0$$E19

where H⁰(0 stands for m=0 condition) is the N×N tridiagonal hamiltonian matrix of the device (N is the total grid-point number in r-direction):

$$H^0=\left[\begin{array}{ccccc}•& \times & 0& \cdots & 0\\ \circ & •& \times & \ddots & ⋮\\ 0& \circ & •& \times & 0\\ ⋮& \ddots & \circ & •& \times \\ 0& \cdots & 0& \circ & •\end{array}\right]$$E20

where symbols are given by the following equations:

$$\circ =-\frac{{\hslash }^{\text{2}}}{2m_r\Delta r^2}+\frac{{\hslash }^{\text{2}}}{4m_r\Delta r^2(j-1)}$$E21

$$•=\frac{{\hslash }^{\text{2}}}{m_r\Delta r^2}+U_{i,j}$$E22

$$\times =-\frac{{\hslash }^{\text{2}}}{2m_r\Delta r^2}-\frac{{\hslash }^{\text{2}}}{4m_r\Delta r^2(j-1)}$$E23

The matrix H₀ defined in (20) is a non-symmetric tridiagonal matrix and the computation of eigenvalues and eigenvectors of such matrix is numerically not optimized. However, using the method proposed in [16] by Zervos and Feiner, we can convert H⁰ into a symmetric matrix S using the relation:

$$S=T{H}^0$$E24

where T is an (NxN) diagonal matrix whose elements are determined as follows. Since S is a symmetric matrix, it must fill the condition:

$$S_{j,j+1}=S_{j+1,j}$$E25

Equation (24) can be rewritten in a discretized form as:

$$S_{j,j+1}=T_{j,j}{H}_{j,j+1}^0$$E26

Combining equations (25) and (26) we obtain the following relation for matrix T:

$$T_{j,j}{H}_{j,j+1}^0=T_{j+1,j+1}{H}_{j+1,j}^0$$E27

Then we can obtain the elements of the matrix T [16]: the first element T₁₁ is set equal to 1 and the subsequent elements are obtained from (27) as:

$$T_{j+1,j+1}=T_{j,j}{H}_{j,j+1}^0/H_{j+1,j}^0$$E28

Matrix S is then symmetric and using equation (19) we obtain:

$$S{\Psi }^0=T{H}^0{\Psi }^0=E^{0}T{\Psi }^0$$E29

In equation (29) Ψ⁰ and E⁰ are the eigenvector and the eigenvalue of the real matrix H⁰, respectively. Then, equation (29) has to be further transformed for obtaining a standard matrix eigenvalue problem. As shown in [16], we define a diagonal matrix L whose elements are the square root of the corresponding elements of T:

$$T_{j,j}=L_{j,j}^2$$E30

Multiplying eq. (29) from the left by L^-1 we obtain:

$$L^{-1}S{L}^{-1}L{\Psi }^0=L^{-1}LL{H}^0{\Psi }^0=E^0{L}^{-1}LL{\Psi }^0=E^{0}L{\Psi }^0$$E31

We define:

$$H^N=L^{-1}S{L}^{-1}$$E32

and

$${\Psi }^N=L{\Psi }^0$$E33

Then, eq. (29) becomes:

$$H^N{\Psi }^N=E^0{\Psi }^N$$E34

which is the new eigenvalue problem. The new Hamiltonian H^N is a symmetric tridiagonal matrix whose eigenvalues are E⁰ and the eigenvectors are Ψ^N. Then, we solve eq. (34) instead of eq. (19) to find eigenvalues and eigenvectors, and to finally calculate the quantum carrier density. The eigenenergies E⁰ corresponding to the wave functions Ψ^N are computed using a QL algorithm [25] with implicit shifts. Next, we obtain the original wave functions Ψ⁰ from eq. (33) as:

$${\Psi }^0=L^{-1}{\Psi }^N$$E35

So far, we explained the resolution of the Schrödinger equation for m=0, which has simplified the equations. To take into account all values of m, eq. (15) must be solved following the same procedure as that used for m=0. The discretization of equation (15) leads to a matrix relation of the same form as the relation (17):

$$H^m{\Psi }^m=E^m{\Psi }^m$$E36

In this equation, the Hamiltonian H^m is always a tridiagonal matrix which has the same elements “$$$$◦” and “×”, respectively given by eqs. (21) and (23). The only term that changes is the diagonal term, “$$•$$”, which becomes:

$$•=\frac{{\hslash }^{\text{2}}}{m_r\Delta r^2}+\frac{{\hslash }^{\text{2}}}{2m_r}\frac{m^2}{j^2\Delta r^2}+U_{i,j}$$E37

Equation (36) is a new eigenvalue problem to be solved separately for each m to find the ladder states for m=0, ± 1, ± 2, ± 3,… The eigenvalues and eigenvectors are determined using the same procedure as previously used for the matrix H⁰. For each m, a distinct set of wave functions and energy levels is found. The energy levels obtained in this way for all values of m will be sorted and classified in ascending order, from the lowest level and to the highest level. The wave function corresponding to each classified level will be also stored. The first lowest levels (and the corresponding wave functions), in a number equal to N_level, will be used next to calculate the electron density, n(r). In the following, we call E_k (respectively Ψ_k), these first lowest levels (respectively their corresponding wave function), where k is the subband index (k=1,..,N_level) which implicitly includes the angular momentum quantum number m.

The procedure described before for the solving of the Schrödinger equation is performed for both m_r=m_r1 and m_r=m_r2. Two distinct sets of energy levels and wave functions are then obtained and are used in the following for the calculation of the electron density, n(r).

The 1-D electron gas density is given by:

$$n(r)={\displaystyle \sum _{r1,r2}{g}_{r1,r2}{\displaystyle \sum _{k=1}^{N_{level}}{n}_k{\left|{\Psi }_k(r)\right|}^2}}$$E38

where g_r1,r2 is the valley degeneracy (2 for m_r=m_r1 and 4 for m_r=m_r2) and n_k is the electron density per unit of length of the k^th subband given by:

$$n_k=N_{1D}{\Im }_{-1/2}\left(\frac{E_{Fn}-E_k}{k_{B}T}\right)$$E39

where N_1D is the effective density of states (DOS), k_B is the Boltzmann constant, T is the temperature, E_Fn is the quasi-Fermi level which will be calculated in the current continuity module (by a self-consistent coupling of the solution of the Poisson-Schrödinger system to the current continuity equation) and $${\Im }_{-1/2}$$ is the Fermi-Dirac integral of order -1/2 [26]. The 1-D density-of-states N_1D is calculated as:

$$N_{1D}=\frac{1}{\pi }\sqrt{\frac{2m_{r1,r2}^{1D}{k}_{B}T}{{\hslash }^2}}$$E40

where $$m_{r1,r2}^{1D}$$ is the 1-D density of states effective mass ($$m_{r1}^{1D}=m_l$$for valleys 1-2 and $$m_{r2}^{1D}=m_t$$for valleys 3-6). Similar to energy levels and wave functions, N_1D and n_k are calculated for the two types of valleys described in Table 1. Then, these values are used in eq. (38) for the first sum in the right side of the equation. Note that the unity of n_k is m^-1 which multiplied by the normalized probability density $${\left|{\Psi }_k(r)\right|}^2$$(eq. 38) gives the volume electron density n(r).

As we said before, the 1-D Schrödinger equation is solved (using the procedure outlined above) in a vertical cut-line in each mesh point of the x axis (which means that it will be solved M times). For each cut-line in a given mesh-point i, the electron density n(r) is calculated using eq. (38), and it will then be assigned to the mesh-point point i; this makes possible to build step by step the total density n(x, r). This density is injected into the module solving the 2-D Poisson equation (eqs. (6)and (8), paragraph 2.2), which will provide a new potential profile V(x, r). This new potential will be introduced later in the solving of the Schrödinger equation, which will give the new carrier density n(x, r) (in the same manner as described above), which in turn will be injected into Poisson's equation and so on. This procedure allows us to self-consistently solve the two Poisson and Schrödinger equations. The loop stops when the convergence criterion is reached, as shown in Fig. 3.

2.4. Current continuity module

As we explained at the beginning of Section 2, the solution of the Schrödinger-Poisson system of equations is coupled with the current continuity equation given by relation (4). In eq. (4), the current density is defined by a model of charge transport, usually obtained in a semi-classical approximation or simplifying the Boltzmann transport equation (ETB). The simplifying assumptions made on this equation lead to several different semi-classical models such as the drift-diffusion model or the hydrodynamic model. Nevertheless, whatever the model used for the calculation of the carrier density (drift-diffusion, hydrodynamic, etc), the condition of current continuity should be guaranteed. Indeed, the continuity equations describe the evolution of carriers in the silicon nanowire film (source-channel-drain) in order to maintain a constant current along the film in which there is no charge accumulation.

In this work we consider the drift-diffusion model for which the current density is given by equation (5). The charge transport in the nanowire MOSFET simulated here involves only the minority carrier conduction. Therefore, the current density used herein concerns uniquely the electrons (for an n-channel transistor, NMOSFET), but a similar equation can be used for holes in the case of a p-channel transistor (PMOSFET) and the procedure of solving the current continuity equation will be rigorously identical with that shown below for electrons. In addition, to simplify the equations, we consider a stationary regime of operation and we assume that no process of generation or recombination of carriers occurs in our simulations. With these assumptions, equation (4) becomes:

$$\nabla J_n=0$$E41

The current density of electrons can be rewritten as a function of the quasi-Fermi level of electrons in the silicon nanowire, Φ_F, as:

$$J_n=-q{\mu }_{n}n\nabla \varphi +q{D}_n\nabla n=q{\mu }_{n}n\nabla {\Phi }_F$$E42

where the electron diffusion coefficient, D_n, is calculated from the Einstein’s relation (for the non-degenerate case):

$$D_n={\mu }_n\cdot \left(\frac{k_{B}T}{q}\right)$$E43

In this work, the electron mobility µ_n is considered constant everywhere in the silicon nanowire. Substituting equation (42) into equation (41) we obtain:

$$\nabla \left(-{\mu }_{n}n\nabla V+D_n\nabla n\right)=0$$E44

Using the properties of vector analysis and after various algebraic manipulations, equation (44) becomes:

$$\Delta {\Phi }_F+\beta \nabla {\Phi }_F\nabla \left(V-{\Phi }_F\right)=0$$E45

where β=k_BT/q and V is the electrostatic potential which comes from the solution of the coupled solving of the Poisson and Schrödinger equations. Equation (45) can be rewritten:

$$\Delta {\Phi }_F(x,r)=-G(x,r)$$E46

where

$$G=\beta \nabla {\Phi }_F\nabla \left(V-{\Phi }_F\right)$$E47

In cylindrical coordinates, $$\Delta {\Phi }_F(x,r)$$ is given by the following expression:

$$\Delta {\Phi }_F(x,r)=\frac{{\partial }^2{\Phi }_F}{\partial x^2}+\frac{1}{r}\frac{\partial {\Phi }_F}{\partial r}+\frac{{\partial }^2{\Phi }_F}{\partial r^2}$$E48

and G is written as:

$$G=\beta \left(\frac{\partial {\Phi }_F}{\partial x}+\frac{\partial {\Phi }_F}{\partial r}\right)\left[\frac{\partial \left(V-{\Phi }_F\right)}{\partial x}+\frac{\partial \left(V-{\Phi }_F\right)}{\partial r}\right]$$E49

Equation (46) is discretized using a three-point finite difference scheme applied to the mesh shown in Fig. 2, as follows:

$${\Phi }_F{}_{i,j}=\frac{G_{i,j}+C_{i,j}}{{\Sigma }_{i,j}}$$E50

where Σ_i,j is given by equation (11) and G_i,j and C_i,j are given by:

$$G_{i,j}=\beta D_{i,j}{E}_{i,j}$$E51

$$C_{i,j}=\frac{{\Phi }_F{}_{i-1,j}+{\Phi }_F{}_{i+1,j}}{\Delta x^2}+{\Phi }_F{}_{i,j-1}\left(\frac{1}{\Delta r^2}-\frac{1}{r_j}\frac{1}{2\Delta r}\right)+{\Phi }_F{}_{i,j+1}\left(\frac{1}{\Delta r^2}+\frac{1}{r_j}\frac{1}{2\Delta r}\right)$$E52

$$D_{i,j}=\left[\frac{{\Phi }_F{}_{i+1,j}-{\Phi }_F{}_{i-1,j}}{2\Delta x}\cdot \frac{\left(V_{i+1,j}-{\Phi }_F{}_{i+1,j}\right)-\left(V_{i-1,j}-{\Phi }_F{}_{i-1,j}\right)}{2\Delta x}\right]$$E53

$$E_{i,j}=\left[\frac{{\Phi }_F{}_{i,j+1}-{\Phi }_F{}_{i,j+1}}{2\Delta r}\cdot \frac{\left(V_{i,j+1}-{\Phi }_F{}_{i,j+1}\right)-\left(V_{i,j-1}-{\Phi }_F{}_{i,j-1}\right)}{2\Delta r}\right]$$E54

The quasi-Fermi level calculated by eq. (50) is then injected into the system of equations Poisson-Schrödinger, where it is used for the calculation of the carrier density. This density is then used to calculate a new potential, which is then reinserted into the module solving the continuity equation that gives a new quasi-Fermi level. In its turn, this new quasi-Fermi level will be injected into the Schrödinger-Poisson and so on. The loop will stop when the convergence criterion is reached. The final values of the electron density, n(x,r), and the quasi-Fermi level, Φ_F(x,r), will be used in the following for calculating the current density and the final drain current of the device.

2.5. Drain current calculation

The drain current density (per unit area) in the drift-diffusion formalism is finally evaluated from equation (42) which can be written as:

$$J_n(x,r)=q{\mu }_{n}n(x,r)\nabla {\Phi }_F(x,r)$$E55

The total current drain as a function of the applied voltages is then calculated by summing the contribution of each current line in the silicon film.

Finally, the code allows us to store the main internal key-parameters such as the carrier density, quantum energy levels, wave functions, potential, quasi-Fermi level, etc. Some of these parameters will be used for the code validation and for the analysis of the device performances in terms of threshold voltage or short-channel effects.

3. Simulation code validation

Before using the code for analyzing the operation of nanowire MOSFETs, we compared the results of the Schrödinger-Poisson module to the same data obtained from a commercial simulator (Silvaco, [24]) which can simulate cylindrical nanowires. We therefore compare in the following the results issued from Cylmos and Silvacoin terms of quantum energy levels, wave functions and drain current. Each time, we will obtain a very good agreement between the two data sets, which demonstrates the validity of our code.

3.1. Quantum energy levels and wave functions

We begin this comparison by looking in detail the internal parameters of the Schrödinger module, namely the quantum energy levels and associated wave functions. We distinguish here between the two radial effective masses, m_r1 and m_r2. Thus, two sets of energy levels and wave function are compared, one set for each effective mass. The results shown in the following are obtained for a nanowire of diameter D=10 nm (i.e. R=5 nm). The simulated MOSFET is a long channel transistor, with 1-nm-thick gate oxide and an intrinsic channel. Other geometric parameters and doping levels are those described in paragraph 2.1.

Table 2 shows the comparison between the energy levels, E_k, calculated by Cylmos and those issued from Silvaco for the radial effective mass. The transistor is biased at V_G = 0 V and a very low polarization is applied on the drain (V_D=10 mV). As the channel is long and the drain bias is negligible, the potential does not vary much along the axial direction x. Therefore, almost identical energy levels will be obtained in any mesh point i of the intrinsic channel region (outside the source and drain regions). The results shown in Table 2 (as well as those reported in Figs. 5-6, Table 3 and Fig. 7) are calculated for a vertical cut-line in the middle of the channel. Table 2 shows that the energy levels for the radial effective mass calculated by Cylmos agree very well with those obtained using Silvaco.

	E2-E1	E3-E1	E4-E1	E5-E1	E6-E1	E7-E1	E8-E1	E9-E1	E10-E1
Silvaco	0.0406	0.0406	0.094	0.094	0.111	0.158	0.158	0.194	0.194
Cylmos	0.0409	0.0409	0.094	0.094	0.112	0.1601	0.1601	0.199	0.199

Table 2.

Energy levels, E_k, calculated by Cylmos and issued from Silvaco for the radial effective mass m_r1 in a circular nanowire of diameter D=10 nm. The transistor is biased at V_G=0 V and V_D=10 mV. The energy levels are calculated for a vertical cut-line in the middle of the channel.

The wave functions corresponding to the energy levels reported in Table 2 are plotted in Fig. 5 as a function of the position on the radial direction (for k from 1 to 6) and in Fig. 6 (for k from 7 to 10). The wave functions calculated by Cylmos correspond very well to those from Silvaco, particularly in the region corresponding to the silicon (0 to 5 nm). However, a rather clear difference is observed in the oxide region, due to the fact that Silvaco, contrary to Cylmos, does not take into account the penetration of wave functions in the oxide. This phenomenon is included in Cylmos, as previously stated, which leads to this discrepancy between the wave functions of Cylmos and Silvaco in the gate oxide.

The same analysis can be conducted for the energy levels and wave functions calculated for the second radial effective mass, m_r2. These results are shown in Table 3 and Fig. 7 for the same transistor and the same polarization as those used above. A very nice agreement is found again between the data calculated by Cylmos and those obtained from Silvaco.

Figure 5.

Wave function as a function of r calculated by Cylmos and Silvaco for the radial effective mass m_r1 in a nanowire of diameter D=10 nm. (a) k=1,2,3; (b) k=4,5,6. Other parameters and biases are the same as reported in Table 2.

Figure 6.

Wave function as a function of r calculated by Cylmos and Silvaco for the radial effective mass m_r1 in a nanowire of diameter D=10 nm and k from 7 to 10. Other parameters and biases are the same as reported in Table 2.

	E2-E1	E3-E1	E4-E1	E5-E1	E6-E1
Silvaco	0.064	0.064	0.148	0.148	0.175
Cylmos	0.067	0.067	0.154	0.154	0.183

Table 3.

Energy levels, E_k, calculated by Cylmos and those issued from Silvaco for the radial effective mass m_r2 in a nanowire of diameter D=10 nm. Other parameters and biases are the same as reported in Table 2.

In this analysis we studied the energy levels E_k and the corresponding wave functions Ψ_k, with k varying between 1 and 10. We recall here that the variation of k includes the variation of the angular momentum quantum number, m. This means that each value of k corresponds to a value of m, and when k varies from 1 to 10, m can have the following values: 0, ± 1, ± 2,... We also recall that the energy levels are sorted in ascending order and that only the first N_level (10 in the present study) lowest levels are stored. The level k=1 always corresponds to m=0 and all energy levels with non-zero m are doubly degenerate, which corresponds to + m and-m or to the clockwise and counter-clockwise polarization of the wave function. This can be seen, for example in Table 2, where E₂ and E₃ levels are twofold degenerated (as well as the pairs of levels E₄-E₅, E₇-E₈ and E₉-E₁₀).

Another remark concerns the variation of the wave functions with the radial position of the axis. For m=0 (which corresponds to k=1 or k=6 for the case studied here), the wave function is always non zero and has zero first derivative at r=0; this is valid for both radial effective masses. This remark can be verified in the results presented in Figs. 5 and 7. For m=1, the wave function is always zero and has zero second derivative at r=0. For m>1, the wave function is always zero and has zero first derivative at r=0.

Figure 7.

Wave functions as a function of r calculated by Cylmos and Silvaco for the radial effective massm_r2 in a nanowire of diameter D=10 nm. (a) k=1,2,3; (b) k=4,5,6. Other parameters and biases are the same as reported in Table 2.

3.2. Drain current

We also compare our code Cylmos with Silvaco in terms of drain current. Unfortunately, Silvaco does not provide the drain current, because the coupling between the Schrödinger-Poisson system of equations and the continuity equation is not yet implemented. Nevertheless, Silvaco couples the Poisson equation directly to the continuity equation (i.e. the Schrödinger equation is ignored), which makes possible the calculation of the drain current in the so-called "classical" case (this means that the quantum confinement effects are not taken into account in this approach). To facilitate comparison with Silvaco, we have implemented the same procedure in Cylmos. Thus, in our code it is possible to calculate, in addition to the "quantum" drain current (issued from the self-consistently solving of the Poisson, Schrödinger and continuity equations), the "classical" drain current (by directly coupling the Poisson equation to the continuity equation). In this “classical” case, the electron density is reevaluated using the well-known Fermi-Dirac or Boltzmann statistics:

$$n(x,r)=N_C\exp \left[-{\beta }_T\left(E_C(x,r)-{\Psi }_F(x,r)\right)\right]$$E56

where E_C(x,r) is the minimum of the conduction band and N_C is the effective 3-D density-of-states for the conduction band [27].

The solving procedure is the same as explained in Section 2. The discretized Poisson equation (eq. (8)) is solved taking into account this time no more the electron density issued from the Schrödinger equation, but the "classical" electron density given by eq. (56). The potential is then calculated and is afterward injected directly into the module solving the continuity equation. The output of this module is the new quasi-Fermi level which is subsequently used in eq. (56) to calculate a new carrier density to be used in the Poisson module. In this way, a new potential is found and will be injected into the continuity equation and so on. The loop stops when the convergence criterion is reached.

Figure 8.

Drain current versus gate voltage calculated by Cylmos and Silvaco in the “classical” case (i.e., without quantum effects) for a nanowire MOSFET with L=40 nm channel length and a silicon nanowire of diameter D=3 nm.

Figure 8 shows the drain current, simulated in this way, for a long channel MOSFET (L=40 nm) with a circular nanowire with a diameter D=3 nm. The drain current is plotted as a function of the gate voltage for V_D=50 mV. This figure illustrates that a good agreement is obtained between the results of Cylmos (“classical” drain current) and those provided by Silvaco.

Figure 9.

Drain current versus gate voltage calculated by Cylmos and Silvaco in the “classical” case (i.e. without quantum effects) for a short-channel nanowire MOSFET with L=D=8 nm.

We also tested a short-channel transistor, for which the gate length is equal to the diameter of the nanowire (L=D=8 nm). This transistor is probably provided with quite important short-channel electrostatic effects and it is necessary to verify that our code correctly reproduces these effects. In addition, the short-channel effects are considerably enhanced when the drain voltage is high, therefore it is important to test the current characteristics for several V_D. Figure 9 shows such a verification on both characteristics plotted in linear scale (for a better view on the regime above the threshold) and in logarithmic scale (to visualize the subthreshold regime). These figures show once again an excellent agreement between Cylmos and Silvaco, which means that the electrostatic effects are properly taken into account in Cylmos.

4. Simulation results

After validating Cylmos by comparison with a commercial simulator, we use in the following our code to analyze the operation and performances of circular nanowire-based MOSFETs subjected to quantum confinement. We will particularly focus on the variation of the drain current as a function of the nanowire diameter and of the gate length and we will discuss the device performances in terms of off-state current, threshold voltage and parasitic short-channel effects. These key parameters will also be compared between the classical and the quantum-mechanical case.

4.1. Drain current characteristics

We start with the variation of the drain current as a function of the nanowire diameter. Figure 10(a) plots the drain current characteristics as a function of the gate voltage, as calculated by Cylmos, for two diameters D=5 nm and D=8 nm. Data obtained in both classical and quantum cases are shown in this figure. A first remark concerns the off-state current (i.e., value of the drain current for V_G=0 V): for D=5 nm this current is lower than for D=8 nm, whatever the computation method, classical or quantum. Similarly, the current in the subthreshold regime is lower for the smallest diameter. This is due to the so-called “volume inversion” phenomenon which takes place in all multiple-gate devices in the subthreshold regime [28]-[29]. Thus, we logically find this phenomenon here in the circular nanowire, for which the gate completely surrounds the silicon film. Because of the volume inversion, the subthreshold current is directly proportional to the nanowire diameter, which implies that the off-state current decreases when the nanowire is thinned.

A second remark concerns the reduction of the drain current in the quantum case compared to the classical case, for both nanowire diameters considered here. As we explained earlier in this chapter, the existence of a strong quantum confinement of carriers leads to the reduction of the charge inversion in the channel (compared to the classic case where this effect is not taken into account). Therefore, for the same gate voltage V_G, the quantum drain current is necessarily smaller than the classical drain current, and this for all nanowire diameters. This behavior is similar to that of all other multiple-gate devices where strong quantum confinement effects (1-D or 2-D) occur [30], mainly due to the geometrical construction of these devices.

Another point concerns the difference between the current calculated in the classical case and the current calculated in the quantum case. This difference is larger for D=5 nm than for D=8 nm, as can be seen in Fig. 10(a), particularly in the subthreshold regime. The explanation comes from the large variation of quantum confinement effects with the nanowire diameter. When the nanowire is thinned, the carrier confinement is more pronounced (which also implies an increase of the allowed energy levels) and the inversion charge density in the channel reduces. This increases strongly the threshold voltage and leads to a significant decrease of the drain current.

Figure 10.

Drain current characteristics versus the gate voltage calculated by Cylmos: (a) L=40 nm and two nanowire diameters D=8 nm and D=5 nm; (b) D=5 nm and two channel lengths L=40 nm and L=5 nm. The drain current is calculated in both quantum and classical cases.

Finally, Fig. 10(b) shows the variation of the drain current as a function of the gate length for D=5 nm. In this figure we can observe that when the gate length is reduced (for the same nanowire diameter) the off-state current strongly increases and the subthreshold slope is much higher than in the case of a long channel (L=40 nm), reflecting a significant degradation of the device performances. This is due to the huge increase of short-channel effects, because in this ultra-short channel having a gate length equal to the nanowire diameter, the gate control on the electrostatic potential of the channel is much less effective. Thus, parasitic electrostatic effects dramatically increase and lead to this large increase in off-state current and subthreshold slope. This behavior appears for both classical or quantum cases.

4.2. Threshold voltage and short channel effects

For a more detailed analysis of these parasitic short-channel effects, we will now focus on the variations of the threshold voltage as a function of the gate length and of the nanowire diameter. The threshold voltage is extracted using the constant current method [31] from the drain current characteristics as a function of the gate voltage. In Fig. 11 we plot the variation of this threshold voltage as function of the nanowire diameter in the case of a long-channel transistor and in the case of a transistor with a channel length equal to the nanowire diameter (L=D). The results reported in this figure come from a quantum-mechanical calculation. We note that the threshold voltage of the long-channel transistor is higher than that of the short-channel transistor with L=D. This is an expected result because the short channel effects that occur in the transistor with L=D greatly reduce the threshold voltage. One can also note that when the diameter of the nanowire is reduced, the threshold voltage remains constant over a certain range of values and next begins to increase with decreasing the nanowire diameter (in both cases, long channel and short channel with L=D). This increase is due, of course, to the fact that the quantum confinement of carriers is more significant when the diameter is reduced, which induces an increase in the threshold voltage (as explained earlier). This increase is higher in the case of long channel than in the case of a short-channel. This indicates that the strong short-channel effects that exist in the channel with L=D affect the quantum confinement, which becomes weaker than in a long-channel transistor.

Another way to highlight this phenomenon is to calculate the impact of short channel effects, reflected in the SCE metric, which is equal to the difference between the threshold voltage of the short-channel transistor with L=D and the threshold voltage of the long-channel transistor (for a given diameter D). SCE is then plotted as a function of the nanowire diameter in Fig. 12. This figure shows that, after an increase and a maximum at D=8 nm, the SCE starts to decrease as the diameter reduces below 8 nm. This behavior is explained by the fact that the strong quantum confinement that occurs in these ultrafine nanowires reduces short-channel effects, which will probably lead to increased performance of these ultra-short and ultra-thin devices.

Figure 11.

Threshold voltage versus the nanowire diameter extracted from quantum calculation of the drain current characteristics.

Figure 12.

SCE metric versus the nanowire diameter for short-channel devices with L=D.

5. Conclusion

In conclusion, we described in this chapter a numerical solver, called Cylmos, that provides self-consistent solving of the 2-D Schrödinger equation and the 3-D Poisson equation in cylindrical coordinates, coupled with the drift-diffusion transport equation. We presented in details the discretization of these equations and the overall solving algorithm used to obtain the final drain current of the device. We also paid a particular attention to justify the approximations and simplifications used in the solving of these equations in cylindrical coordinates. A meticulous validation step was performed by comparing the results of our code with those issued from a commercial code. We demonstrated that our code properly takes into account particular phenomena specific to MOSFET devices based on circular nanowires with surrounding gate, such as volume inversion, short channel and quantum-mechanical effects. Simulated drain current versus gate voltage characteristics have also been successfully compared to data obtained from simulation with this commercial code. Our code was finally used to analyze electrical parameters such as drain current, off-state current, threshold voltage and short-channel effects in ultra-short radially-symmetric nanowire MOSFETs. We have shown that quantum-mechanical confinement is very important in circular nanowires, but the presence of this phenomenon reduces the impact of parasitic short-channel effects. This is a very encouraging result for the operation of future integrated circuits based on ultra-short and ultra-thin circular nanowires.