Multilevel 3-D Device Simulation Approach Applied to Deeply Scaled Nanowire Field Effect Transistors

Three silicon nanowire (SiNW) field effect transistors (FETs) with 15-, 12.5- and 10.6-nm gate lengths are simulated using hierarchical multilevel quantum and semiclassical models verified against experimental <inline-formula> <tex-math notation="LaTeX">${I}_{D}$ </tex-math></inline-formula>–<inline-formula> <tex-math notation="LaTeX">${V}_{G}$ </tex-math></inline-formula> characteristics. The tight-binding (TB) formalism is employed to obtain the band structure in <inline-formula> <tex-math notation="LaTeX">$\mathit {k}$ </tex-math></inline-formula>-space of ellipsoidal NWs to extract electron effective masses. The masses are transferred into quantum-corrected 3-D finite element (FE) drift-diffusion (DD) and ensemble Monte Carlo (MC) simulations, which accurately capture the quantum-mechanical confinement of the ellipsoidal NW cross sections. We demonstrate that the accurate parameterization of the bandstructure and the quantum-mechanical confinement has a profound impact on the computed <inline-formula> <tex-math notation="LaTeX">${I}_{D}$ </tex-math></inline-formula>–<inline-formula> <tex-math notation="LaTeX">${V}_{G}$ </tex-math></inline-formula> characteristics of nanoscaled devices. Finally, we devise a step-by-step technology computer-aided design (TCAD) methodology of simple parameterization for efficient DD device simulations.

T HE 3-D multigate architectures are ruling industry lines since fin-like field-effect transistors (FinFETs) were adopted as the new standard architecture. The faster and less power-hungry FinFETs facilitate continuous scaling of semiconductor transistors [1] by increasing transistor density on a chip and reducing costs [2]. However, this architecture cannot provide the required electrostatic integrity when scaled to future sub-5-nm technology nodes [3]. Hence, gate-allaround nanowire (NW) FETs are considered as a replacement [4] thanks to their superior electrostatic integrity against FinFETs [5] and nanosheet FETs [6].
Modeling of semiconductor devices via technology computer-aided design (TCAD) undoubtedly plays an essential role in the scaling of transistors [7]. TCAD tools based on quantum-mechanical corrected classical and semiclassical transport methods can provide a good tradeoff between the simulation precision and computational cost. The quantum corrections can be included via the density gradient (DG) approach, limited by the fitting of DG effective masses [8], or through the solution of the Schrödinger (SCH) equation, which is calibration-free [9]. However, when studying devices with channel cross section dimensions below 10 nm, geometry starts to play a substantial role in the band-structure [10]. The quest for an accurate band-structure leads to a couple of approaches to use. The effective mass approximation [11] is relatively simple and accurate near the conduction band (CB) minima [12] but will struggle to correctly estimate the quantization levels in nanostructures [10]. On the other hand, the sp 3 d 5 s * nearest-neighbor empirical tight-binding (TB) model [13] gives a full Brillouin zone description of the bandstructure, allowing to atomistically discretize a realistic device at nanometer scale [14]. The effects of the quantum confinement band-structure in silicon NWs (SiNWs) were studied before [15]- [17], however, previous works: 1) calculated only a low-field mobility and/or 2) did not compare the results against experimental data obtained from a SiNW FET.
In this work, we use multilevel 3-D material and device simulations to obtain the I D -V G characteristics for three ellipsoidal 110 SiNW FETs with 15-, 12.5-, and 10.6-nm gate lengths, using bulk and quantum-mechanically confined electron effective masses of silicon. These gate lengths were chosen following the International Roadmap for Devices and Systems (IRDS) predictions for sub-5-nm industry nodes [18]. The 15-nm gate length simulated results are initially compared to experimental data [19]. Then, we analyze the role of quantum-mechanical confinement on the silicon bandstructure for the three SiNW FETs. Finally, we present a step-by-step methodology of a relatively simple parameterization for speedy drift-diffusion (DD) quantum-corrected simulations.

