A New Approach to Modeling Ultrashort Channel Ballistic Nanowire GAA MOSFETs.

We propose a numerical compact model for describing the drain current in ballistic mode by using an expression to represent the transmission coefficients for all operating regions. This model is based on our previous study of an analytic compact model for the subthreshold region in which the DIBL and source-to-drain tunneling effects were both taken into account. This paper introduces an approach to establishing the smoothing function for expressing the critical parameters in the model's overall operating regions. The resulting compact model was tested in a TCAD NEGF simulation, demonstrating good consistency.

1. Introduction

As electronic device dimensions shrink, nanoscale devices with multigate structures are considered to improve immunity to short-channel effects (SCEs), such as drain-induced barrier lowering (DIBL) and threshold voltage roll-off [1,2,3,4,5,6,7,8,9]. In recent years, several novel structures of devices [10,11,12,13,14,15,16,17], such as nanowire (NW) and nanosheet (NSH) metal-oxide semiconductor field-effect transistors (MOSFETs), have attracted more attention and become potential structures that most likely led to downscaling into the sub-7 nm regime. Then, because of the reduction in scattering events in the short-channel nanowire GAA MOSFETs, carriers traveling through the channel without scattering increase ballistic transport [18,19,20,21,22]. Thus, the properties of such devices in which the quantum transport of carriers should be taken into account were intensively investigated and modeled [23,24,25,26,27,28,29].

Moreover, since the DIBL effect cannot be neglected in an ultrashort channel gate-all-around (GAA) MOSFET, it is more meaningful to develop a compact model incorporating ballistic transport, quantum confinement, and quantum tunneling effects. Then, the ballistic transport and quantum mechanism in ultrashort channel nanowire GAA MOSFETs were widely recognized as significant factors in determining drain current characteristics by using a numerical calculation method such as Schrödinger–Poisson and nonequilibrium Green’s function (NEGF) simulations [30,31,32,33]. The results of the current characteristics from the TCAD simulation and the experiment both demonstrate higher and better subthreshold swing (SS) characteristics in a GAA MOSFET than those in a Fin field-effect transistor (FinFET) in a wide range of temperature. They also indicate that GAA MOSFETs can provide more robust solutions to next-CMOS devices beyond FinFETs and low-supply-voltage applications [29,34,35]. However, the computations of the simulations could be time-consuming, even for estimating the characteristics of one single device. A compact model of nanowire GAA MOSFETs would be more efficient in describing the drain-current characteristics for a reduction in computational time when taking the ballistic transport and quantum effects into account.

Several compact models evaluating GAA-MOSFET properties were reported [10,36,37,38,39,40,41,42,43,44,45,46]. Some models focused on the expressions of the drain-current characteristic with a DIBL effect, incorporating the transport modes on the basis of the drift-diffusion [39,44] and ballistic theories [41]. Although the DIBL effect is mentioned in some of those models, the quantum confinement [39,45] and the source-to-drain tunneling effects [10,41,42] also need to be considered in the models. Consequently, the drain current would be represented by combining tunneling and thermionic ballistic transport expressions in the subthreshold and inversion regions [47]. However, a fully analytic compact model that can be introduced to SPICE-class simulators using Verilog-A language to implement a circuit-level simulation is expected to be mathematically developed to predict the behavior of a nanowire GAA MOSFET as a final purpose. According to our previous work, the potential distribution of a sub-band energy level and transmission coefficient of the source-to-drain tunneling can be expressed and determined approximately by using one parabolic function and the Wentzel–Kramers–Brillouin (WKB) approximation, respectively. Although source-to-drain tunneling was incorporated in the ballistic device with an ultrashort channel according to our previous study [47], it is increasingly necessary to simplify the calculation process in a compact model when evaluating the drain current by only using the Landauer formula in all operating regions.

In this paper, we propose one numerical compact model of the nanowire GAA MOSFET with an ultrashort channel to describe the ballistic drain current by introducing an approach of the smoothing function to determine the significant parameters through connecting the expressions of the subthreshold and inversion regions mentioned in our previous study [24]. Hence, we improved expressing the unknown parameter, the barrier top locations, and the transmission coefficients for all operating regions using the smoothing functions mentioned above. Then, the drain current for the ultrashort channel GAA MOSFET is evaluated from the Landauer formula, which is in line with NEGF simulation results. The rest of this paper is organized as follows. We introduce the model structure, definition of parameters, and the ballistic current model in Section 2. Section 3 illustrates how we use the smoothing functions to approximately express the unknown parameter, the barrier top location, and the transmission coefficients being valid for all operating regions. Lastly, the calculation results and discussion are in Section 4.

2. Device Structure and Model Formulas

2.1. Model Structure

