A Multi-Method Simulation Toolbox to Study Performance and Variability of Nanowire FETs

An in-house-built three-dimensional multi-method semi-classical/classical toolbox has been developed to characterise the performance, scalability, and variability of state-of-the-art semiconductor devices. To demonstrate capabilities of the toolbox, a 10 nm gate length Si gate-all-around field-effect transistor is selected as a benchmark device. The device exhibits an off-current (IOFF) of 0.03μA/μm, and an on-current (ION) of 1770 μA/μm, with the ION/IOFF ratio 6.63×104, a value 27% larger than that of a 10.7 nm gate length Si FinFET. The device SS is 71 mV/dec, no far from the ideal limit of 60 mV/dec. The threshold voltage standard deviation due to statistical combination of four sources of variability (line- and gate-edge roughness, metal grain granularity, and random dopants) is 55.5 mV, a value noticeably larger than that of the equivalent FinFET (30 mV). Finally, using a fluctuation sensitivity map, we establish which regions of the device are the most sensitive to the line-edge roughness and the metal grain granularity variability effects. The on-current of the device is strongly affected by any line-edge roughness taking place near the source-gate junction or by metal grains localised between the middle of the gate and the proximity of the gate-source junction.


Introduction
Gate-all-around nanowire field-effect transistors (GAA-NW FETs) are one of the main contenders for future CMOS technologies [1] since they provide a better electrostatic control of the channel when compared to fin field-effect transistors (FinFETs) [2], the current architecture adopted by the semiconductor industry. In addition, nanowire based transistor architectures extend beneficial properties of multi-gate devices required in digital circuits such as quasi-1D current transport, largely confined electrical fields, and immunity of threshold voltage from substrate bias [3].
On the other hand, GAA-NW FETs, like all deeply scaled semiconductor devices, are greatly affected by variability issues [4], related to either the fabrication process or material variations that can limit their performance and reliability [5]. Previous studies have shown that in the sub-threshold region, GAA-NW FETs are less resilient to intrinsic sources of variability than FINFETs [6,7]. In addition, the significant degradation observed in the on-region performance of GAA-NW FETs due to line edge roughness variations could be a critical issue for the scaling of these devices [6].
Nowadays, technology computer-aided design (TCAD) tools play a key role in the advancement of the semiconductor industry [8]. The TCAD tools are able to quickly characterise semiconductor devices, not only fabricated but also foreseen, and allow to investigate the impact of changes in materials, designs or fabrication processes. Currently, three-dimensional (3-D) simulations are necessary to appropriately model devices such as FinFETs or GAA-NW, due to the two-dimensional (2-D) nature of the quantum confinement, which increases the computational cost of a study [9]. There are different approaches that can be used in simulations of state-of-the-art semiconductor devices, ranging from the relatively simple and low computationally demanding drift-diffusion method [10], to extremely complex quantum mechanical approaches, such as pseudopotential-based electron quantum transport [11] or the non-equilibrium Green's functions (NEGF) [12] formalism, that can also be coupled to empirical tight-binding models [13]. The use of fully quantum simulators is computationally prohibitive for statistical studies, being essential a trade-off between the simulation's accuracy and the calculation time.
In this work, an in-house built finite-element multi-method semi-classical/classical simulation toolbox acronymed VENDES (Variability Enabled Nanometric DEvice Simulator) is used to characterise nano-scaled semiconductor devices including their operational performance and variability. To demonstrate capabilities of the VENDES, the performance and variability of a 10 nm gate length Si GAA-NW FET scaled down from an experimental device [14] is studied and assessed. The paper is organized as follows: the simulation techniques available in VENDES used in this study are described in Section 2. Section 3 analyses the performance and resilience to variability of the 10 nm gate length Si GAA-NW FET. Finally, Section 4 draws the main conclusions of this work.

