scieee AI-readable full text Open interactive document viewer

Multi-Subband Ensemble Monte Carlo simulations of scaled GAA MOSFETs

Donetti, Luca,Sampedro Matarín, Carlos,García Ruiz, Francisco Javier,Godoy Medina, Andrés,Gámiz Pérez, Francisco Jesús

Abstract

We developed a Multi-Subband Ensemble Monte Carlo simulator for non-planar devices, taking into account two-dimensional quantum confinement. It couples self-consistently the solution of the 3D Poisson equation, the 2D Schrödinger equation, and the 1D Boltzmann transport equation with the Ensemble Monte Carlo method. This simulator was employed to study MOS devices based on ultra-scaled Gate-All-Around Si nanowires with diameters in the range from 4 nm to 8 nm with gate length from 8 nm to 14 nm. We studied the output and transfer characteristics, interpreting the behavior in the sub-threshold region and in the ON state in terms of the spatial charge distribution and the mobility computed with the same simulator. We analyzed the results, highlighting the contribution of different valleys and subbands and the effect of the gate bias on the energy and velocity profiles. Finally the scaling behavior was studied, showing that only the devices with D = 4 nm maintain a good control of the short channel effects down to the gate length of 8 nm .

Full text

Multi-Subband Ensemble Monte Carlo simulations of scaled GAA MOSFETs L. Donetti∗, C. Sampedro, F.G. Ruiz, A. Godoy, F. Gamiz Departamento de Electr´onica and CITIC, Universidad de Granada, 18071 Granada, Spain. Abstract We developed a Multi-Subband Ensemble Monte Carlo simulator for non-planar devices, taking into account two-dimensional quantum confinement. It couples self-consistently the solution of the 3D Poisson equation, the 2D Schr¨ odinger equation, and the 1D Boltzmann transport equation with the Ensemble Monte Carlo method. This simulator was employed to study MOS devices based on ultra-scaled Gate-All-Around Si nanowires with diameters in the range from 4 nm to 8 nm with gate length from 8 nm to 14 nm. We studied the output and transfer characteristics, interpreting the behavior in the sub-threshold region and in the ON state in terms of the spatial charge distribution and the mobility computed with the same simulator. We analyzed the results, highlighting the contribution of different valleys and subbands and the effect of the gate bias on the energy and velocity profiles. Finally the scaling behavior was studied, showing that only the devices with D=4 nm maintain a good control of the short channel effects down tho the gate length of 8 nm. Keywords: Gate-all-around MOSFET, nanowire, Monte Carlo simulation, multi-subband, short-channel effects. 1. Introduction Non-planar MOS transistors with multiple gates represent the most promising solution to the ultimate scaling of CMOS technology [1] because of their superior immunity to short channel effects. Indeed, they are not only a possible future option under scrutiny, but they have already entered mass production [2, 3]. When the size of a device reaches the nanometric scale, quantum effects play an important role; quantum confinement is relevant for reduced lateral size (i. e. in the direction perpendicular to transport) and even more so in nonplanar structures where confinement is two-dimensional. An approximately correct electron distribution taking into account quantum confinement can be obtained with corrections to the classical potential, employing algorithms which, however, need proper calibration for each considered structure and crystal orientation (see for example [4, 5]). As a consequence, to properly take into account quantum confinement effects in the cross section of the device it is necessary to solve the Schr¨ odinger Equation (SE). The Multi-Subband Ensemble Monte Carlo (MSEMC) approach employs the semi-classical Monte Carlo (MC) method to solve the Boltzmann transport equation, coupled with the solution of the SE in the perpendicular direction. The MC method, compared to common transport framework such as DriftDiffusion, takes directly into account non-equilibrium carrier transport. On the other hand, compared to full quantum approaches, MC allows a simpler implementation of scattering mechanisms and a relatively reduced computational effort. In this work, we describe an MS-EMC simulator for 3D devices developed by our group, and show the possibilities it ∗Corresponding author Email address: [email protected] (L. Donetti) offers by studying the scaling properties of ultimate Gate-AllAround (GAA) MOS transistors with ultra-thin Si nanowires. Section 2 is devoted to the description of the simulator, while in the following one (Section 3) we define devices under study and report the obtained results and insights into the device behavior. Finally, conclusions are drawn in Section 4 2. Simulator description To model a non-planar device, such as GAA MOSFETs, we employ a simulator based on the MS-EMC approach. This method has been widely and successfully employed for the simulation of planar semiconductor devices [6, 7, 8, 9], and only recently it has been applied to 3D devices [10, 11]. The simulator is based on the space-mode approach [12], where the SE is solved in several cross sections perpendicular to the transport direction z, for each considered conduction band valley. In this way, the eigen-energies, Eν,i(z), and the wave functions, ξν,i(x,y,z), are obtained for different values of zalong the device, where νand iare the valley and subband indices, respectively. After that, the energy levels are corrected for the effects of non-parabolic band structure as described in [13] while wave functions are left uncorrected. Carrier transport is simulated through the MC method: once electrons are assigned to a subband (specified by its valley νand index i), its motion is restricted to one dimension in the zdirection and the driving force is computed employing the derivative of subband energy levels, that is F=−∂Eν,i(z)/∂z. Then, as usual in the MC scheme, the subband population, nν,i(z) is computed by counting the simulation particles in subband ν, i at different cross sections. By multiplying the population by the corresponding distribution function |ξν,i(x,y,z)|2, the total Preprint submitted to Solid-State Electronics August 19, 2017 ‚Fix VG;VD“0 ‚Self-consistent loop: §Poisson equation §Schr¨odinger eq. §Fermi distribution Self-consistent Schr¨ odinger-Poisson (SP) ‚Init MC particles ‚Self-consistent loop: §Poisson equation §Schr¨odinger eq. §Monte-Carlo Self-consistent equilibrium MC (SC0) ‚Loop: §Increase VD §Poisson equation §Schr¨odinger eq. §Monte-Carlo Monte Carlo with increasing VD (VD) ‚Self-consistent loop: §Poisson equation §Schr¨odinger eq. §Monte-Carlo Self-consistent Monte Carlo (SC1) ‚Self-consistent loop: §Poisson equation §Schr¨odinger eq. §Monte-Carlo §Compute current Current computation w. self-cons. MC (SC2) converged converged requested VDconverged current converged, new VD Figure 1: Simulation scheme for a given value of VG. electron density is obtained. The resulting charge density is employed in the 3D Poisson Equation (PE), which is solved in a loop with the SE and the MC simulation in order to achieve a self-consistent solution. To improve the stability and the convergence properties of the self-consistent loop, the non-linear PE is employed [14]. For the simulation, the device structure is described employing a 3D finite element mesh with tetrahedral elements. Because of the need to discretize the SE in different cross sections, this 3D mesh is constructed by extruding a 2D triangular mesh. The finite element mesh allows a good representation of complex geometries (e.g. round nanowires, rounded corners, leaning sidewalls in FinFETs) and a natural formulation of the equations (PE, SE) near material boundaries. The MC simulation includes carrier scattering by acoustic and optical phonons [15], taking into account Pauli exclusion principle [16]. To improve the MC statistics, especially in the sub-threshold regime, we employ a variance reduction technique based on non-uniform super-particle weight [17]. Here the weight is computed according to the total energy of the particle injected into the device through the contacts, as in the following equation: w(E)=w0 1+exp qE−EF 2kBT!−1 (1) where EFis the Fermi energy (known at source and drain contacts where electrons are injected), kBis the Boltzmann constant, Tis absolute temperature, and w0a normalization constant, chosen and updated during the simulation in order to obtain a given number of (super-)particles. Notice that the expression of Equation (1) is similar to the Fermi-Dirac distribution function with an added factor 1/2 in the exponent and a normalization factor. The general simulation strategy for a given value of the gate bias VGis shown in the diagram of Figure 1. The algorithm starts with VD=0 V and a self-consistent solution of SE and PE in equilibrium (without transport) is found, employing a predictor-corrector method [18]. This allows the initialization of the MC simulator employing the electron density given by the equilibrium Fermi-Dirac distribution for each subband. Then a loop with the MC simulation and the solution of PE and SE is started. In a first phase (SC0), the drain bias VDis kept equal to 0 V while a self-consistent solution is obtained. Then, a second phase is started (VD), in which the boundary conditions at the drain are changed increasing the voltage by a small ∆VDstep in each iteration of the loop, until the desired value of VDis reached. After that, the boundary conditions are kept fixed and the MC-PE-SE loop (SC1) is repeated until the variation of the potential between one iteration and the following becomes lower than the prescribed tolerance, that, is when a self-consistent solution is obtained. Once this is achieved, the simulator starts computing the current (SC2) and an estimation of the corresponding statistical error, until the latter is below a specified threshold or a maximum number of iterations is reached. Then, the procedure is repeated from the VD stage with the following target value of the drain voltage. Alternatively, it is possible to compute the low-field mobility, µn, by a similar procedure. Only the channel of the device is considered and small values of a uniform electric field are applied in the longitudinal direction. As in the full device simulation, an equilibrium self-consistent state with VD=0 V is first obtained. Then, different small values of the drift electric field, F, are applied by shifting the subband energy levels, while an infinite channel length is emulated by employing periodic boundary condition. Finally, mobility is extracted by fitting the obtained values of the velocity, vnwith the linear relation vn=µnF. To improve the performance of the simulator, a high level of parallelism is employed. In the MC code, the independent super-particle flights are simulated in a parallel fashion through the use of OpenMP; particular care is needed for synchronization to inject and remove particles at contacts and to gather particle statistics. On the other hand, specialized and optimized sparse matrix parallel routines are employed for the solution of the PE and the SE. 3. Results We simulate GAA field-effect transistors based on cylindrical Si nanowires with channel along the h100idirection, with diameter Dranging from 4 nm to 8 nm. The coordinate reference system is chosen so that the xand yaxes are in the cross section plane and the zaxis is in the transport direction. The gate oxide (SiO2) thickness is Tox =1 nm. The channel is considered undoped (NA=1×1013 cm−3) and a midgap metal is assumed for the gate (with work function 4.56 eV). 2 0 0.2 0.4 0 5 10 VG= 0.5 V VG= 0.6 V VG= 0.7 V VD(V) ID(µA) Figure 2: ID-VDcurves for the device with D=8 nm and LG=14 nm. The source and drain doping density is NS D =1×1020 cm−3, with an underlap of Lsp =2 nm and Gaussian distribution with σ=0.8 nm. The simulated devices have gate length LGranging from 14 nm down to 8 nm, while the total length of the source and drain regions, including extensions, is LS D =14 nm each. As stated in Section 2, we employ a 3D mesh obtained by extrusion of a 2D mesh, as the latter is used to discretize the SE in the cross sections. This 2D mesh makes use of a variable number of triangles depending on the nanowire diameter: the minimum is ∼1500 for D=4 nm and the maximum is ∼3500 for D=8 nm. The mesh spacing in the transport direction is finer in the channel region (∆zmin =0.5 nm) and gradually grows in the source and drain regions (up to ∆zmax =2 nm). The number of nodes and tetrahedra in the full 3D mesh is approximately 32 500 and 182 500, respectively, for the smallest considered device; 98 000 and 560 000 for the largest one. All simulations have been performed at T=300 K. Figure 2 shows the simulated output characteristics of the largest device, with diameter D=8 nm, and gate length LG= 14 nm. Figure 3 shows the simulated transfer characteristics for all the considered devices: as expected the current is larger for nanowires with wider cross section and for shorter gate length. The statistical noise inherent to the MC method is greatly suppressed, thanks to the variance-reduction technique employed in the simulator: current fluctuations can be observed only for IDvalues smaller than 1 nA. In any case the residual statistical noise in the ID-VGcurve gives rise to larger fluctuations in the second derivative of IDwith respect to VG: therefore we cannot locate the threshold voltage Vth with the maximum of ∂2ID/∂V2 G. Therefore, we extract Vth by employing a fixed value of ID. To take into account the difference in the gate length of the considered devices, we employ ID,th(LG=14 nm) = 0.08 µA for the longest device and scale it inversely with LG; that is ID,th(LG)=0.08 µA×(14 nm/LG). The results are shown in Figure 4. The obtained values of Vth are larger for D=4 nm for all cannel lengths, due to quantum confinement. Moreover, a strong Vth roll-offis observed for the wider nanowires with D=6 nm or D=8 nm, while only a small drift of Vth is present in the case of D=4 nm. If we restrict our analysys to the longest devices (i.e. with 0.2 0.4 0.6 10−4 10−3 10−2 10−1 100 101 8 nm 6 nm 4 nm VD= 0.05 V VG(V) ID(µA) LG= 8 nm LG= 10 nm LG= 12 nm LG= 14 nm Figure 3: ID-VGcurves for all simulated devices. 8 10 12 14 0.3 0.4 0.5 D= 4 nm D= 6 nm D= 8 nm LG(nm) Vth (V) Figure 4: Threshold voltage Vth (computed at VD=50 mV) as a function of gate length LG, for different values of nanowire diameter D. LG=14 nm) and plot the drain current vs. the gate voltage overdrive, VG−Vth, the curves corresponding to devices with different diameters collapse (see Figure 5) except for very large applied biases. This means that the behavior of IDbelow the threshold voltage and slightly above it is independent of the diameter. To investigate this fact, we plot in Figure 6 the linear electron density ninv, that is the charge density integrated in a cross section near the middle of the channel, versus the gate voltage overdrive. Only very small differences can be observed for values of VG−Vth larger than approximately 0.2 V, indicating that linear density does not depend on the nanowire diameter D. This, in turn, means that the average spatial charge density must be inversely proportional to the cross-section area, for the same gate overdrive. Such behavior can be explained by considering the spatial charge distribution, n, along a cross section in the middle of the gated region as it is shown in Figure 7. In narrow nanowires, such as those considered in this paper, the charge concentration is almost always peaked around the center of the nanowire, as shown in Figure 7 (a), (b), and (c): for VG=Vth +0.2 V the shape of the electron distribution is very similar, with increasing width and decreasing peak concentration as the diameter grows from 4 nm to 8 nm. Only for the widest device (D=8 nm) and larger gate bias (VG=Vth+0.3 V) the maximum value of nis located at a position different from 3 −0.2 0 0.2 0.4 10−3 10−2 10−1 100 101 VD= 0.05 V LG= 14 nm VG−Vth (V) ID(µA) D= 4 nm D= 6 nm D= 8 nm 0 2 4 Figure 5: Drain current IDas a function of gate overdrive VG−Vth for devices with LG=14 nm. −0.2 0 0.2 0.4 101 102 103 104 105 106 107 LG= 14 nm VD= 0.05 V VG−Vth (V) ninv (cm−1) D= 4 nm D= 6 nm D= 8 nm 0 2·106 4·106 6·106 Figure 6: Linear electron density ninv in the middle of the channel as a function of gate overdrive. the geometrical center of the nanowire. In particular, due to the anisotropy of the Si conduction band valleys, four maxima appear along the xand yaxes. If we compare Figures 5 and 6, however, we can notice that the differences among the IDcurves for a gate overdrive larger than 0.1 V are larger than the differences in the ninv curves. Therefore, the differences in the drain current cannot be attributed exclusively to the inversion charge and must stem from other causes. To shed light on this issue we also computed the mobility, µn, shown in Figure 8. Here, we can observe that µnis severely degraded for decreasing D, which justify the residual differences in the IDcurves. The simulator also allows for measuring and comparing internal quantities such as, for example, those related to the subband profiles and populations. In a h100inanowire with circular cross section, the six equivalent ∆valleys of Si conduction band split into nonequivalent sets. With our choice of coordinates, the doubly degenerate ∆xand ∆yvalleys, present their longitudinal direction in the cross-section plane. For both valleys, the confinement masses are ml=0.926 m0in one direction and mt=0.19 m0in the perpendicular one (where m0is the free electron mass): due to the symmetry of the circular shape, −4−2 0 2 4 −4 −2 0 2 4 (a) 02·1019 −4−2 0 2 4 −4 −2 0 2 4 (b) 01·1019 −4−2 0 2 4 −4 −2 0 2 4 (c) 03·1018 6·1018 −4−2 0 2 4 −4 −2 0 2 4 (d) 05·1018 1·1019 Figure 7: Electron density in the middle of the channel for VG=Vth +0.2 V and D=4 nm (a), D=6 nm (b), D=8 nm (c); VG=Vth +0.3 V and D=8 nm (d). The dashed line indicates the Si/SiO2interface. Units of dimensions in the cross sections are nm while the units of electron density are cm−3. the corresponding energy levels are equal, with wave functions with different symmetry axes. On the other hand, the ∆zvalleys with their longitudinal effective mass along the transport direction give rise to subbands with higher energy: in this case the confinement effective mass is isotropic in the cross-section plane and equal to mt. Figure 9 shows the energy profiles of the two lowest subbands of each valley for the largest device, with D=8 nm and LG=14 nm. Notice that, as explained before, the lowest energy subband corresponds to the ∆xand ∆y valleys. However, the confining potential is different inside the channel and in the source and drain regions, so that the subband with second lower energy belongs to different valleys according to the position z. Moreover, since in the ∆zvalleys the transport mass is larger than in the ∆xand ∆yvalleys (mlvs. mt), the population of the corresponding subbands can be larger even if the energy levels are higher, as shown in Figure 10. Here we can see that the population of the first subbands of ∆x,∆yand ∆zare approximately equal in the gated region. We now turn to the analysis of the profiles as a function of the gate voltage VG, in the same device (D=8 nm, LG= 14 nm). Figure 11 shows the fundamental subband energy, Ex,1 (or, equivalently, Ey,1), the total linear electron density, ninv, and the average electron velocity in the transport direction, hvzi, as a function of the position zalong the channel, for VD= 0.05 V and different values of VG. In the subband profile (Fig4 01·1072·1073·107 300 400 500 600 D= 4 nm D= 6 nm D= 8 nm ninv (cm−1) µn(cm2/V s) Figure 8: Phonon limited electron mobility for Si GAA nanowires as a function of inversion charge ninv for different values of nanowire diameter D. −20 −10 0 10 20 −0.1 0 0.1 0.2 VD= 0.05 V VG= 0.4 V z(nm) Eν,i (eV) Ex(y),1Ez,1 Ex(y),2Ez,2 Figure 9: Energy profile of the first two subbands of each valley for the device with D=8 nm and LG=14 nm. Subbands of ∆xand ∆yvalleys have the same energy due to the symmetry of the cross-section. −20 −10 0 10 20 103 105 107 ninv VD= 0.05 V VG= 0.4 V z(nm) nν,i (cm−1) Ex(y),1 Ex(y),2 Ez,1 Ez,2 Figure 10: Population of the first two subbands of each valley for the device with D=8 nm and LG=14 nm. Subbands of ∆xand ∆yvalleys have the same population due to the symmetry of the cross-section. The black solid curve indicates the total electron density. −0.1 0 0.1 0.2 VG= 0.3 V to 0.7 V in 0.1 V steps VD= 0.05 V (a) Ex(y),1(eV) 104 106 108 VG (b) ninv (cm−1) −20 −10 0 10 20 0 0.2 0.4 0.6 0.8 VG (c) z(nm) hvzi(107cm/s) Figure 11: Fundamental subband profile (a), electron density (b), and average electron velocity in the transport direction (c) as a function of position along the device, for VD=0.05 V and different values of VG. All curves correspond to the device with D=8 nm, LG=14 nm. ure 11(a)), we can see that the maximum energy (that is the peak of the source-to-drain barrier) is located in the center of the device in the sub-threshold regime (for VG=0.3 V) and gradually moves towards the source end of the channel as VG increases. This displacement is also associated with a change in the shape of the profile inside the channel: rounded for small VGand almost straight for high VG, sloped towards the drain. In the ninv curve (Figure 11(b)), we can notice that, even in subthreshold regime (when the electron density in the channel is four orders of magnitude smaller than in the source and drain), the profile is very smooth, and the statistical noise produced by the MC method is kept under control thanks to the aforementioned variance reduction technique. Next, the shape of the curves presents a trend similar to the one observed in the subband profile, in the sense that the shape gets more flat when the gate bias increases. However, there is a noticeable difference: the position of the minima in the sub-threshold regime does not correspond to the center of the channel but it is displaced towards the drain. The position of such minima corresponds to maxima of the velocity curves Figure 11(c), as the current is constant in the whole device in a stationary state. As expected, in the channel region the average velocity is higher than in the source/drain regions, because of the reduced electron density. In general, the average velocity increases for larger values of VG, following the increase of the drain current. However, this trend is reversed in a small region inside the channel, around the position of the average velocity maximum: the peak velocity decreases for increasing VG. We can also check the correct implementation of the scattering mechanisms and of the Pauli exclusion principle as mentioned in section 2, by analyzing the carrier distribution as a 5 0 0.2 0.4 10−6 10−3 100 Fermi-Dirac Boltzmann VG= 0.4 V VD= 0 V E−EF(eV) fx,1(kz, z) z= 21 nm z= 8 nm z= 0 nm Figure 12: Occupation function for electrons in the lowest energy subband as a function of energy at different zpositions, for the device with D=8 nm and LG=14 nm. Boltzmann and Fermi-Dirac distribution are also plotted for comparison. 8 10 12 14 0 50 100 150 200 D= 4 nm D= 6 nm D= 8 nm LG(nm) DIBL (mV/V) Figure 13: DIBL, computed employing VD=0.05 V and VD=0.5 V, as a function of the gate length for different values of nanowire diameter. function of energy. To be able to compare to an analytical expression we consider the equilibrium case, with VD=0 V. In this case, there exists a constant Fermi level, EF, in the whole device, aligned with the Fermi level of source and drain, and the energy distribution of electrons must follow the FermiDirac distribution. To implement the Pauli exclusion principle, the electron distribution function is computed by particle counting for every subband and every device cross-section, obtaining the function fν,i(kz,z). Then, to perform a comparison with the Fermi-Dirac distribution we compute fν,i(E,z)= fν,ikz(E),z+fν,i−kz(E),z, where kz(E) is the positive value of wave-vector kzcorresponding to the total energy Efor subband v,i. In Figure 12, we show fx,1(E,z), the distribution function of the fundamental subband (the first subband of ∆xvalley) as a function of energy at three different positions: in the center of the channel (z=0 nm), just outside the channel towards the drain region (z=8 nm), and in the drain region near the contact (z=21 nm). The result show that Fermi-Dirac distribution is correctly reproduced. Finally, we turn to the analysis of the short channel effects (SCEs) in the simulated devices and their dependence on the 8 10 12 14 60 80 100 120 D= 4 nm D= 6 nm D= 8 nm VD= 0.05 V LG(nm) subthreshold swing (mV/dec.) Figure 14: Sub-threshold swing as a function of gate length for different values of nanowire diameter. gate length LG. We can expect the narrower devices to possess better electrostatic control of the channel, especially for larger gate lengths. However, the MS-EMC simulator allows us to obtain quantitative results which properly take into account lateral quantum confinement. A first indication of the SCEs is given by the Drain Induced Barrier Lowering (DIBL), a measure of the variation of the threshold voltage shift caused by the drain bias: DIBL =Vth@VD2 −Vth@VD1/(VD2−VD1). In this case, the DIBL values shown in Figure 13 are computed employing the following drain voltage values: VD1=0.05 V and VD2=0.5 V. In the figure, we can observe that for D=4 nm the DIBL hardly increases for decreasing LG, with a largest value as small as 58 mV/V for the shortest device (LG=8 nm). For the nanowires with with D=6 nm and D=8 nm, gate length scaling produces a larger increase of the DIBL, which exceeds 100 mV/V at LG=8 nm and LG=10 nm, respectively. Similar conclusions can be drawn by observing the subthreshold swing (SS), represented in Figure 14. The longest devices, with LG=14 nm, show the same SS value for every D(within the numerical accuracy): SS ≃65 mV/dec., which is close to the ideal MOSFET limit at room temperature. However, the scaling behavior is quite different: while the narrowest device keeps SS values lower than 70 mV/dec.for all the considered gate lengths, the device with D=8 nm shows a linear increase reaching SS ≃110 mV/dec.at LG=8 nm. 4. Conclusion This paper describes a numerical simulator that solves in a self-consistent way the 3D Poisson equation, the 2D Schr¨ odinger equation, and the 1D Boltzmann transport equation through the Ensemble Monte Carlo method. This software can be used to study the static characteristics of ultra-scaled MOSFET devices, also in the sub-threshold region thanks to the implementation of a variance reduction technique based on variable superparticle weights. It also allows us to inspect the carrier distribution, taking into account quantum confinement effects, and to compare the phonon limited mobility and drain current including the short-channel and high-field effects, both computed 6 with the MC method. We applied this simulator to narrow Si nanowire MOSFETs with reduced gate length. We correlated the differences in the curves above threshold with the different shape of the charge distribution in the cross section of the device and with the mobility degradation for small diameters. We analyzed the electron density and velocity distribution inside the channel of the device as a function of the applied gate bias. Finally we studied the short channel effects and we found that the narrowest cylindrical devices with D=4 nm can keep a good electrostatic control down to LG=8 nm. Acknowledgment The authors would like to thank the financial support of Spanish Government under project TEC2014-59730-R and EU H2020 program under projects REMINDER (grant agreement No 687931) and WAYTOGO FAST (ECSEL-2014-2-662175). References [1] J.-P. Colinge, FinFETs and Other Multi-Gate Transistors, 1st Edition, Springer, 2007. [2] C. Auth, C. Allen, A. Blattner, D. Bergstrom, M. Brazier, M. Bost, M. Buehler, V. Chikarmane, T. Ghani, T. Glassman, R. Grover, W. Han, D. Hanken, M. Hattendorf, P. Hentges, R. Heussner, J. Hicks, D. Ingerly, P. Jain, S. Jaloviar, R. James, D. Jones, J. Jopling, S. Joshi, C. Kenyon, H. Liu, R. McFadden, B. McIntyre, J. Neirynck, C. Parker, L. Pipes, I. Post, S. Pradhan, M. Prince, S. Ramey, T. Reynolds, J. Roesler, J. Sandford, J. Seiple, P. Smith, C. Thomas, D. Towner, T. Troeger, C. Weber, P. Yashar, K. Zawadzki, K. Mistry, A 22nm high performance and low-power CMOS technology featuring fully-depleted trigate transistors, self-aligned contacts and high density MIM capacitors, in: 2012 Symposium on VLSI Technology (VLSIT), 2012, pp. 131–132. doi:10.1109/VLSIT.2012.6242496. [3] S. Natarajan, M. Agostinelli, S. Akbar, M. Bost, A. Bowonder, V. Chikarmane, S. Chouksey, A. Dasgupta, K. Fischer, Q. Fu, T. Ghani, M. Giles, S. Govindaraju, R. Grover, W. Han, D. Hanken, E. Haralson, M. Haran, M. Heckscher, R. Heussner, P. Jain, R. James, R. Jhaveri, I. Jin, H. Kam, E. Karl, C. Kenyon, M. Liu, Y. Luo, R. Mehandru, S. Morarka, L. Neiberg, P. Packan, A. Paliwal, C. Parker, P. Patel, R. Patel, C. Pelto, L. Pipes, P. Plekhanov, M. Prince, S. Rajamani, J. Sandford, B. Sell, S. Sivakumar, P. Smith, B. Song, K. Tone, T. Troeger, J. Wiedemer, M. Yang, K. Zhang, A 14nm logic technology featuring 2nd-generation FinFET, air-gapped interconnects, self-aligned double patterning and a 0.0588 µm2 SRAM cell size, in: 2014 IEEE International Electron Devices Meeting, 2014, pp. 3.7.1–3.7.3. doi:10.1109/IEDM.2014.7046976. [4] F. M. Bufler, L. Smith, 3D Monte Carlo simulation of FinFET and FDSOI devices with accurate quantum correction, Journal of Computational Electronics 12 (4) (2013) 651–657. doi:10.1007/s10825-013-0518-z. [5] M. A. Elmessary, D. Nagy, M. Aldegunde, J. Lindberg, W. G. Dettmer, D. Per´ ıc, A. J. Garc´ ıa-Loureiro, K. Kalna, Anisotropic Quantum Corrections for 3-D Finite-Element Monte Carlo Simulations of Nanoscale Multigate Transistors, IEEE Transactions on Electron Devices 63 (3) (2016) 933–939. doi:10.1109/TED.2016.2519822. [6] J. Saint-Martin, A. Bournel, F. Monsef, C. Chassat, P. Dollfus, Multi subband Monte Carlo simulation of an ultra-thin double gate MOSFET with 2d electron gas, Semiconductor Science and Technology 21 (4) (2006) L29. doi:10.1088/0268-1242/21/4/L01. [7] E. Sangiorgi, P. Palestri, D. Esseni, C. Fiegna, L. Selmi, The Monte Carlo approach to transport modeling in deca-nanometer MOSFETs, Solid-State Electronics 52 (9) (2008) 1414–1423. doi:10.1016/j.sse.2008.04.007. [8] C. Sampedro, F. G´ amiz, A. Godoy, R. Val´ ın, A. Garc´ ıa-Loureiro, F. G. Ruiz, Multi-Subband Monte Carlo study of device orientation effects in ultra-short channel DGSOI, Solid-State Electronics 54 (2) (2010) 131– 136. doi:10.1016/j.sse.2009.12.007. [9] C. Sampedro, F. G´ amiz, A. Godoy, On the extension of ET-FDSOI roadmap for 22 nm node and beyond, Solid-State Electronics 90 (2013) 23–27. doi:10.1016/j.sse.2013.02.057. [10] C. Sampedro, L. Donetti, F. Gamiz, A. Godoy, F. Garcia-Ruiz, V. Georgiev, S. Amoroso, C. Riddet, E. Towie, A. Asenov, 3D multisubband ensemble Monte Carlo simulator of FinFETs and nanowire transistors, in: 2014 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 2014, pp. 21–24. doi:10.1109/SISPAD.2014.6931553. [11] L. Donetti, C. Sampedro, F. G´ amiz, A. Godoy, F. J. Garc´ ıa-Ru´ ız, E. Towiez, V. P. Georgiev, S. M. Amoroso, C. Riddet, A. Asenov, Multi-Subband Ensemble Monte Carlo simulation of Si nanowire MOSFETs, in: 2015 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 2015, pp. 353–356. doi:10.1109/SISPAD.2015.7292332. [12] R. Venugopal, Z. Ren, S. Datta, M. S. Lundstrom, Simulating quantum transport in nanoscale transistors: Real versus mode-space approaches, Journal of Applied Physics 92 (7) (2002) 3730–3739. doi:10.1063/1.1503165. [13] S. Jin, M. V. Fischetti, T.-w. Tang, Modeling of electron mobility in gated silicon nanowires at room temperature: Surface roughness scattering, dielectric screening, and band nonparabolicity, Journal of Applied Physics 102 (8) (2007) 083715. doi:10.1063/1.2802586. [14] P. Palestri, N. Barin, D. Esseni, C. Fiegna, Stability of self-consistent Monte Carlo Simulations: effects of the grid size and of the coupling scheme, IEEE Transactions on Electron Devices 53 (6) (2006) 1433– 1442. doi:10.1109/TED.2006.874758. [15] A. Godoy, F. Ruiz, C. Sampedro, F. G´ amiz, U. Ravaioli, Calculation of the phonon-limited mobility in silicon Gate All-Around MOSFETs, Solid-State Electronics 51 (9) (2007) 1211–1215. doi:10.1016/j.sse.2007.07.025. [16] P. Lugli, D. K. Ferry, Degeneracy in the ensemble Monte Carlo method for high-field transport in semiconductors, IEEE Transactions on Electron Devices 32 (11) (1985) 2431–2437. doi:10.1109/T-ED.1985.22291. [17] A. Pacelli, U. Ravaioli, Analysis of variance-reduction schemes for ensemble Monte Carlo simulation of semiconductor devices, Solid-State Electronics 41 (4) (1997) 599–605. doi:10.1016/S0038-1101(96)001980. [18] A. Trellakis, A. T. Galick, A. Pacelli, U. Ravaioli, Iteration scheme for the solution of the two-dimensional Schr¨ odinger-Poisson equations in quantum structures, Journal of Applied Physics 81 (12) (1997) 7880–7884. doi:10.1063/1.365396. 7