II. SIMULATION METHODOLOGY
To model the SiNW FETs, we use a multiscale simulation approach. Initially, we perform empirical TB band-structure calculations for the SiNWs to provide bandstructure parameters for electron transport simulations. These parameters are introduced into an in-house-built V ariability Enabled N anometric Device Simulator (VENDES) [20] that include 3-D finite-element (FE) self-consistent classical and semiclassical transport models. The classical transport approach is the DD method, which is a workhorse of many commercial TCAD tools [21]. The semiclassical transport model is the ensemble Monte Carlo (MC) technique, which is widely recognized as an accurate predictive simulation technique for semiconductor devices operating in a highly nonequilibrium transport regime [22]. Note that transport simulations are selfconsistently run with solutions of the Poisson equation, using the linear scheme, accounting for long-range electron-electron interactions [7].
Since quantum-mechanical confinement will play a role in the carrier transport in the SiNW FETs, both transport approaches incorporate quantum confinement corrections via the solution of 2-D FE SCH equation. The quantum corrected DD model of carrier transport is used when the diffusion transport dominates, in the subthreshold/OFF-region of device operation. The ensemble MC transport technique is used when a highly nonequilibrium carrier transport dominates, in the ON-region of device operation. This quantumcorrected ensemble MC technique was proven to reproduce the highly nonequilibrium carrier transport in weakly quantummechanically confined systems in [23] and [24]. The advantages of the quantum corrected 3-D FE MC simulations are: 1) the solutions of 2-D FE SCH equation accurately determine the carrier density distribution along the channel [25]; 2) the device domain described by the 3-D FE mesh precisely reproduces nonequilibrium electron transport in the device active region including interface roughness (IR); and 3) the accurate injection of electrons into the channel from a heavily doped source/drain (S/D) thanks to self-consistent Fermi-Dirac statistics used for the electron scattering with ionized impurities.

A. Tight-Binding Model
A 20-band sp 3 d 5 s * nearest-neighbor TB model [13] is selected in order to obtain realistic dispersion relations of the SiNWs. The parameterization by Niquet et al. [26] was chosen because of its accurate reproduction of both the experimental electron and hole effective masses along all directions in kspace, crucial for accurate determination of the confinement energies. This parameterization yields a fundamental bandgap energy of 1.17 eV for bulk Si. Thus, all the CB levels obtained from the TB calculation were downshifted by 0.05 eV to account for finite temperature, and the gap values we report hereafter reflect this correction. In the calculation, surface dangling bonds were passivated by H atoms, with the Si-H bonds parameterized according to the values provided in [27]. Fig. 1 shows the bands for two SiNWs with ellipsoidal cross  5.06 × 6.36 nm 2 cross section SiNWs. The energy levels have been downshifted 0.05 eV to set the bulk CB edge at 1.12 eV. Blue (light gray) circles represent Si (H) atoms. The levels originating from the bulk Δ Z pocket and the folding of the bulk Δ X , Δ Y pockets are labeled. Also, Δ 1 , Δ 2 and Δ 3 indicate the bands that will be treated explicitly by the MC simulator. In particular, Δ 1 for the large SiNW captures the behavior of many Γ subbands.
sections, indicating the minima originating from the Z pocket and the folding of the X , Y pockets.

B. Ensemble Monte Carlo Simulations
The MC engine simulating kinetics of electrons in the semiconductor device utilizes a nonparabolic anisotropic model within the effective mass approximation [11], [28] and it includes electron scattering mechanisms assuming first-order perturbation theory (Fermi Golden Rule). It considers electron scattering with acoustic and nonpolar optical phonons (intravalley and intervalley) [22], [29], ionized impurity scattering using the third body exclusion [30], and IR scattering using Ando's model [31]. The electron screening in ionized impurity scattering incorporates a self-consistent calculation of Fermi energy and electron temperature from the electron density and the average electron energy [32] at every real position of the electron in a 3-D simulation domain. More details on the 3-D FE MC engine can be found in [9] and [24].

C. Schrödinger Equation Quantum Corrections
The quantum corrections are incorporated by solving the 2-D FE SCH equation repeatedly across 2-D slices perpendicular to the transport direction, equidistantly placed along the device. The SCH equation is solved in a minimum of the conduction valley of Si assuming longitudinal and transverse electron effective masses while still considering the penetration of wave-functions into a surrounding dielectric layer [9], [23]. The SCH equation quantum corrections assume Boltzmann statistics but are anisotropic since the longitudinal and the transverse electron effective masses depend on the valley orientation with a degeneracy factor of 2. Quantum electron density, n Sch (y, z), is thus obtained separately for each of the three valleys in Si [33]. For each valley, this 2-D electron density is interpolated into a 3-D FE domain to calculate n Sch (r). This electron density n Sch (r) is then used to calculate the quantum-corrected electrostatic potential V QC (r) [34].