Figure 1a,b show the schematic structures of the modeled nanowire GAA MOSFETs in the channel and transverse directions, respectively. In this structural model, we consider the source and drain electrodes as the ideal reservoirs of electrons since they were produced with highly doped n-type silicon that would provide sufficient electrons into the intrinsic silicon channel surrounded by gate oxide and metal gate with a circular cross-section. Moreover, the backscattering from the source and drain electrodes was neglected. Channel length, radius, and gate oxide thickness are denoted as , R, and , respectively. Correspondingly, the cylindrical coordinate with length, radius, and angular components along with z, r, and directions, as shown in Figure 1b, is used in the following formulations. However, on account of the perfect angular symmetry in our model, it was sufficient to represent the electrostatic potential distribution only along the channel and radius directions without considering the component.

Figure 1

Schematic structure (a) in the channel direction z and (b) in the transverse direction r of the nanowire GAA MOSFET model.

2.2. Potential Profile

The profiles of the conduction band edge and the higher sub-band energy levels along the z and r axes at and are shown in Figure 2a,b, respectively, where it was assumed that the short-channel effects to the channel along the transverse direction in this model were consistent. Then, as shown in Figure 2a, is the total energy of an electron, and the maximal value of the sub-band energy within the channel was located at , where and are the angular and radial quantum numbers, respectively. The energy value around the source edge within the same sub-band is denoted as . Furthermore, and represent the Fermi energies at the source and drain, respectively, while the voltage across, expressing the difference of the electrostatic potential between the source and drain in Fermi level , is . The electrostatic potentials in the transverse direction at the center and surface of the cylindrical channel are defined as and , respectively, and the conduction band edges at the same locations are given by and , respectively, as illustrated in Figure 2b. Correspondingly, the difference between and was described using an unknown parameter, . Moreover, the higher sub-band energy levels measured from are . According to our previous study, can be written in terms of along the z axis as follows [6,22]: where ℏ, , and are the Planck’s constant, the confinement electron effective mass in r axis, and the constant described as the -th zero point of the first-kind Bessel function with order , respectively, is the gate-source voltage, refers to the work function difference between the gate and the channel material, and is the conduction band edge at the flat-band condition corresponding to the relationship of . Moreover, is the perturbation matrix constant that was calculated in our previous work [24], denotes a geometric constant defined by , is the dielectric constant of the channel, and is the gate oxide capacitance per unit length. along the z axis is given by where A and B are the coefficients related to the voltage between the gate/drain and the source, respectively; is the scaling coefficient depending on mentioned above that can be given by

Figure 2

Rough sketch of the potential energy profiles along the channel and the transverse directions, and mechanisms governing the carrier transport in ballistic tunneling and thermionic modes. (a) Representation of energy-level distribution along the z direction at the channel center (). (b) Schematics of confinement potential energy distribution along the r direction at the barrier top () in the cross-section. The elementary charge is denoted by e.

Accordingly, the above coefficients are given by: where is the junction built-in potential between the source and channel that represents the electrostatic potential energy level of a sub-band at the source measured from [6]. In addition, the corresponding position to the barrier top of the sub-band energy level along the channel can be derived with:

2.3. Landauer Formula

Considering that the device channel length was shorter than 7 nm, electrons traveling through the device channel comprised the ones whose energies were higher or lower than the sub-band energy level at the barrier top. Thereby, the Landauer formula could determine the drain current contributed by the electrons in the ballistic transport between the source and drain with both the thermionic and tunneling modes. Then, the ballistic current between the drain and source in the channel can be written as follows [48]: where is the valley index of the energy band, is the transmission coefficient of electrons tunneling through the source-to-drain barrier in the sub-band with quantum numbers of and at a given , which is the electron kinetic energy in the z direction, and refers to the Fermi–Dirac distribution. Here, denotes the transmission probability derived from the WKB approximation: where is the electron effective mass along with the channel direction, and and are defined as the classical turning points at which is equal to the potential barrier energy.

3. Calculation Method

In this study, electrons carrying energy higher than the maximal potential barrier were transmitted from the source to the drain by thermionic emission. In comparison, electrons with lower energy can traverse the channel only by quantum tunneling through the potential barrier between the source and drain. However, it was assumed that electrons transporting from the drain could be neglected when calculating both the thermionic ballistic and the source-to-drain tunneling current calculations under the condition of high drain bias as , regardless of the electrons traveling over and through the potential energy barrier. Moreover, it was assumed that the ballistic current consists of electrons only excited to the lowest sub-band when and as the quantum numbers for all operating regions since no electrons are excited to the second or higher sub-band energy levels in the nanowire GAA MOSFET. Then, the carrier quantum reflection was considered to be negligible on the barrier in the channel. Nevertheless, in the channel of purely ballistic transport, the device behavior model could be simplified to describe the electrons’ transmission over and through the source-to-drain potential energy barrier in one consistent expression.

3.1. Parabolic Energy Level Approximation