Simulation Framework
VENDES, a 3D finite-element (FE) based device simulator, has been developed jointly at Universidade de Santiago de Compostela (Spain) and at Swansea University (United Kingdom). Figure 1 shows the basic flowchart of the VENDES toolbox.
The starting point is the generation of the FE mesh via the open source software Gmsh [15]. The FE method allows not only an accurate description of complex simulation domains, as in the case of elliptic cross-section shaped GAA-NW FETs [16], but also the possibility of introducing realistic deformations to the device dimensions. This capability of accurate geometrical description is crucial in the modelling of variability effects because a correct distribution of potential and carrier density is essential to predict the experimentally observed behaviour. Note that, for devices deeply scaled into the nanometre regime, the size of device variations and deformations can be comparable to critical device dimensions. The sources of variability included in VENDES either alter the structure dimensions or modify some physical properties affecting the device nodes and are described in detail in Section 3.2.
The classical electrostatic potential, V cl , is obtained from the Poisson equation solution on every node of the 3D FE tetrahedral mesh: where r = (x, y, z) is the spatial coordinate, ε(r) is the dielectric constant of the material, n(r) and p(r) are the electron and hole densities and N + D (r) and N − A (r) are the effective doping concentrations of donors and acceptors, respectively. Quantum corrections are incorporated in VENDES via two different techniques: i) the 3D density gradient (DG) equation and ii) the 2D Schrödinger (SCH) equation. In the first case, the DG quantum potential for electrons, V dg (r) [17], is obtained as: where Here, φ n (r) is the quasi-Fermi potential for electrons, n i (r) is the intrinsic carrier concentration of electron and holes, k B is the Boltzmann constant, T is the lattice temperature,h is the reduced Planck constant, r n is a dimensionless parameter that models statistical phenomena [18] and m x , m y and m z are the DG electron effective masses in the x-, yand z-directions, respectively [17]. It is important to remark that these effective masses operate as fitting parameters [19] and are not related to the material transport effective masses. The DG effective masses in the transverse directions (m y and m z ) will account for the strength of the quantum-mechanical confinement of the carriers in the device channel through a threshold voltage shift [20]. The DG effective mass in the transport direction (m x ) can account for the source-to-drain tunnelling by lowering the barrier of classical electrostatic potential which occurs between the source and the drain when the transistor is operating in the sub-threshold region [20]. The main drawback of the DG based quantum corrections is that they require calibration against either experimental data (when available) or more complex simulation techniques (such as Monte Carlo or Non-Equilibrium Green's Functions) [21]. On the other hand, the quantum correction method based on the solution of the Schrödinger equation [22] is calibration free. This technique assumes longitudinal and transverse electron effective masses in a minimum of the conduction valley of silicon and accounts for wave-functions penetrating into a surrounding dielectric layer [16,23]. The SCH equation is solved on two-dimensional (2D) slices placed across the device channel using a non-uniform distribution dependent on the gradient of electron density. The 2D quantum-mechanical electron density, n sc (y, z), is obtained from the SCH equation eigen-states, ψ i (y, z; E i ), and their corresponding eigen-energies, E i , as follows: where E F n is the electron quasi-Fermi level, g the degeneracy factor, and m * the electron effective transport mass. The Equation (4) considers Boltzmann statistics and assumes six equivalent valleys for Si (g = 6). The SCH quantum correction can be considered 'isotropic', when the electron effective mass in silicon is taken to be average of longitudinal and transverse electron effective masses, but also 'anisotropic', when longitudinal and transverse electron effective masses that are dependent on the valley orientation are considered. In that case, Equation (4) is solved separately for each of the three ∆ valleys, taking into account the different sub-band edges (i.e., appropriate energy levels) for the different valleys, obtaining a different n sc (y, z) for each valley (g = 2) and m * will be dependent on the channel orientation which can be 100 or 110 as shown in [24]. The electron density, n sc (y, z), calculated on the 2D slices, is interpolated to a 3D device density domain to obtain n sc (r). The resulting SCH quantum correction potential, V sc (r) [22], is as follows: Note that, in the anisotropic SCH quantum correction, a separate V sc (r) is obtained for each valley.
To simulate the transport inside the channel of the device, VENDES has implemented two different carrier transport methods: i) the drift-diffusion (DD) approach and ii) an ensemble Monte Carlo (MC) technique. The DD approach couples the electrostatic potential obtained from the quantum corrected solution of Poisson equation with the current continuity equation for electrons in order to obtain the electron current density, J n (r), as: div(J n (r)) = qR(r), where µ n (r) is the electron mobility and R(r) is the recombination term (set to zero by default). Note that, the DD method only accounts for the local relationship between the velocity and the electric field and it is unable to correctly represent non-equilibrium transport effects [25]. However, some of the non-equilibrium phenomena can be partially mimicked via appropriate mobility models. To model the carrier transport behaviour in GAA-NW FETs, VENDES uses Caughey-Thomas doping dependent low-field electron mobility model [26] coupled with perpendicular and lateral electric field models [27] which better describe carrier transport at large electric fields. When using these mobility models, the main calibration parameters are a low-field carrier mobility, a critical electric field, E cn , and a saturation velocity, v sat . The limitations of the DD approach can be overcome by using a semi-classical transport model, the MC technique where an ensemble of particles representing carriers evolves through free flights governed by Newton equations and undergoes scattering events with a probability which is determined quantum-mechanically. The best way to initialise the distribution of carriers in the real space, is to use the quantum corrected potential from the solution of the Poisson's equation, which results in speeding up the simulation. MC uses the analytical non-parabolic anisotropic approximation [28] for the silicon band structure taking into account three valley minima, X, L, and Γ, using Herring-Vogt transformation to transform ellipsoidal surfaces to spherical ones in order to simplify a calculation of free flights and scattering events. The MC technique considers carrier scattering in a quantum-mechanical way by using typically Fermi Golden Rule [29] to obtain the transition rates. The following electron scattering mechanisms, important for silicon based devices, are included in VENDES: i) electron interactions with intra-and inter-valley acoustic and non-polar optical phonons [28,30], ii) electron interactions with ionised impurities using Ridley's third body exclusion [31,32] and static screening [29], and iii) electron interaction with interface roughness using Ando's 2D potential approach [33]. VENDES uses Boltzmann statistics when solving 3D Poisson equation and determining a final state after electron scattering but, the electron scattering with ionised impurities uses Fermi-Dirac statistics to calculate the static screening by a self-consistent calculation of the Fermi energy and the electron temperature from the average electron density and kinetic energy in a whole real space device simulation domain at each scattering event [34]. The inclusion of Fermi-Dirac statistics into electron scattering with ionised impurities turns to be sufficient to correctly simulate injection of carriers into the channel from a heavily doped source/drain when comparing the results from quantum corrected 3D finite element Monte Carlo device simulations with experimentally measured I-V characteristics in nanoscale FinFETs [29] and nanowire FETs [16].