D. Parameter Extraction for Drift-Diffusion
For fast turnaround in performance screening, the DD transport model is often desirable, despite assuming a nearequilibrium carrier transport in nanoscale transistors for all biases. We generate a workflow to turn the dispersion relations from the TB model in Section II-A into parameters suited for DD simulations as follows.
1) Conduction subband minima, E ci , and effective masses, m * i , are computed by diagonalization of the TB Hamiltonian and extraction from the subband curvatures.
2) The density-of-states (DOS) effective mass in 1-D is obtained by a proper weighting of the individual subband DOS effective masses as where k B is Boltzmann's constant, T is the temperature, i runs over the CB subbands, g i is the degeneracy of a subband excluding spin (g i = 1 for a spin-degenerate 1-D subband), and E ci − E c0 is the energy distance between the i th subband and the lowest CB energy.
where A is the NW cross section. 4) The average scattering time ρ and the transport effective mass m * DD,t,1−D are calculated from an experimental or computed diffusive mobility, μ diff , which is the long channel mobility, obtained when L λ, being L the effective channel length and λ the average mean free path [35], by solving the equation 1/m * which allows for a regime between diffusive and ballistic transport. 5) The unidirectional thermal velocity [35] is calculated as 6) λ is then obtained from the relation 7) The apparent mobility, μ app , combines the diffusive and ballistic mobilities [35] as 8) The critical length [36], L crit , in the saturation regime is given by where the prefactor is obtained by fitting the data in [36, Fig. 8(a)]. This prefactor shows a very little dependence on the SiNW diameter or the type of carrier. 9) The transmission coefficient in the saturation regime is 10) The saturation velocity is finally computed as [36] v

A. Benchmark Devices
The theoretical models and numerical methods included in VENDES are initially validated against experimental I D -V G characteristics for a 15-nm gate length SiNW FET that has an elliptical cross section of 7.15 × 8.99 nm 2 [19]. A scheme of the device structure is illustrated in Fig. 2, with the dimensions and doping profiles summarized in Table I. The 15-nm gate FET doping profile is extracted by a reverseengineering process by simulating the device characteristics in the subthreshold [23]. The SiNW FETs are assumed to have a uniformly p-type doped channel and a Gaussian n-type doping in the S/D regions. The doping profile including a maximum doping and a lateral straggle of the Gaussian S/D doping is adjusted until the subthreshold slope of the simulated I D -V G  [18] following the same methodology as in [23].

B. DD and MC Parameters for the Benchmark Devices
The application of the procedure in Section II-D to extract the relevant DD simulation parameters for the benchmark devices produces the DD DOS electron effective mass, electron saturation velocity, apparent electron mobility, perpendicular critical electric field, and the critical length in saturation collected in Table II. VENDES uses the Caughey-Thomas doping and temperature-dependent low-field carrier mobility model [37] that defines the low-field electron mobility, μ low , as where N D and N A are the donor and the acceptor doping concentrations. μ min and μ max indicate the mobility for high or low doped materials, respectively. In our simulations, for a temperature of 300 K, we take μ diff from [36], and we consider that μ min = μ app (see values in Table II) [38]. To describe the carrier transport between low-field and high-field transport conditions, the electron mobility dependence on the electric field (E) [37], [39] is included as being E and E ⊥ the parallel and perpendicular electric fields, respectively, β a fitting parameter equal to 2 and, E c the perpendicular critical electric field (see Table II). The results from the TB calculations inform the MC simulator and the SCH equation solver of energy gaps and electron effective masses (see Table III). Since the VENDES explicitly treats three bands, the following assumptions are made to incorporate many TB bands into VENDES. The 15-and 12.5-nm gate length NW FETs have many bands originating from the point which fall inside the energy window relevant for the electron transport [see Fig. 1(a)]. Since a large density of subbands indicates energy spectrum close to the bulk, the longitudinal effective masses in the transport direction are set to the bulk values ( 1 in Table III), and the effects of level quantization are accounted for by using a 3-D DOS effective mass with a value reduced with respect to bulk, calculated from the individual band curvatures and including the occupation probability, given by (1) and (2). On the other hand, the two lowest bands with minima close to the Brillouin zone edge Z are treated explicitly ( 2 and 3 ), with the 3-D DOS effective mass considering a possible near-degeneracy. To prevent the electron movement in the transverse-to-transport direction, both transverse (m t ) effective masses are set to a value of 1.9m 0 for all three valleys. This large transverse mass will: 1) prevent the electron motion along the transverse directions due to quantum confinement and 2) assure that the occupation in phase space does not artificially alter the scattering rate. Note that, the chosen 1.9m 0 reproduces the experimental drain current for the 15-nm gate length device. Indeed, the drain current is sensitive to m t , observing a 22/38% reduction/increase in the ON-current when m t is increased/reduced by 26%.
The 10.6-nm gate length device is qualitatively different. The quantization of the energy levels at the point is so strong that only a few bands are within the energy window available to electron transport [see Fig. 1(b)]. In this case, we treat the two bottom bands at ( 1 and 2 ) and the bottom band close to Z explicitly ( 3 ) using longitudinal effective masses obtained from the band curvatures, and the DOS masses reflecting a possible near-degeneracy. The transverse effective masses for the 10.6-nm gate length device are also set to a value of 1.9m 0 as in the larger cross section devices.