Since WKB approximation was applied to calculate the transmission probability for the source-to-drain tunneling, it had to be evaluated by numerically integrating the classical turning points from Equation (12). Therefore, it could render the model analytically calculable by assuming that the sub-band energy profile along the device channel to be a parabolic function, which is valid in a short-channel device [47]. Then, when the channel length is sufficiently short, the lowest sub-band energy level along the z direction could be assumed to be approximately described by a parabolic function of z, denoted as [47]. This leads to the result below: where and are given by Equation (1) at the positions when and , respectively, as follows:

Then, by substituting and into Equation (3), and can be derived in terms of A and B, related to the bias conditions:

3.2. Unknown Parameter for All Operating Regions

This study introduces and implements a mathematical approach for establishing the smoothing function by connecting two kinds of expressions in their different regions. In addition, this method requires two different expressions of an arbitrary function and a constant in the expressions’ respective ranges of values. Compared with the previous work on the unknown parameter at the top of the potential barrier in all operating regions [6], an alternative approach to represent was achieved. Implementing the smoothing function to represent in Equation (18), it was set as zero for a large gate voltage bias, while its expression was valid in the subthreshold region incorporating the DIBL effect. It could thus be rewritten as for all operating regions with the following formula: where is a fitting constant to be adjusted for the best accuracy; this value was fixed to be for wide devices and bias conditions. To obtain a single valid expression for all operating regions, we blended with by following the same assumption as that in [24] as follows:

As a result, by substituting Equation (21) into Equation (1), we could rewrite the value of the lowest sub-band energy level at the barrier top for all operating regions as follows:

3.3. Potential Barrier Top for All Operating Regions

Furthermore, from analyzing Equation (8), the expression of was valid only in the subthreshold region and took meaningless values in the inversion region. According to the same inference in Section 3.2 to describe the smoothing function for a whole range of values, it could be assumed that , representing the distance between the locations of the source edge and the potential barrier top, decreased to nearly zero with the ever-increasing gate voltage. Thus, two expressions were needed to be given by Equation (8) and zero in the subthreshold and inversion regions, resectively. Then, the barrier top position in all operating regions was approximately evaluated as follows: where is a fitting parameter without any physical meaning; this value was fixed to be for wide devices and bias conditions to achieve reasonable accuracy.

3.4. Transmission Coefficient for All Operating Regions

According to Equation (12), the tunneling transmission coefficient through a sub-band energy barrier could be expressed in two conditions corresponding to a given . Then, representing the sub-band energy level profile as a parabolic function, the tunneling transmission coefficient applying the WKB approximation to substitute Equation (13) into Equation (12) was determined as follows:

On the other hand, the transmission coefficient of the energy sub-band here was unity and corresponding to the electrons with the total energies of and , respectively. Here, varied between 0 and 1. Considering that the smoothing function requires the expressions given by and one in their ranges of values taking with and , the transmission coefficient for all operating regions could then be represented as follows: where denotes the fitting parameter, and the value was set to be 1.5 to obtain a good agreement for all operating regions in the next section.

3.5. Drain Current for All Operating Regions

As illustrated in Figure 3, the transmission coefficient of the lowest sub-band here was unity and , corresponding to the electrons with energies higher or lower than the potential barrier top. For a high drain bias, the second term on the right-hand side of Equation (9) including could be neglected. Therefore, only the electrons injected from the source to drain could contribute to the ballistic drain current in this model. Since we could represent and in all operating regions using Equations (21) and (23), respectively, we demonstrate the drain current in all operating regions according to our previous work. Then, the formula of the drain current is as follows:

Figure 3

Schematic representation of electrons transported between the source and drain, corresponding to their total energy. Illustrations of the quantum tunneling and thermionic ballistic transport mechanism of electron transport across the potential barrier from the source to the drain along the channel. Sub-band energy profile obtained from Equation (1) is represented. Turning points ( and ) and are shown in literal values. Here, is the electron wave vector component in the z direction corresponding to the parabolic dispersion relation for . Schematic representation for the electrons transporting between the source and drain according to the total energy of electrons. Several possible cases for electron transport are demonstrated.

4. Results and Discussion

This section compares the results from our compact model and those from a TCAD simulation using SILVACO ATLAS [49]. Numerical computations were performed on the basis of the nonequilibrium Green’s function (NEGF) formalism under the effective mass approximation coupled with Poisson’s equation. Moreover, Wolfram Mathematica software was used to present the compact model. The cylindrical coordinate was used for quasi-three-dimensional calculation, no scattering was considered, and only the lowest sub-band was taken into account. The supply voltage was set to be 0.5 V following IRDS [50], and the value was extracted from the NEGF simulation.