Performance and Variability of GAA-NW FETs
In this work, VENDES has been applied to study state-of-the-art nanoscale GAA-NW FETs designed for future digital technology node generations [35]. Section 3.1 presents the GAA-NW FET description and main figures of merit. Section 3.2 shows a thorough analysis of the impact that different sources of fluctuations have on this architecture.

Benchmark Device
The device under study is a 10 nm gate length Si GAA-NW FET with an elliptically shaped cross-section that has been scaled [16] following the ITRS guidelines [35] from an experimental 22 nm gate length device from IBM [14]. Ref. [14] includes TEM images of the fabricated structure and the I D -V G characteristics that we have used to validate our device. The elliptically shaped cross-section of the transistor body as well as a lateral shape is a result of the advanced fabrication process in which the shape formation is mostly affected by etching. Figure 2 shows a comparison of experimental I D -V G characteristics of the 22 nm gate length GAA-NW FET versus the simulation results provided by VENDES SCH-MC. The drain bias is 1.0 V. The empirical doping values were not included in [14], so they were reversed engineered following the methodology described in [16]. Note that, the MC device simulations in VENDES are able to accurately reproduce the experimental results in all the active regions of the device, except for at a very low gate bias of 0.0 V, where the MC statistical noise is too high. At the very low gate bias, the DD device simulations are typically used.  Table 1. The gate work-function (WF) was set to 4.4 eV. For this device, Figure 3 shows the I D -V G characteristics at a high drain bias of 0.7 V in both linear and logarithmic scales for DG-DD, SCH-DD and SCH-MC simulations. Note that SCH-MC simulations are calibration free, whereas the DG and DD models need to be properly fitted (see the main calibration parameters in Table 1) in order to achieve the agreement shown in Figure 3. The main figures of merit (FoM) that characterise the I D -V G characteristics are shown in Table 1.
FoMPy [36,37] is a python-based open source post-processing tool implemented in VENDES (see Figure 1) that automatically extracts the main FoMs of a I-V characteristics. This tool is very useful when performing statistical studies, where a large ensemble of devices needs to be analysed. In this work, the threshold voltage (V T ) has been obtained using the second derivation method, the off-current is obtained at a 0.0 V gate voltage, and the on-current has been extracted at a gate bias equal to V T + V DD , being V DD the supply voltage (set to 0.7 V). The analysed device has a low off-current of 0.03 µA/µm, acceptable for applications in mobile low power devices with a long battery life, and an on-current of 1770 µA/µm, that has been achieved by increasing the maximum S/D doping from 5 × 10 19 cm −3 , used in the 22 nm gate length experimental device, to 10 20 cm −3 . This increase in the doping has allowed to raise the device on-current by 40% (as previously shown in [38]), at the cost of a slight deterioration in the sub-threshold slope (SS). The device SS is 71 mV/dec, not far from the ideal limit of 60 mV/dec. Therefore, the device doping is one of the key design parameters that needs to be considered when designing a device for a specific application. For transistors aimed at high performance (HP), standard performance (SP) and low power (LP) applications, the I ON /I OFF ratio is also a key parameter because it provides a global characterisation of the device operation. The observed I ON /I OFF (6.63 × 10 4 ) is 27% larger than that of a similar gate length Si FinFET [7] of 10.7 nm.
The MGG is modelled by altering the work-function of the device gate so it matches metal grain distributions either observed empirically via KPFM [42], or generated using the Voronoi approach, where the experimental shapes and values of different grain orientations are mimicked [43]. Figure 4a shows an example of a Voronoi TiN metal profile applied to the device gate. The WF values are 4.4 eV and 4.6 eV and their respective probabilities of occurrence 40% and 60%. The average grain size is 5 nm.
RD are introduced in the n-type doped S/D regions using a rejection technique from the doping profile (shown in Table 1) of the ideal device. Initially, dopants with their associated charge are distributed on an atomistic grid defined by the location of the atoms. Then, this charge is mapped to the device tetrahedral mesh using the cloud-in-cell technique, in order to generate an atomistic electron density distribution [44], as shown in Figure 4b.
GER and LER are modelled similarly, the device gate (in case of the GER) or the edge of the nanowire (in case of the LER) are deformed according to the shape of a given roughness profile created via the Fourier synthesis method [45]. Two parameters are used to characterise these deformations: i) the root mean square (RMS) height, that sets the amplitude of the roughness, and ii) the correlation length (CL), that accounts for the spatial correlation between the deformations in the different points of the device.  For each variability source, ensembles of 300 device configurations were created and simulated at a high drain bias of 0.7 V. Figure 5 shows the impact of the aforementioned sources on the 10 nm GAA-NW FET threshold voltage variability. The statistical sum of the four sources of variability (COMB) has also been included as comparison. Results show that GAA-NW FETs are heavily influenced by the LER variability in the sub-threshold region, with σV T values 1.4 and 2.0 times larger than those of the MGG and the RD, respectively. The GER is the least influential source of variability having its σV T a 86% lower than that of the LER. Note that the combination of the four sources of variability leads to a V T standard deviation of 55.5 mV, a value 85% larger than the one observed (σV T = 30 mV) in a similar gate length Si FinFET [7]   Variability studies are highly computational demanding because they require the simulation of hundreds or thousands of device configurations in order to obtain results with statistical significance. Table 2 shows, for the different simulation methodologies implemented in VENDES, the total times for the solution of one I D -V G bias point at a high drain bias of 0.7 V on a single core for two different CPUs. Note that the simulation time for a SCH-MC simulation is around 70 times longer than a quantum-corrected DD study. For that reason, VENDES performs sub-threshold region variability studies (see the flowchart in Figure 1) using either DG or SCH-DD simulations. The reason for this is twofold: i) in the sub-threshold, the electrostatics dominate and quantum-corrected DD simulations are able to provide accurate results and, ii) MC results can be extremely noisy at very low gate biases and lead to incorrect off-current or sub-threshold slope values. However, in the on-region regime, VENDES performs the variability studies via the SCH-MC simulations since the DD approach is unable to capture a non-equilibrium carrier transport even if it is properly calibrated, leading to large over-or under-estimation of the variability [38]. However, it is important to remark (as seen in Figure 3), that both SCH-DD and MC-DD simulations match perfectly at the threshold.
In the quest for the reduction of the computational time, several alternatives have been investigated: i) the parallelisation of the simulation code using a message passing interface (MPI), explained in detail in [46], to take advantage of increasingly available computational infrastructures, such as clusters and supercomputers and, ii) the development of a methodology to predict the impact of the variability sources. A sequential simulation of one I D -V G bias point using the SCH-DD method is 960 s, see Table 2 for runs using the Intel(R) Xeon(R) CPU. When using a parallel version of the code with 2 and 4 processors, this time is reduced to 613 s (78 % parallel efficiency) and 363 s (66 % parallel efficiency), respectively. On the other hand, the Fluctuation Sensitivity Map (FSM) approach [47,48] is a methodology that we developed to predict the impact of the variability sources. This post-processing tool (see the VENDES flowchart in Figure 1) is based on the creation of a map that provides information of the sensitivity of the different regions of a device to a particular source of variability. Once the map is created, it can be used to predict the statistical variability, with a very small error, under different input parameters as shown in Ref. [47]. In a typical variability study, we simulate at least 300 device configurations per variability source and characteristic parameter. This characteristic parameter can be the grain size in the MGG study, or the CL and RMS in a study of LER and GER. For instance, a full MGG variability study will require simulations of several grain sizes (a minimum of three). Therefore, the total computational cost of this full study using the SCH-DD or SCH-MC methods will be 300 h and 31500 h, respectively (for the Intel(R) Xeon(R) CPU in the sequential case). These times can be reduced by 66% using the FSM because once one set of 300 device configurations is simulated and used to create the FSM (which can take up to 2 min), the map can be used to predict the variability results for the remaining grain sizes without any further statistical computations.
The FSM can be also used as an assistance in the design of variability-resistant device architectures since it pinpoints to parts of the device the manufacturers should concentrate their efforts on. Figure 6 shows an example of the on-current FSM obtained when a single LER deformation is applied to a specific location of the device, narrowing its width. Using this synthetic deformation, it is possible to sweep all the locations along the device, measuring the changes in the FoM that enable us to the spatial sensitivity to LER variations. Note that, in Figure 6 (bottom), a negative (positive) sensitivity indicates an increase (decrease) in the on-current. The effect of the synthetic deformation on the on-current depends on its position along the transistor. Any change in the NW width happening near the source-gate junction will heavily impact the device on-current whereas if the deformation occurs near the source or drain ends, its impact on the on-current will be minimal. Similarly, the FSM can be applied to other sources of variability, like the MGG, taking into account that the generated fluctuation map will now be two-dimensional (2D), in order to characterise the whole device gate. Figure 7a shows a scheme of the device that has a fixed WF value of 4.6 eV in all the gate except from a narrow strip in which the WF is 4.4 eV. This narrow strip is swept along the gate and for each position, the device configuration was simulated at a high drain bias of 0.7 V and the corresponding on-current was extracted. The resulting 2D on-current FSM due to the MGG is shown in Figure 7b. This map allows us to establish that any variation in the WF located between X = 0.0 nm (middle of the gate) and X = −2.5 nm, will have the largest impact on the device performance. However, the GAA-NW FET will be practically insensitive to WF variations happening in the proximity of the source (X = −5.0 nm) or the drain (X = 5.0 nm) junctions.