IV. SINW FETS I-V CHARACTERISTICS
In this section, we present the I D -V G characteristics obtained using the quantum-corrected MC and DD simulations, compare them against experimental data, and analyze the impact of the TB extracted parameters on the results.  Table III are plotted in Fig. 3 for the 15.0-nm gate length SiNW FET. There is an excellent agreement between the experiment and the simulations at both low and high drain biases. The external series resistance has been accounted for by increasing the length of the S/D regions up to 80 nm (see Table I). Fig. 4 shows, on both linear and logarithmic scales, a comparison of the experimental I D -V G characteristics against the 3-D   SCH equation quantum-corrected DD simulations (DD-SCH) (see blue circles) that incorporate the TB-generated parameters presented in Table II. The DD simulations are able to precisely capture the experimentally observed subthreshold region (note the good match in the subthreshold slope and in the vicinity of the threshold voltage) thanks to the inclusion of SCH equation quantum corrections. However, in the ON-region, the DD simulations predict a 1.19 times larger ON-current than the experimentally seen. This failure to properly capture the nonequilibrium effects happening in the ON-region of the device is a well-known limitation of the DD transport model [40]. However, this difference can be used as a fitting parameter to fully calibrate the I -V curve, as done in Fig. 4 (pink hexagons), if the SCH equation quantumcorrected DD simulations were to be used in an ON-region study.   Fig. 5 presents, for the three studied devices, a comparison of the MC-SCH simulated I D -(V G − V T ) characteristics when using either the TB-generated or the bulk parameters. For the bulk Si bandstructure parameters, in a 110 channel crystallographic orientation, we consider a longitudinal electron effective mass of 0.916m 0 , transverse electron effective masses of 0.190m 0 , a 3-D DOS electron effective mass of 0.553m 0 , and a band-gap energy of 1.12 eV. As expected, the use of bulk parameters provides an accurate characterization of the I -V curve in the subthreshold region and only a minor overestimation (around 4%) of the current in the ON-region for the device with a 15.0-nm gate length. The SiNW FET scaled to a 12.5-nm gate length has very similar I -V characteristics in the subthreshold region when using TB or bulk parameters.

B. Bulk Versus TB Parameters
The current starts to differ only from (V G − V T ) > 0.4 V with a difference up to 8%. When scaled to the 10.6-nm gate length, the differences in current appears at lower gate biases of (V G − V T ) > 0.25 V with a 87% overestimation in the ON-current when the bulk parameters are used instead of the TB generated ones.

C. DD-SCH for Strongly Confined SiNW FETs
The MC simulation results serve as the reference to validate the DD-SCH simulations, since no experimental data are available for the smaller SiNW FETs. The MC simulations were previously verified against experimental data for different architectures: a 25-nm gate length SOI FinFET [24], a 22-nm gate length gate-all-around NW FET [23], and a 12-nm gate length nanosheet FET [6] with an exceptional agreement. Fig. 6 shows I D versus (V G − V T ) characteristics for the 12.5-and 10.6-nm gate length SiNW FETs comparing results from the DD-SCH and the MC-SCH simulations using the TB-extracted parameters. The DD-SCH simulations overestimate the ON-current when compared to the results from the MC-SCH. When using the TB-generated parameters, this overestimation is by a factor 1.19 (as previously seen in Fig. 4 for the 15-nm gate length device) indicating a systematic error. This overestimation can be mitigated by calibrating simulations of a large gate length device against available experimental data and extracting a correction factor that can be applied for smaller gate length dimensions with no experimental measurements. The introduction of the TB-extracted parameters into the quantum-corrected DD and MC simulations will thus allow us to perform accurate analyses of performance and variability for advanced architectures at an affordable computational cost.

V. CONCLUSION
We have presented a multilevel physics-based methodology to simulate semiconductor devices in the nanometer regime, where quantum-mechanical confinement effects play an essential role at equilibrium and in highly nonequilibrium carrier transport. This methodology couples TB calculations (to accurately capture the impact of device cross section on band-structure) with quantum-corrected DD and MC techniques. The bulk Si parameters and band-gap energy values are perfectly adequate to reproduce the experimental I D -V G characteristics for a 7.15 × 8.99 nm 2 cross section SiNW FET. However, for smaller cross section devices of 6.20 × 7.80 nm 2 and 5.06 × 6.36 nm 2 , the results obtained from the quantum-corrected MC using bulk parameters overestimate the ON-current by 8% and 87%, respectively, when compared to the results obtained from the MC-SCH simulations using the TB-extracted electron effective parameters. In addition, the quantum-corrected DD simulations can still provide accurate results by using a very small number of fitting parameters obtained from our multilevel physics-based methodology to extract them.