As illustrated in Figure 4, the lowest energy level profiles between the NEGF simulation and our compact model using Equation (13) were compared. Three lines correspond to the results with 0, 0.2, and 0.5 V from the top. The channel length used in this compact model was empirically expanded into the source and the drain regions by 0.5 and 1 nm throughout the comparisons between the calculations of numerical NEGF simulation and the compact model considering Equation (13) for good agreement, respectively, referring to [6].

Figure 4

Comparison of the results of the lowest sub-band profile plot between the compact model using Equation (13) and NEGF simulation under different gate bias conditions.

Figure 5 shows the calculation results of using a solid black line corresponding to Equation (21) in determining the unknown parameter at the barrier top for all operating regions, denoted as . The results from and , corresponding to in the subthreshold and inversion regions, are also shown as the red and green dashed lines, respectively. When the DIBL effect was neglected, in the subthreshold region was 0. On the other hand, took a negative value in the subthreshold region, while it became 0 in the inversion region deriving from Equation (20). When 2.5 nm and 4 nm, , which should have been 0 in our mentioned assumption, took a negative value, even in the inversion region. Although this case could be fixed by resetting as 2 or a greater number, an inappropriate value of affected the accuracy of the calculations for , and characteristics, as shown in the following discussion.

Figure 5

Comparison of the results of , and .

Figure 6 demonstrates comparisons of the dependencies for -, corresponding to the potential barrier top positions in the device channel at different gate voltages resulting from the NEGF simulation and our compact model. The dashed lines, solid curves, and open circles show the calculations from the compact model by using Equations (8) and (23), and the NEGF simulation, respectively. The value of calculated from (8) was valid only in the subthreshold region, and took an impossible value in the inversion region. Thereby, the necessity of deriving Equation (23) from (8), describing as a non-negative value, speaks for itself. Moreover, when exceeded 0.5 V, the compact model could not capture the results of the barrier top positions from the NEGF simulation. This limited the valid range of the operating gate bias voltage in the compact model.

Figure 6

Gate bias dependence of the potential barrier top position along the channel calculated with the NEGF simulation (open circle) and our compact model using Equations (8) (dashed line) and (23) (solid curve).

Figure 7 illustrates the gate voltage dependence of the lowest sub-band energy level at the barrier top using Equation (22). The comparisons between the NEGF simulation and our compact model with the lowest sub-band are demonstrated as open circles and solid lines, respectively, maintaining good consistency.

Figure 7

Comparison of the results of the lowest sub-band energy level at the barrier top between NEGF simulation (open circle) and our compact model (solid line) for the lowest sub-band.

Figure 8 shows comparisons of dependencies for -, corresponding to electrons having the total energy of , between our numerical compact model using Equation (27) and NEGF simulator calculations. The open circles and solid lines relate to the calculations from the compact model and the NEGF simulation, respectively, demonstrating reasonable agreement between these results at each gate voltage level, from 0 to 0.4 V in 0.1 V steps.

Figure 8

Comparison of the transmission coefficient calculations results between the NEGF simulation (solid line) and our numerical compact model (open circle) for the lowest sub-band.

Figure 9 illustrates the drain current versus gate voltage characteristics calculated at various sets of the channel lengths and wire radius by considering Equation (28) in all operating regions. The open circle and solid line represent the comparison results of the compact model and NEGF simulation calculation, respectively. As shown in this figure, our compact model could capture the drain current increase caused by the source-to-drain tunneling in the subthreshold region while maintaining acceptable accuracy in on-current characteristics. This model was not valid when the sub-band energy level at the barrier top was lower than that at the source, which determined the maximal value of the gate voltage bias to sustain the validity of the model by solving using Equations (16) and (17). In our computations, these gate bias voltages were larger than the given supply voltage. Nevertheless, in order to obtain reasonable agreement for all operating regions, the diameter and length values of the channel needed to meet a ratio of less than unity ( < 1) in this study through their comparisons between the NEGF simulator and our compact model. Furthermore, the tunneling current should also be considered in the inversion region.

Figure 9

Comparison results between the compact model (open circle) considering the Landauer formula and NEGF simulation (solid line) for drain current characteristic calculations in all operating regions.

5. Conclusions

This paper presented a new approach to express the numerical compact model of the drain current in the ultrashort channel nanowire GAA MOSFET using the Landauer formula. Assuming that a parabolic function could approximately determine the potential distribution along the ultrashort channel, the WKB approximation was capable of evaluating the transmission coefficients for source-to-drain tunneling. In addition, it could approximately evaluate the position of the potential barrier top and the transmission coefficient for all operating regions by using smoothing functions. Lastly, the drain current characteristics could be calculated for all operating regions. This numerical compact model was tested with a NEGF simulation, demonstrating reasonable agreement in all operating regions. Although the compact model contained several empirical parameters, they were fixed throughout the comparisons. However, a compact model with fully analytic expressions is expected to be developed to carry out a circuit simulation in a future study.