Conclusions
VENDES, an in-house-built 3D multi-method device simulator, has been used to characterise the performance and resistance to variability of a 10 nm gate length Si GAA-NW FET scaled down from an experimental transistor [14] following ITRS guidelines [35]. The off-current of the device is 0.03 µA/µm, and the on-current of 1770 µA/µm, delivering the I ON /I OFF ratio of 6.63 × 10 4 . The device SS is 71 mV/dec, not far from the ideal limit of 60 mV/dec. σV T due to the statistical combination of LER, GER, MGG and RD is 55.5 mV, a value significantly larger than that of a similar gate length Si FinFET of 10.7 nm (30 mV). This larger threshold voltage variability indicates that the variability effects may be another limiting factor for the adoption of the GAA-NW FETs in the future technological nodes.
Finally, the FSM allowed us to determine which regions of the device are the most sensitive to the LER and MGG variations and influence the device characteristics the most. In the case of the LER, the changes in the device width occurring near the source-gate junction will heavily impact the device on-current whereas any deformation happening near the source or the drain, will have a negligible influence on the on-current. In the case of the MGG, the most sensitive region of the device is localised between the middle of the gate and the proximity of the gate-source junction. The information provided by the FSM can be very useful as an aid for the creation of fluctuation-resistant device architectures. Funding: This research was supported in part by the Spanish Government under the projects TIN2013-41129-P, TIN2016-76373-P and RYC-2017-23312, by Xunta de Galicia and FEDER funds (GRC 2014/008) and by the Consellería de Cultura, Educación e Ordenación Universitaria (accreditation 2016-2019, ED431G/08).

Acknowledgments:
The authors thank Centro de Supercomputación de Galicia (CESGA) for the computer resources provided.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: