A review of quantum transport in field-effect transistors

While the Lilienfeld patent is generally considered to be the start of the FET, it is not all that clear from the description in the patent, and it was a later patent by Heil that more clearly described the surface-oriented device [4]. It is generally believed that the work at Bell Telephone Laboratories was attempting to make the device in Ge, but the lack of a stable oxide led them to the point contact transistor (two metallic points attached to the Ge layer) [5], developed by Bardeen and Brattain [6, 7]. Shockley would shortly follow this with the junction transistor (diffused or grown layers of various doping) [8].

Bringing the FET to reality would take somewhat longer. Atalla stabilized the Si surface with a thermal oxide [9], and this led he and Kahng to develop a MOSFET around 1960. Atalla got irritated at Bell's disinterest in the MOS technology and left in 1962; perhaps this is the reason why Kahng is the sole author of the patent [10]. Fortuitously, others in California recognized the advantages of MOS technology that eventually led to its use in integrated circuits through two patents filed in 1959 [11, 12].

Complementary n-channel and p-channel devices appeared in 1962 [13, 14], as low power devices which eventually dominated the continued development of integrated circuits and became the well-known CMOS. The advantage possessed by these devices lay in their planar layout (the plane of the current flow is in the layer plane, whereas it was normal to the plane in a junction bipolar transistor), which was much favored for large scale integration, and their low power dissipation. Growth in the technology is largely measured by Moore's law; the idea that the dimensions would decrease from one generation to the next with the result that the number of transistors on a single chip would increase exponentially [15], although this law is basically an economic law rather than a physical law [16, 17]. But, these integrated circuits were, by and large, digital, which means that the transistors were switching transistors. Such a switching transistor needs well defined stable states, which define logic levels, and these are given by the ground and bias voltage. Transition between these two states is not of much concern, as long as it happens quickly and reliably.

Other requirements exist for FETs used in amplifiers. These include a desire for a linear behavior in small, or large, signals, and FETs have provided a great many analog applications, especially in microwave applications. Often these are not Si-based devices, but are III–V devices, such as the high-electron-mobility transistor (HEMT), which still is a FET. The highest frequency performance to date have come from InP MOSFETs (metal-oxide-semiconductor FET [18]) and InP-based, InAs channel HEMTs [19]. The development of GaN-based devices has led to excellent power HEMTs [20].

The structure of the MOSFET has undergone a great deal of change since the first MOS integrated circuit. In particular, the down-sizing of the device has led to the introduction of new structures, new materials, and new configurations. This evolution has occurred in both the switching devices and the analog devices. For example, as the size of the switching devices has decreased, a growing problem appeared as the mobility continued to drop. This drop was due to the effect of surface-roughness scattering at the oxide-semiconductor interface [21], as well as an increase in the impurity scattering due to the increased doping in the substrates. To overcome this decrease in mobility, the first change was to introduce strain in the channel [22]. The problem lay in the fact that one wanted tensile strain in the n-channel in order to split the six-fold degenerate conduction band valleys, and compressive strain in the p-channel in order to warp the valence band, both of which effectively lowered the effective mass and increased the mobility [23]. This was accomplished by putting a SiN layer over the gate stack in the n-channel devices, and using a SiGe alloy for the source and drain regions in the p-channel devices. Eventually, one had to go further and address the problem of the oxide becoming too thin, which was addressed through the use of a new gate structure featuring HfO₂ [24]. Finally, continued evolution forced one to turn the normally horizontal FET on its edge, leading to the FinFET [25], which became mainstream around 2011. What the future holds is open for discussion, but gate-all-around quantum wire FETs seem to be one way forward [26, 27], particularly with IBM's announcement of 2 nm nanosheet technology in 2021.

As the FET progresses over time, it is worthwhile to discuss what properties are needed in the evolving device. It was already mentioned that switching devices need to be different than analog devices. For example, the switching device needs to make the on-to-off transition, and the off-to-on transition rapidly and the linearity of the transition between these two points is not particularly relevant. What is needed is a low voltage in the on state and an extremely low current in the off state. Leakage current in the off state can create significant problems with heat dissipation in the circuit [28]. Even an off current of 1 nA is too large for a circuit with 10¹⁰ transistors! The transition from the planar FET to the FinFET was due to the need to control this off current by being able to pinch the device off with two opposing voltages from either side of the fin [25]. Similarly, the likely transition to gate-all-around quantum wires will be for a similar reason.

In the analog device, whether it be in Si or GaAs, the transition between on and off states is not nearly as important as the linearity of the gain over the range of voltages that will be input to the device. Here, the device is nearly always on, and biased to provide as large a range of voltage swing as possible while yielding as linear a gain as possible. These are distinctly different requirements than those for the switching transistor. Quite often, this is further complicated in large signal devices by the need to deliver significant power, so that the gain compression at high power is also a strong consideration in design. This means that the microwave power device design is quite different from the low power linear amplifier design [18], and mixed-mode devices become even more difficult [29].

To aid in the design of FETs, there are a great many commercial software packages available, and many private simulation packages within the scientific community. In the early years, most evolution within the world of fabricating FETs did not rely upon such packages, but they have become important as the device size has continued to be reduced, not least because of the increased role of parasitics to the design. That is, through most of the evolution of Moore's law, transistors have been downsized through the use of a strict scaling relationship [30]. In this scaling, the electrostatics are maintained. All dimensions are reduced by the same factor and voltages and doping densities are adjusted to maintain the electrostatics. But, this approach ignores the effects of parasitics, which become more important the smaller the device becomes and the closer other devices appear. For example, while the electric fields can be scaled, capacitances begin to be dominated by edge effects and so do not continue to scale properly. Further, the device can no longer really be considered in isolation; it sits in an environment of other devices and circuit elements. This background provides an environment from which the device cannot be isolated [31]. Whether this is interface roughness scattering [32, 33], which is of course already well known, remote phonon scattering [34], or device-device interactions arising from unintended coupling between devices [35], it is not fully apparent how much these effects, particularly the latter, have been included within widely used simulation and modeling packages.

There is also the question of granularity, particularly with respect to dopants. If one considers a small device in a semiconductor region 20 nm long, 20 nm high, and 5 nm thick, doped to 10¹⁸ cm⁻³, there are on average only two dopants in the volume. And, for a FET, it critically depends exactly where those dopants are located. They are far more effective when close to the source than when they are close to the drain. Even with more impurities, there is the problem that an electron can be interacting with several impurities at one time [36]. Difficulty here can arise from the fact that it might be difficult to compute an accurate screening function, or even to establish that there is any time between collisions; multiple scattering can dominate [37]. While this is primarily with impurity scattering, other factors, such as rough interface scattering can also be affected by this granularity; basically an inhomogeneity that spreads throughout the device. This can be especially problematic when it couples with quantum effects.

There have been a great many discussions about the role of quantum effects and quantum transport in FETs. However, the basic idea boils down to the problem of just how big is a quantum electron (or hole). Classically, the size of the electron has been pondered over the years, but the consensus may put it at ∼10⁻¹⁸ m. But this will not work for a quantum particle, and the size of the electron depends upon how big a wave packet representation of the electron is quantum mechanically [38]. Several considerations of this suggest the size is in the range of 3–7 nm. This is quite significant when the idea of the 5 nm node (for integrated circuits according to Moore's law) is discussed. It already becomes a problem in a simple MOSFET, for how can an ∼5 nm electron be stuffed into the classical size of the inversion layer, which is a ∼1 nm potential well? In fact, it just cannot be done and a self-consistent solution of Poisson's and Schrödinger's equations must be pursued [39].

Pursuing the self-consistent potential mentioned for the case of electrons in an n-channel Si MOSFET leads to the splitting of the six-fold degenerate valleys of the conduction band. The two valleys whose heavy mass direction is normal to the oxide-semiconductor interface form one set of quantum levels which are lower in energy than the levels of the other four valleys that have the light mass normal to the interface. From such calculations, nearly all MOSFETs will show this quantization. But, the quantization is not observed in every day operation for several reasons. First, the transport direction is parallel to the interface, not normal to it, so the quantization is a second-order effect. Secondly, such a calculation is a one electron effect appropriate for having no inversion density present. When an inversion density is present, there is another quantum effect, the many-body interaction [40]. While the confinement raises the quantum levels, the self-energy of the carriers lowers these energy levels. Thus, a full consideration is quite complicated. Nevertheless, the use of tensile strain to increase the separation between the two-fold valleys and the four-fold valleys is a common situation in today's devices, and the mobility is higher in the two-fold (and lower energy) valleys.

While not apparent, many modern simulation packages account for the size of the electron; not directly, but through clever manipulation. The total energy of the electron (or hole) population is a summation over the density, weighted by the potential at each point. Mathematically, the assumed Gaussian shape of the electron wave packet can be moved from the electron to the potential [38]. This leads to what is called the effective potential, which greatly simplifies simulations of these devices [41]. It is the effective potential that appears in these packages, but it arises from the electron's size.

From these considerations, it becomes apparent that quantization can set in whenever confinement effects become comparable to the electron (or hole) wave packet. In today's ultrasmall devices, this means that quantum effects play a considerable role. A further example of this appears in the FinFETs that are popular at this time. The FinFET turns the channel on its side, so that the inversion layer can extend up either side of the fin, and even over the top. With the gate potential on either side of the fin, this makes turnoff more effective by these two potentials working together. However, the normal state in which there are two inversion layers, one on either side of the fin, can change into a single inversion layer in the bulk of the fin. When the fin is sufficiently thin and the density is not too high, this central inversion layer is the preferred state [42]. Thus, the quantum effects can change the basic nature of the FET, from a surface-oriented device to a bulk-oriented device, which has the added benefit of less scattering from the rough interfaces.

Another effect which has been around for some time is tunneling, a true quantum effect. Concern in the past has dealt with tunneling through the gate oxide. In the presence of high electric fields across the gate oxide, the barrier between the gate metal and the FET channel is no longer the simple rectangular barrier of the textbooks, but becomes trapezoidal or triangular in nature [43]. Such tunneling has been known to be a problem in thin dielectrics almost as long as we have had MOSFETs [44]. What has become more interesting with the continued down-sizing of the device, is direct tunneling between the source and drain [45]. While this early work observed resonant tunneling between source and drain, likely due to a localized state, the trend of ever smaller devices has led to direct tunneling between the two. The adoption of high-κ dielectrics was a direct response to the problem of gate tunneling in the very thin oxides required in the scaled devices. It is not clear as yet how the problem of source-drain tunneling will be addressed, or solved.

Nevertheless, it is clear that quantum effects are appearing in ever greater quantity as we continue to evolve with Moore's law. Moreover, it is not at all evident that they can be treated in isolation, one at a time, as quantum mechanics is notoriously a nonlocal theory, so that we may expect many of the quantum effects to interact with each other as well as with our understanding of how the FET operates.

Probably, almost everyone has heard of Feynman's quote '...nobody understands quantum mechanics...' [46]. If one is trying to explain quantum mechanics in Bohr's view that nothing exists until it is measured, then you can understand Feynman's comment. Bohr was very interested in his desires to do away with causality and determinism [47]. But, there are alternatives, and most engineers and device scientists tend to follow these alternatives whether they are aware of it or not [48]. It is a form of quantum mechanics in which determinism exists, particles are real and follow real trajectories, which is what makes semiconductor devices work. The difference then, between quantum transport and classical transport, lies in that '...additional quantum potential...' described by Kennard [2]. We know that this quantum potential derives from the wave function itself, and introduces both interference and coherence effects into the transport. The interference also encompasses entanglement, that magical force unique to quantum mechanics and essential for quantum computing [49].

The quantum potential has a history in the hydrodynamic expansion of the Schrödinger equation, first done by Madelung [50] and Kennard, and much later by Bohm [51]. Like the effective potential, the quantum potential already has been used in device simulation [52], and appears in a number of commercial simulation packages. This provides a definite path forward. However, again like the effective potential, it is not a full correction, as it does not account for all the nonlocal and interference interactions. More extensive quantum transport equations must be used, with the result that a great deal more complexity enters the problem [53]. And, this brings us to this review, in which the more extensive transport equations, and their use in simulation/modeling of FETs, is the heart of the topic.

One example of this is transport in the presence of very weak scattering — quasi-ballistic transport. Classically, the electrons do not behave any differently from a strong scattering regime. However, quantum mechanically, the electron can interfere with itself as it passes a single impurity. This interference can affect the resulting self-consistent potential [54]. The need to properly handle such interferences requires the use of quantum transport approaches.

In section 2, the nature of quantization with the FET will be discussed. Principally, this begins with a more in depth discussion of some of the topics already mentioned in section 1.2, but goes on to thoroughly address the concepts of interference and coherence, especially those that arise in the physically short channels that are emerging today.

In section 3, the nature of quantum transport, its equations, and its philosophy will be developed. This will extend from the simple approaches already mentioned to the complexities of Wigner functions and non-equilibrium Green's functions. The principle concept of their difficulties is discussed as well, and how one addresses dissipation and decoherence, e.g. the role of scattering.

Finally, in section 4, some concluding remarks and summary will be presented. In addition, some thoughts about the future will be presented.

In a classical treatment of the inversion layer in a standard MOSFET, the carrier density peaks at the interface between Si and the oxide, and then decays exponentially with the surface potential away from the surface. This behavior is completely opposite from that expected in the quantum treatment. While the positive gate voltage attracts the charge to the surface (an n-channel device is considered here), the nature of the quantum effect is that the charge has to be nearly zero at the oxide interface. The local charge is usually treated as $\rho \left( z \right) = - e{\left| {\psi \left( z \right)} \right|^2}$, where $\psi \left( z \right)$ is the wave function. Since, the wave function must nearly vanish due to the large offset potential of the oxide, this means the charge must similarly vanish. In addition, the potential formed between the oxide potential and the band bending within the semiconductor leads to quantization of the motion normal to the surface. This forms one or more sub-bands, which are quasi-two-dimensional energy-momentum relationships for motion parallel to the interface. The actual shape of this potential, and the corresponding wave functions must be found self-consistently by solving both Schrödinger's equation and Poisson's equation [39].

An example of the quantization is shown in figure 1, where the potential and two sub-band energies are shown [61]. The self-consistent potential is labeled as V_eff(z). The thick curves include the quantum exchange-correlation self-energy correction, while the thin curves do not include this correction. Here, it was assumed that the p-type substrate was doped to 2.8 × 10¹⁵ cm⁻³ and the inversion density was 4 × 10¹² cm⁻². The results are obtained by including the normal Hartree potential found from the Poisson equation, an exchange-correlation correction, and an image force [39]. The exchange-correlation potential is a many-body correction accounting for the carrier-carrier interaction, and the form used was taken from Hedin and Lundqvist [62]. One can see that the wave function peaks approximately 1.5–2.0 nm from the surface and then decays.

Zoom In Zoom Out Reset image size

The two major observables from the role of quantization is the set back of the peak in the wave function, and thus the charge density, from the surface and the quantization energy of the lowest sub-band, which appears to be about 45 meV, from the minimum of the conduction band. These two changes affect the gate capacitance of the MOSFET [63].

An important impact of the quantization is that the energy levels depend inversely upon the effective mass of the carriers normal to the interface in each valley. Hence, in a (001) surface of the Si, carriers in the pair of (001) valleys exhibit the heavy ∼0.98m₀ longitudinal mass, while the other four valleys exhibit the light ∼0.19m₀ transverse mass. This leads to splitting of the six conduction band valleys, with the four transverse valleys sitting higher in energy. Thus, the sub-bands shown in figure 1 are those for the two valleys with their major axis in the (001) direction. As these two valleys exhibit the lighter mass in the transverse directions (down the channel), this will lead to higher mobility. This effect will be exaggerated when strain is introduced, as discussed below.

This quantization has been known for a great many years [39]. Over the intervening decades, there has been a search for methods of treating the quantum effects without having to solve Schrödinger's equation. As remarked above, one approach is to use an effective quantum potential [2]. If the wave function is written as $\psi \left( r \right) = A\left( r \right){\text{exp}}\left( {iS\left( r \right)/\hbar } \right)$, then it has been suggested that the Bohm quantum potential is of the form [2, 50, 51]

$$\begin{align}{V_{QB}} = - \frac{{{\hbar ^2}}}{{2mA}}\frac{{{\partial ^2}A}}{{\partial {r^2}}},\end{align} \tag{ 1 }$$

where r is the direction in which the force is to be determined. If we use the previous relationship above between the density and the wave function, then this equation can be rewritten as (this follows directly with the substitution of the amplitude written as $A = \sqrt n$):

$$\begin{align}{V_{QB}} = - \frac{{{\hbar ^2}}}{{2m\sqrt n }}\frac{{{\partial ^2}\sqrt n }}{{\partial {r^2}}},\end{align} \tag{ 2 }$$

where n is the local density. This form is often called the density-gradient potential, and appears to have been first used in FET simulations by Grubin and Kreskovsky [64]. The compact form of this quantum potential means that it can be incorporated into a wide range of transport models for devices of many sorts, including MESFETs (metal-semiconductor FET [52]) and the MOSFET [65–67].

Wigner also introduced a form of effective potential, in his study of the impact of quantum mechanics on thermodynamics [68]. He found a pair of correction terms to the kinetic energy E in the form of a potential given by

$$\begin{align}{V_{QW}} = - \frac{{{\hbar ^2}}}{{12m}}{\nabla ^2}E + \frac{{{\hbar ^2}}}{{24m}}{\left( {\nabla E} \right)^2}.\end{align} \tag{ 3 }$$

If the normal thermodynamic form

$$\begin{align}n\sim {{\text{e}}^{ - E/{k_{\text{B}}}T}},\end{align} \tag{ 4 }$$

where k_B T is the thermal energy, is used, (3) becomes

$$\begin{align}{V_{QW}} = - \frac{{{\hbar ^2}}}{{8m}}{\nabla ^2}\left[ {{\text{ln}}\left( {n/{n_0}} \right)} \right],\end{align} \tag{ 5 }$$

where ${n_0}$ is a reference density. With a little manipulation, this latter form can be shown to be only a factor of 2 different from the density-gradient form (2) [69].

Still another form of quantum potential modification arose from the study of the non-zero size of an electron in a semiconductor [38]. In this case, it was shown that by considering the total energy of a potential (weighted by the density) in the Hamiltonian, the shape of the electron wave packet could be moved to the potential arising from Poisson's equation, which led to it being smoothed by [69]

$$\begin{align} {V_{QF}}& = \frac{1}{{\sqrt {2\pi } \alpha }}\int {{\text{d}}r^{\prime} } V\left( {r^{\prime} } \right){\text{exp}}\left( { - \frac{{{{\left| {r - r^{\prime} } \right|}^2}}}{{{\alpha ^2}}}} \right) \nonumber\\ & = \frac{1}{{\sqrt {2\pi } \alpha }}\int {\text{d}} r^{\prime} V\left( {r - r^{\prime} } \right){\text{exp}}\left( { - \frac{{{{\left| {r^{\prime} } \right|}^2}}}{{{\alpha ^2}}}} \right), \end{align} \tag{ 6 }$$

where [70]

$$\begin{align}{\alpha ^2} = \frac{{{\hbar ^2}}}{{2m{k_{\text{B}}}T}},\end{align} \tag{ 7 }$$

and the last line of (6) arises following a simple change of variables. As mentioned, this smooths the potential, especially in the region where the quantum well is located, and provides both the charge set-back and the quantization energy of the lowest subband [41]. This is shown in figure 2. It also can be shown to be equivalent to the above quantum potentials. To do this, the potential is expanded in a Taylor series around the position r, as follows [69]:

$$\begin{align} {V_{QF}}& = \frac{1}{{\sqrt {2\pi } \alpha }}\int {\text{d}} r{\prime} \left[ V \left( r \right) + r{\prime} \frac{{\partial V}}{{\partial r}} + \frac{{r{{\prime} ^2}}}{2}\frac{{{\partial ^2}V}}{{\partial {r^2}}} + \ldots \right]\nonumber\\ & \quad \times {\text{exp}}\left( { - \frac{{{{\left| {r{\prime} } \right|}^2}}}{{{\alpha ^2}}}} \right). \end{align} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

The first term is just the normal potential, and the next term vanishes upon integration, so that the leading term is the quantum correction which gives, in three dimensions

$$\begin{align}{V_{QF}} = \frac{{{\hbar ^2}}}{{2m}}\frac{{{\partial ^2}V}}{{\partial {r^2}}} = \frac{{{\hbar ^2}}}{{2m}}\frac{{{\partial ^2}{\text{ln}}\left( {\frac{n}{{{n_0}}}} \right)}}{{\partial {r^2}}},\end{align} \tag{ 9 }$$

which is again within a constant of (5) when using (4). The point is that almost all of the various quantum potentials wind up being of the same form for small quantum effects, that is for small nonlocality. In figure 3, an example is shown comparing the Bohm potential (2) and the total effective potential (6) [71]. The structure is similar to a MESFET, in that the short barrier that creates the quantum point contact at the bottom of the structure is similar to a very short gate. In panel (a), the Bohm potential is created after solving for the full quantum waves in the structure in the absence of a self-consistent potential, and particle transport through the structure is indicated by the light stream-lines. In panel (b), both the Schrödinger and Poisson equations are solved self-consistently, and the Bohm potential is found, and particle transport is then determined. In panel (c), only the Poisson equation is used for the self-consistent potential, and the effective potential is computed from (6) and used to study the particle motion. In this latter method, there are far more vortices beginning to form in the transport, as there is more confinement observable in this method, as it is non-perturbative.

Zoom In Zoom Out Reset image size

When strain is added to a device, with the intent to modify the transport properties by modifying the band structure, it becomes a quantum effect [22, 23]. In the case of the n-channel MOSFET, the strain is tensile and often accomplished by putting a Si₃N₄ layer over the gate stack. This tries to stretch the channel. It was remarked in the previous section that quantization in the inversion layer separated the six-fold conduction band into a two-fold set of elliptical energy surfaces whose major axis is normal to the interface and a four-fold set whose minor axis is normal to the interface. Unfortunately, this separation is not particularly large, especially in elevated temperatures. The use of tensile strain increases this separation so that the conduction is dominated by the two-fold pair of sub-bands. A 1% uniaxial strain along the [110] transport direction can give as much as 70 meV splitting between the two sets of valleys [72]. Here, the effective mass for transport in the channel becomes the transverse mass, (∼0.19m₀ in Si). This results in an approximately 40% reduction in the transport mass relative to the normal conductivity mass in Si, although the strain will affect this mass [73]. In addition, the carrier scattering between the two-fold and four-fold sets of valleys is greatly reduced, and these two effects together lead to a much higher mobility for the carriers in the inversion channel. Because these valleys are anisotropic, the actual motion in the channel can be quite complicated [74].

In figure 4, the effect of strain on the mobility is plotted for an n-channel device on the (001) surface of Si with the uniaxial strain along the [110] direction [23]. The strain deforms the normally circularly symmetric shape of the two-fold valleys, as can be seen in the figure. Several different levels of strain are given by different curves, and one can see that the mobility, though larger, becomes anisotropic in the valleys. Here, the effect of strain on the bands was calculated with an empirical pseudopotential method, and the transport was simulated with an ensemble Monte Carlo method. The results are comparable to what is seen experimentally in such devices. The result of warping of the ellipsoidal energy surfaces is quite expected in the understanding of these devices. The elongation of the normally circular energy surface in the plane of the device, causes a decrease in the effective mass in the channel direction.

Zoom In Zoom Out Reset image size

Note that figure 4 actually plots the mobility as a function of the direction, and not the energy surfaces. Strain along a [110] direction elongates the crystal in that direction. Since momentum is proportional to the inverse of the lattice constant, this compresses the normal energy circle, so that the long axis of the constant energy line would lie along $\left[ {\bar 110} \right]$, which is normal to the strain. Then, if the resulting constant energy line is written as

$$\begin{align}E = \frac{{{\hbar ^2}}}{2}\left( {\frac{{k_{\left[ {110} \right]}^2}}{{{m_{\left[ {110} \right]}}}} + \frac{{k_{\left[ {\bar 110} \right]}^2}}{{{m_{\left[ {\bar 110} \right]}}}}} \right),\end{align} \tag{ 10 }$$

the fact that the major axis in the ellipse is normal to the strain direction implies that

$$\begin{align}{m_{\left[ {\bar 110} \right]}} > {m_{\left[ {110} \right]}}.\end{align} \tag{ 11 }$$

As a result, the mass in the strain direction is reduced, while the normal direction has an increased mass, a result in keeping with the theoretical calculations [23]. This leads to a larger mobility in the strain direction, as is indicated in figure 4.

While the two valleys in the lowest sub-band are normally degenerate, quantum interactions between carriers in these two valleys can lead to their splitting, which has been observed experimentally [75]. It has been found that the strain applied along the [110] direction, and the warping of the energy surfaces in the plane can lead to significant splitting of the these two valleys, while maintaining a coupling. The resulting energy surfaces are quite dramatic, as shown in figure 5 [76]. Here, the structure was computed using a $k \cdot p$ method, including the strain, and the dispersion is seen to become very non-parabolic [77]. Another important aspect is that the effects can become stronger in ultrathin layers of Si, such as in SOI or nanosheets [78, 79].

Zoom In Zoom Out Reset image size

More than 30 years have passed since SiGe alloys were suggested for use in semiconductor devices [80]. If the SiGe alloy is relaxed, then a subsequent thin Si layer will be under tensile strain. This was thought to be a good method of improving the transport [81]. This would lead to its use in hole transport.

In opposition to the strained Si above, growing a strained SiGe layer on Si would put this layer in compressive stress, and using a subsequent Si layer on top would bury the SiGe layer as a quantum well [82]. The top Si cap layer moves the holes away from the surface as a means to reduce the surface-roughness scattering, while the strain in the SiGe layer splits the light and heavy-hole bands at Γ, providing a reduced transport effective mass. A drawback is that the SiGe layer has significant alloy scattering which tends to lower the mobility, so that there are competing effects. To generally describe the valence band structure, in both the quantum effects and with strain, requires a form of full-band calculation for the complicated splitting, warping and crossing of the various valence bands [83].

The approach to p-channel devices more or less settled upon the use of SiGe in the source and drain regions, where the larger lattice constant put the channel under compressive strain [22]. This uniaxial compression raises the light-hole band extremum above the heavy-hole band, so that the effective mass of holes becomes much smaller and the mobility is raised. This creates a problem, in that the two bands will now cross at a not too high value of carrier energy. They do actually cross, but hybridize in a manner to avoid this crossing; the result nevertheless is a very non-parabolic behavior with the initial light mass getting very large as the carrier energy increases. To understand this better, the normal unstrained valence bands in Si are shown in figure 6, along with an equal energy surface for the two top bands [84]. One can understand the phrase 'warping' from the two constant energy surfaces. The maximal extension of the heavy-hole surface is along the (111) directions, while that of the light-hole surface is along the (100) directions.

Zoom In Zoom Out Reset image size

The energy surfaces in figure 6 are at very low energy. As the energy increases, the warping becomes much more pronounced, as can be seen in figure 7 for the unstressed bands (the upper row is the upper band after hybridization while the lower row is for the lower band). As the channels are normally in the [110] direction, the holes will see a much reduced effective mass and therefore an enhanced mobility.

Zoom In Zoom Out Reset image size

The effect of the strain and the gate field can be observed in figure 8, where it can be seen that the gate field really has little effect on the constant energy surface for the dominant band [85]. Here, the lowest energy sub-bands are shown for both the strained and unstrained energy surfaces, as calculated with a $k \cdot p$ method. As the surfaces are not affected by the gate field, the masses similarly are not affected. Similar calculations for the valence band structure in the presence of the uniaxial strain along [110] have been done by Shifren et al [86], Wang et al [87], and Kotlyar et al [88], among others.

Zoom In Zoom Out Reset image size

The problem with fluctuations in the number of impurities under the gate was discussed already by Keyes [28] half a century ago. The actual distribution of the impurities cannot be controlled, and fluctuations in the actual number lead to fluctuations in the threshold voltages. In large devices, this fluctuation is small, but in small devices, it becomes much more significant. For example, if we consider a 10 × 10 × 10 nm region with a doping of 10¹⁹ cm⁻³, there are only 10 ± 3 dopant atoms in this volume. Hence, as devices become small, the fluctuations become large (30% in this example).

Attention was called to this again when Wong and Taur included the atomistic nature of these impurities in their MOSFET simulations [89]. The importance of this was quickly realized and picked up in the simulations of MESFETs and HEMTs [90, 91]. In figure 9, simulations of 16 different MESFETs of nominally the same dimensions and doping profiles and 18 HEMTs are shown, where the drain current is plotted as a function of the actual number of dopants under the gate, although the doping level is kept constant. It can be seen that the variation of the current with dopant number is not monotonic and fluctuates dramatically.

Zoom In Zoom Out Reset image size

The inclusion of the actual positions of the random doping distributions was adopted widely in the simulation of MOSFETs as well [92–95]. In figure 10, the variation caused by the atomistic doping is illustrated by plotting the drain current as a function of gate voltage for a great many different implementations of a 0.1 μm gate length MOSFET, doped to 8 × 10¹⁷ in the substrate [93]. The fact that the current varies by more than an order of magnitude in the sub-threshold region points to the importance of this effect. It also became evident that impurities near the source-gate entry point into the channel affected the current far more than impurities further from this region [96].

Zoom In Zoom Out Reset image size

The rise in these fluctuations has led to the ideas that the best course in past few years has been to not add any doping in the substrate, but let the channel appear similar to a n⁺ –i–n⁺ structure. With operation below 1 V, adequate voltage separation from the substrate can be achieved especially with double-gate FinFET type devices. A name of junctionless devices was coined some time ago [97].

While the fluctuations are important, it should also be recognized that an individual impurity can have a significant effect on the current due to quantum transport behavior. In figure 11, simulations of the behavior of a single [98] electron nearing a single impurity in the MOSFET channel are shown for various times. Here, the electron is represented as Gaussian wave packet (the size of the electron and these wave packets were discussed above in section 2.1). The right hand column is for a classical Boltzmann solution. The wave packet expands somewhat, which is normal for Gaussian packets, and moves deeper into the channel as it is repelled from the impurity. On the left hand panel though, the wave function becomes quite broken up due to the quantum interference, that leads to decoherence of the packet.

Zoom In Zoom Out Reset image size

But, the impurity can have a larger effect. The reader should be familiar with Young's two-slit optical experiment [99]. It can be replicated with electrons. A charged wire in a transmission electron microscope (TEM) creates what is known as a bi-prism, and electrons split as they pass this wire and form an interference pattern [100] just as that of photons in the Young experiment. Importantly, a single electron in the TEM, or a single photon [101], will lead to the interference pattern. That is, the single electron (or photon) must have both wave and particle properties, as recognized by Einstein [102]. For the purpose here, these observations mean that a single electron passing by a single impurity can create an interference pattern representing its wave-like self-interactions in diffracting around the impurity.

Two electrons and two impurities provide a complicated wave interaction process that illustrates the interference principle. Consider a simple nanowire structure with two repulsive impurities embedded in it, as shown in figure 12 [54]. Here, the nanowires are 40 nm wide, and the area of interest is some 60 nm long. The impurities are located at the two green circles, which are constant energy lines. The electrons are continually injected at y = 0, x = 20 nm, and absorbed completely at y = 60 nm. Consecutive injections are considered as independent, identically distributed statistical experiments, giving rise to the stationary distribution of a single electron. These electrons are simulated by a Wigner function wave packet (Wigner functions are covered in part 3 below). The continual injection allows the study of the steady-state inference with impurities. The electrons are injected at an energy of 0.14 meV, which is above the second transverse sub-band in the nanowire. This, plus the nonlocal interaction with the impurities leads to the injected carriers forming the two beams apparent in the figure in the region before the impurities (the carriers travel from bottom to top in the figure). The carrier 'density' from ${\left| {\psi \left( {x,y} \right)} \right|^2}$ is plotted in the figure. It is clear that the electron is diffracted by the impurities, and then interferes with itself downstream from them. This would be present even if only a single impurity were present. The interference peak evolves all the wave to the exit of the nanowire. These interference peaks are certainly reminiscent of the two-slit experiment.

Zoom In Zoom Out Reset image size

An analytical study of similar scattering with a single impurity has been done by Barker [103]. In this work, it is clear that the interference arises from the matrix element itself for the scattering interaction. Yet, the result shows almost the same interference pattern as figure 12, but with a few amplitude differences. For example, in Barker's work, the amplitudes of the different beams decay from the peak for the forward beam. In figure 12, the second beam to either side is weaker and this is likely the result of having two impurities in which their second beams tend to interact negatively to reduce this amplitude. Nevertheless, it is apparent from these, and other sources, that the quantum interference around an impurity is a real event that will influence nano-scale devices. Indeed, a study of the electron wave on its own when near its donor atom shows a quite complicated wave function that is dependent upon the presence of strain in the system [104].

The interference that appears with one or a few impurities becomes much more complex when many impurities are present [55]. The impurities lead to a random potential, which if sufficiently large can localize a number of the band states [105]. But, the random potential has another effect. Transport through the random potential is quite sensitive to the position of the Fermi level and the presence of any magnetic field, even the self-magnetic field of the current in the FET channel. Small variations in these quantities can lead to large variations in the conductance through the channel, and even chaotic behavior. These fluctuations are typically referred to as universal conductance fluctuations, and have been observed in a large Si MOSFET at low temperatures [106]. In small modern devices, this can lead to significant current fluctuations in the drain characteristics if there is insufficient scattering in the channel [107]. This will be seen below in figure 15(a).

Even when there is adequate scattering, current filaments can form in the channel [108], as shown in figure 13 for an n-channel, 50 nm MOSFET. In the figure, the channel runs from 50 to 100 nm, and the current is in A cm⁻². The doping in the channel was 5 × 10¹⁷ cm⁻², while the source and drain were doped to 2 × 10¹⁹ cm⁻². The individual dopants are shown as dots or small circles, while the size of the circle corresponds to the depth below the surface of the impurity. The impurities are treated in real space by a molecular dynamics approach in which the impurity potential is split into a long-range part which appears in the Poisson equation and a short-range part treated by the direct molecular dynamics forces between the electrons and impurities [109]. An interesting effect is that the current filaments do not really break up and decohere until almost 40 nm into the drain region!

Zoom In Zoom Out Reset image size

These effects are often grouped under the heading of excess noise in experimental studies [110, 111]. As the localized states can be interpreted as traps, the fluctuations often appear as trapping/detrapping effects [112–114]. And, often the fluctuations are coupled to those from the random interface potential due to surface roughness [115].

One of the oldest problems in FETs is drain-induced barrier lowering [116]. In this process, a potential applied to the drain of the device affects the injected current at the source-gate barrier, thus leading to an increase in current with drain potential in a region where saturation is supposed to occur. It is not well appreciated that this effect will be dramatically increased when the transport becomes ballistic, that is when there is insufficient scattering in the channel of the FET.

Ballistic transport in semiconductors also is a relatively old idea. In the strictest sense, ballistic transport means the lack of any scattering. It is often discussed in mesoscopic structures, where the mean free path is comparable to the device size, in connection with the Landauer formula [117], but the ideas of ballistic transport are older, and derive from the earliest discussion of vacuum diodes. The Langmuir–Child law describes the ballistic transport of electrons in a thermionic diode, with space charge built up near the cathode (corresponding to our source in a MOSFET) [118, 119]. More recently, Shur and Eastman suggested that ballistic transport in ultra-short channel length semiconductor devices would have the same space charge and current relationship as the diode [120]. Drain-induced barrier lowering is one of the first indications of ballistic behavior [121], and leads to qualitatively similar curves.

However, one cannot argue just from the shape of the curves that ballistic transport is present, but must carry out the comparison with, and without, scattering as done in [122]. The impact of ballistic transport can be seen from a normal derivation of the current in a MOSFET, where

$$\begin{align}\int\limits_0^L I {\text{d}}y = W{C_{ox}}\int\limits_0^L v \left( {{V_G} - {V_T} - {V_y}} \right){\text{d}}y,\end{align} \tag{ 12 }$$

where C_ox is the gate oxide capacitance, L is the channel length, V_G and V_D are the gate and drain potentials, W is the gate width, v is the velocity along the channel, and V_y is the surface voltage along the channel. Normally, one would now introduce the mobility, but in the case of ballistic transport this has no meaning. Rather, the velocity is a function of the local potential as $v = \sqrt {2e{V_y}/{m^*}}$. In the ballistic case, the potential can be described through [119]

$$\begin{align}\frac{{{\text{d}}{V_y}}}{{{\text{d}}y}} = \alpha \sqrt I V_y^{1/4},\,{\alpha ^2} = 8\pi \sqrt {\frac{{2e}}{{{m^*}}}} .\end{align} \tag{ 13 }$$

Introducing this into (12), one finds that saturation occurs in the same manner, but that

$$\begin{align}I = K{\left( {{V_G} - {V_T}} \right)^{3/2}}\,\,,\,\,K = {\left[ {\frac{{4\alpha W{C_{ox}}}}{{90\pi L}}} \right]^{2/3}}.\end{align} \tag{ 14 }$$

Hence, the appearance of saturation has little to do with the properties of the transport in the channel, although the voltage dependence of the current does change somewhat [123]. But, this result arises entirely because it is assumed that pinchoff does occur in the channel. The Langmuir–Child law for diodes does not have any restrictions on the carrier motion, such as occur with pinchoff. If the pinchoff restriction is removed, then one might well expect diode/triode like behavior in the MOSFET, especially under quantum conditions.

In materials like InAs, the transport length for the transition from ballistic to resistive transport can be 15–20 nm at room temperature [124]. Yet, this is sufficiently long to affect the device characteristics as is evident in figure 14, which is a quantum simulation of a 30 nm gate length InAs nanowire MOSFET [125]. The cross-section of the channel is 9 × 8.5 nm, and the channel region is undoped. It is clear that this device has some strange behavior in the gate characteristics near turn-on, and does not saturate in the normal manner. Studies of the transport itself confirm that it is almost completely ballistic in nature. This can explain the diode/triode like behavior of the drain current.

Zoom In Zoom Out Reset image size

But ballistic transport is only one of the effects that arise in short channels. When the transport is fully quantum, then the transitions from the low-dimensional channel to the larger three-dimensional source and drain begin to play a role. When crossing one of these latter interfaces, there is a discontinuity in the carrier momentum, and since this occurs at both source and drain, the possibility arises for quasi-bound states, often called resonances, in the longitudinal wave function for carriers in the channel. That is, at certain biases, the drain current can have large peaks that are almost like resonant tunneling peaks. These may be seen in figure 15, which is a quantum simulation for ballistic transport in a 9.8 nm channel length Si MOSFET [126]. The channel cross-section is 18.5 nm wide and 6.5 nm in depth. The channel area is doped p-type with 5 × 1018 cm-3 concentration. In panel (a), the gate characteristics are shown for six different device doping configurations; the dopants are randomly placed according to the doping concentration and the grid size in the simulation. The peaks are the random fluctuations discussed above for atomistic doping when only a few dopants are present. There are, on average, only six dopants in the channel, so that their exact positions will introduce considerable variability. In panel (b), the drain characteristics are shown for a single device. While the small fluctuations are dopant caused, the large peaks are the result of longitudinal resonant levels arising from the mis-matches at the source and drain transitions. In a sense, these can be considered as source-drain resonant tunneling.

Zoom In Zoom Out Reset image size

The long road to FinFETs began many decades ago when it was realized that the danger of scaling would be seen in a situation in which the standby 'off' current of a chip was larger than the operating current. The road to FinFETs can largely be said to be a road in which control of the off current was the leading concept. The off current arises from the bipolar like nature of the n⁺ –p–n⁺ doping of a standard n-channel MOSFET (the complement is true for the p-channel). Much of the leakage current that provides the off current is a weak bipolar behavior in which the substrate plays an important part. The first step along the road was the concept of silicon-on-insulator (SOI) [128], as a means to suppress the substrate contribution to the off current. Indeed, the first chip appeared with this technology shortly after [129], and it was shown that these transistors did have a better sub-threshold slope (reduced leakage) [130].

The next step to control was the double-gate MOSFET, with top and bottom gates to better control the pinchoff [131]. But, this was a complicated and expensive fabrication process, so others began to look at different gate configurations that offered similar control [132, 133]. Thus was the arrival of the FinFET in which the channel is turned to the vertical direction with gates on either side. A variant was called the tri-gate transistor, shown in figure 16 [127], although gate materials have changed since this time [134].

Zoom In Zoom Out Reset image size

Turning the channel vertical allows the surrounding gates to work together to dramatically reduce the off current in these devices. Nevertheless, new roughness appears due to fluctuations in the body thickness, essentially, the fin thickness [135]. Despite many perceived problems and various approaches, it was shown that the trigate design was a more scalable transistor than the conventional planar MOSFET [136]. In 2009, integrated circuits using the FinFETs and trigate FETs began to appear [137], although there were still experiments using alternative materials to Si [138]. Nevertheless, Intel introduced the trigate as a mainstream device at the 22 nm node in 2011, and the industry has not looked back.

With the gate potential on either side of the fin, turnoff of the transistor becomes more effective as these two potentials work together to shut off the channel. The normal state in which there are two inversion layers, one on either side of the fin, can change into a single inversion layer in the bulk of the fin. When the fin is sufficiently thin and the density is not too high, this central inversion layer is the preferred state [42]. Thus, the quantum effects can change the basic nature of the FET, from a surface-oriented device to a bulk-oriented device, which has the added benefit of less scattering from the rough interfaces. In figure 17, the charge density is plotted as a function of the fin width, for a total inversion density of 3 × 10¹² cm⁻². For smaller fins, one gets bulk inversion as opposed to surface inversion. Of course, larger gate voltages also pull the charge to the surface, enhancing surface inversion.

Zoom In Zoom Out Reset image size

Even though one has better gate control in the FinFet, it should not be assumed that random doping effects will go away. This is just not in the cards [139], although here it seems that the devices are most sensitive to dopants in mid-channel. However, it does appear that the device is less sensitive to short-channel effects, and have low 1/f noise [140].

Once the concept of multiple gates was brought forward, then thoughts turned to what could be done with quantum wires, and a gate-all-around (GAA) technology. The first quantum wire device was simply a very narrow channel MOSFET, of 10 nm width [141]. The announced purpose of the device was to begin to study quantum effects in ultra-small MOSFETs, although the SOI device had a gate length of 250 nm, and the Si layer was 7 nm. The channel doping was only 1015 cm^–3, so that there was typically no dopants in the channel, although they studied a range of channel widths (1.25–43.75 nm). They observed an increase in threshold voltage as the channel width was reduced below ∼10 nm. Simulations of these devices soon appeared [142, 143].

The GAA Si nanowires seem to have appeared both theoretically [144] and experimentally [145] around 2004. The fabricated wires were about 12 nm high and 20 nm wide, nominal gate lengths of 100 nm, and were surrounded with oxide, although the gate was closer to a trigate than a GAA. These showed excellent sub-threshold slopes, comparable to a large planar device. In simulation, it was shown that these devices had the same bulk-to-surface channel effects, discussed above for the FinFETs, with varying diameter wires [144]. In addition, these simulations did not show any increase in mobility in these nanowires. More advanced forms of simulation continued to appear [146].

It was pointed out Moore's law results from an economic law rather than a technology law [16, 17]. The economics boils down to the cost of Si itself, and leads to the fact that, for Moore's Law to proceed, the gate periphery must be larger than the Si area in order to minimize Si cost, as the technology proceeds. This certainly supports the transition to FinFETs. But, it was shown that nanowires laid on the surface of Si could not satisfy this economic requirement, and could not compete with FinFETs [147]. This is because you basically cannot pack a single layer of circular wires dense enough to overcome the height advantage of the FinFET. An alternative would be to have the wires vertical [16], but that is a more difficult manufacturing technology. Of course, a better way was found, and this is stacked nanowires [148, 149]. The nanowire stack was good, but it was not that much preferred to the FinFET. This changed with the nano-sheet FET, which was a stack of squashed wires—the vertical thickness was much less than the width of each wire [150], although the phrase nano-sheet would come later [151], in a joint effort of IBM, Samsung, and Global Foundries.

By utilizing other materials, such as SiGe and Ge, a range of strain possibilities arise in these GAA stacked devices [152, 153]. Generally, the substrate under the stack can undergo punch-through by the source-drain bias which leads to leakage, so that some form of punch-through stopping layer is used [154]. Later, the use of a bottom oxide was used to provide this behavior and to give better control [155]. In figure 18, a nanosheet stack is shown [156]. Although this stack has four nano-sheets, the usual and preferred, is to use only three nano-sheets. Use of the oxide led to better sub-threshold slope and less drain-induced barrier lowering.

Zoom In Zoom Out Reset image size

This has led to the announcement this year of IBM's 2 nm node technology using their nano-sheets and oxide spacer layers [157]. In this approach, the nanosheets are 40 nm wide and 5 nm thick, with an effective gate length of 12 nm. It is expected that this technology will appear in mainline production in 2024. But, other technologies are being examined for the nanosheets, since these are formulated by various epitaxial and deposition methods. Indeed, nano-sheets have been studied in the III–Vs, as mentioned above [138]. More recently, the monolayer transition-metal di-chalcogenides (TMDC) have been studied in this association [158]. Materials such as MoS₂ and WS₂ offer the ultimate in nano-sheets, as the layer is a two-dimensional material less than a nm in thickness. These materials have a sizable band gap and reasonable velocities and mobilities [159, 160].

Perhaps as important is the problem of keeping the n-channel and p-channel transistors of the CMOS pair close together. One suggestion has led to a variety of approaches that lead to perhaps stacking the n and p devices on top of one another as part of the nanosheet stack [161, 162]. One approach to this is the forksheet technology which utilizes a self-aligned common gate between the two devices.

The Spin-FET was theoretically predicted over 30 years ago [163] (for relevant reviews see, e.g. [56, 164–166]). This particular type of transistor uses the spin properties of electrons for switching. In essence, a Spin-FET is based on two ferromagnetic contacts (source and drain) connected by a semiconductor channel. Spin-polarized electrons are injected via the source contact into the semiconductor region. The manifesting channel current is modulated by the gate-voltage-dependent spin–orbit interaction [167] which results in electron spin precession during transport. The drain contact's magnetization acts as a filter, as only spin-aligned electrons can pass through. All ballistic electrons have the same spin rotation at the end of the channel, which is linked to the spin–orbit field dependence on the electron momentum. As a consequence, the spin–orbit interaction strength dictates the minimum required length of the semiconductor channel for sufficient spin precession. In principal, two dominant spin-orbit interaction mechanisms are known: Rashba (geometrically induced structural asymmetry [167]) and Dresselhaus (bulk inversion symmetry breaking [168]). In silicon thin films spin–orbit interaction is Dresselhaus-dominated [169, 170], and is due to interfacial-induced breaking of the inversion symmetry (for transport modeling see, e.g. [171, 172]). The spin–orbit interaction strength depends almost linearly on the effective electric field and has been reported to be in the order of 2 μeV nm [173], making confined silicon structures candidates for Spin-FET channels ([100] oriented thin silicon film channels seem to be those best suited [174, 175]). However, among the challenges is the fact that a channel length in the order of micrometers is necessary. This requirement is contradictory to the ongoing geometric device down-scaling. Modern semiconductor devices offer roughly two orders of magnitude shorter channels: a severe competitive disadvantage of silicon Spin-FETs. Moving forward, required channel lengths can, in principal, be reduced by increasing the gate-voltage-dependent spin–orbit field strength (e.g. III–V materials). Another challenge is electron-phonon scattering which introduces randomization and acts against the spin precession coherence. Tackling this challenge requires cryogenic operational temperatures, further limiting potential applicability of Spin-FETs. The original limitation of requiring ferromagnetic contacts by injecting spin-polarized electrons through electrostatically created point contacts has been overcome [176]. More recent work on Spin-FETs has investigated alternative channel (e.g. 2D materials) and electrode materials (e.g. cobalt) as well as multi-gate and multi-functional logic devices and systems (for recent reviews see, e.g. [165, 166]).

In 2015, researchers showed the first demonstration of a Spin-MOSFET with a high on/off ratio operating at room-temperature [177]. Two metallic ferromagnetic contacts (source and drain) are connected by a non-magnetic semiconductor channel, allowing for charge and spin transport. Parallel/antiparallel magnetization alignment between source and drain leads to a current increase/decrease at the drain contact respectively. The ability to change the contact magnetization of the contacts by an external magnetic field and/or by the current (spin-transfer torque) provides opportunities for reprogrammable logic [178]. A particularly important feature of Spin-MOSFETs is the fact that the contact magnetization is preserved without external power, partly enabling non-volatile logic devices. However, in contrast to Spin-FETs, spin orientation is solely determined by the injecting ferromagnetic contact orientation.

As a concluding remark on the matter of Spin-FETS, let us point out that from an efficiency point of view, Spin-FETs only hold true advantages over conventional transistor designs when no current flow is required for the fundamental transistor switching mechanism. Indeed, only if this switching mechanism is realized solely via spin manipulation, will Spin-FETs be able to advance to a high-impact transistor technology. For reviews of recent applications of Spin-FETs see in particular [165, 166].

While the idea of tunneling in semiconductors is quite old, the idea of putting a tunnel junction into an FET seems to have appeared around 2007 [179]. At this time, the problem of poor sub-threshold slope and leakage current was becoming ever more important, and of course led to the rise of FinFETs and now nanowire FETs, as discussed previously. But, the idea of greatly improving the sub-threshold slope by using tunneling from the source into the channel was quite promising, and was quickly pursued by some [180, 181]. In this approach, a resonant tunneling structure is placed between the source and the channel, to create a large resistance in the sub-threshold region; but a smaller resistance through resonant tunneling as the device turns on. Interband tunneling from a p-source into an n-inversion channel did improve the sub-threshold slope. But there was a correlated problem with the device, and that was that the tunneling barrier lowered the available 'on' current in the device. This problem seems to be intrinsic to its very concept [182]. While there have been many approaches to try to raise the on current, including III–V materials [183] and monolayer compounds [184, 185], it does not appear that this problem has gone away. Almost immediately, it was suggested that graphene would be a suitable material for this device [184]. However, the device is still under study today, following ideas such as anisotropic insulators [186]. Yet introducing tunneling generally lowers device performance, in particular by limiting the available current that is intrinsic to a tunnel barrier, and this may be a problem not easily overcome.

The normal approach to solving transport using the Schrödinger equation is straight from the textbook. This involves matching the wave function and its derivative at each slice of the structure (across the interface, as shown in figure 21). But, this is generally unstable when a large number of such interfaces are present, due to the presence of evanescent modes that must be accurately treated. This leads to both increasing and decreasing exponentials, and that is the source of the instability. Stability can be restored by modifying the recursion to one based upon the scattering matrix, which has long been a staple of microwave systems and entered quantum mechanics through the Lippmann–Schwinger equation [205]. The stability of this approach lies in the fact that modal solutions maintain their orthogonality through the scattering process [205]. In small devices, the connection to microwave theory becomes stronger as the transport becomes dominated by the modes introduced by the lateral confinement, much as in a microwave waveguide. The use of scattering states provides a method of building up an orthogonal ensemble, even with weighting of the states by e.g. a Fermi–Dirac distribution [206]. This is particularly useful for the initial conditions of the wave function in a particular semiconductor device [207]. In the end, this approach allows to determine the transmission from one end to the other and evaluate the conductance via the Landauer formula.

Normally, the transport would be governed by the size of the Hamiltonian. In the scattering matrix approach however, the matrix is twice that of the Hamiltonian due to the existence of waves traveling in both directions. These are separated to avoid the exponential problem of wave function matching. Thus, the modes and their momentum wave numbers are computed at the far left slice in figure 21, via the matrix [208–210]

$$\begin{align}{T_0} = \left[ {\begin{array}{*{20}{c}} {{U_ + }}&{{U_ - }} \\ {{\lambda _ + }{U_ + }}&{{\lambda _ - }{U_ - }} \end{array}} \right].\end{align} \tag{ 18 }$$

This is a scattering matrix form. Prior to this, the Hamiltonian for this first slice is solved to find the eigen-functions and the eigen-energies. The U matrices are then the similarity matrices for the forward and backward solutions for the wave functions, while the λ matrices are the eigen-value matrices. Normally, we think of the similarity matrices being the transformations on the Hamiltonian as

$$\begin{align}{U^\dagger }HU = EI,\end{align} \tag{ 19 }$$

where I is a unit matrix. In other words, this matrix allows one to diagonalize the Hamiltonian, and find the eigen-values for the particular slice. Since, these eigen-values, for propagating waves, can correspond to forward or backward propagation for that particular slice, the system is doubled in dimension in order to separate the two types of eigen-values [208]. When the Hamiltonian is inverted to find the eigenvalues, a result is the similarity transformation matrix U, which is used to diagonalize it. This matrix is also computed during the process, and it (or its adjoint, depending upon the details of how one defines it—the order of matrices in (19)) contains valuable information. The columns of this matrix are the modes and their values at each site on the grid used to solve Schrödinger's equation. The rows are the values of the various modes at a particular grid point. This matrix is now a mode-to-site transformation matrix.

The matrix (18) changes the mode wave functions into the site wave functions in this zero-slice initialization of the process. Process through the structure is generated by solving the required normalization to a unit matrix of the following set of matrices [55, 57, 211]:

$$\begin{align} \left[ {\begin{array}{*{20}{c}} {C_{j + 1}^{\left( 1 \right)}} & {C_{j + 1}^{\left( 2 \right)}} \\ 0&I \end{array}} \right]& = \left[ {\begin{array}{*{20}{c}} 0&I \\ I& {T_l^\dagger \left( {{H_j} - EI} \right)} \end{array}} \right]\nonumber\\ & \quad \times \left[ {\begin{array}{*{20}{c}} {C_j^{\left( 1 \right)}}& {C_j^{\left( 2 \right)}} \\ 0&I \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 0&0 \\ {P_j^{\left( 1 \right)}}& {P_j^{\left( 2 \right)}} \end{array}} \right],\end{align} \tag{ 20 }$$

where the I are unit matrices. $C_j^{\left( 1 \right)}$ and $C_j^{\left( 2 \right)}$ are the amplitudes of the forward and backward modes respectively. The initial values of $C_0^{\left( 1 \right)}$ are determined as the value of the Fermi function in the far left slice (between 0 and 1), while $C_0^{\left( 2 \right)} = 0$. The P matrices will be found from the basic scattering matrix requirement that

$$\begin{align}\left[ {\begin{array}{*{20}{c}} {{\psi _j}} \\ {{\psi _{j + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0&I \\ {T_l^\dagger {T_l}}&{T_l^\dagger \left( {{H_j} - EI} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\psi _{j - 1}}} \\ {{\psi _j}} \end{array}} \right].\end{align} \tag{ 21 }$$

Equating these last two equations gives

$$\begin{align}P_j^{\left( 1 \right)}& = C_{j + 1}^{\left( 1 \right)} = - P_j^{\left( 2 \right)}C_j^{\left( 1 \right)}\nonumber\\ P_j^{\left( 2 \right)}& = C_{j + 1}^{\left( 2 \right)} = {{\left[ {C_j^{\left( 2 \right)} + T_l^\dagger \left( {{H_j} - EI} \right)} \right]}^{ - 1}}.\end{align} \tag{ 22 }$$

These are now propagated to the Nth slice, which is the end of the device and H_r region, and then onto a terminating slice. For various points in the device, T_l will become T_dl, as indicated in figure 21. At this point, the inverse of the mode-to-site transformation matrix is applied to bring the solution back to the mode representation, so that the transmission coefficients of each mode can be computed. The density at each point in the device is determined from the wave function squared magnitude at that point, which is back propagated from the terminal so that the wave function at site i in slice j is given as

$$\begin{align}{\psi _{i,\,j}} = P_j^{\left( 1 \right)} + P_j^{\left( 2 \right)}{\psi _{i,j + 1}}.\end{align} \tag{ 23 }$$

This solution technique for the direct Schrödinger equation was used for the MOSFET depicted in figures 14 and 15 above.

If there is no scattering, and the device transverse dimensions remain constant from one end to the other, then the various modes will be uncoupled, so one can simply use the Landauer formula. This provides a type of ballistic transport for a small device and has been used to study wave function penetration into the gate oxide [212], and to study a double-gate FET in TMDCs as a biological sensor [213]. A slightly different formulation, termed the quantum transmission boundary method [214], has been used to study the effect of defects in graphene nanoribbon FETs [215].

It is not necessary to avoid including scattering, as it can naturally be added via a self-energy correction [216] (the self-energy is discussed further below). Primarily, the scattering is represented by the imaginary part of the self-energy, and this was calculated for all common scattering processes [216]. With these scattering processes, the crossover from ballistic to collision-dominated transport (usually referred to as the ballistic-to-diffusive crossover) was studied for Si and for a variety of temperatures and materials [124, 216, 217]. Thus a direct solution method provides a viable simulation technique for FETs that can be done on a standard desktop computer and provides accurate representations of the quantum physics.

The second approach, which is one of the most straight-forward approaches, creates a matrix of modes. Generally, the Schrödinger equation is solved by assuming an expansion of the wave function in a suitable basis set so that each basis function is an energy eigen-function according to

$$\begin{align}H{\varphi _n} = {E_n}{\varphi _n},\end{align} \tag{ 24 }$$

just as is done for the regular solutions of the previous section. Here, E_n is the energy level corresponding to the particular basis function ${\varphi _n}$. For example, in a quantum well, the basis functions are the set of wave functions corresponding to the bound states of the quantum well. In a quantum wire, these are the bound states arising from the quantization resulting from the wire width (and height). In both of these cases, the basis functions are the mode wave functions. They have values on the sites of the grids used in the computation. Then, the total wave function can be written as

$$\begin{align}\psi \left( {r,t} \right) = \sum\limits_n {{c_n}} {\hat \varphi _n}\left( r \right){{\text{e}}^{ - i{E_n}t/\hbar }}.\end{align} \tag{ 25 }$$

For the general situation, one may then define the density matrix, as a single time function, to be

$$\begin{align} \rho \left( {r,r^{\prime} ,t} \right) = \sum\limits_{m,n} {{c_{mn}}} \hat \varphi _m^\dagger \left( r \right){{\hat \varphi }_n}\left( {r^{\prime} } \right) \times {{\text{e}}^{i\left( {{E_m} - {E_n}} \right)t/\hbar }}. \end{align} \tag{ 26 }$$

The off-diagonal elements are rapidly oscillating in the difference frequency, so that they are usually ignored in transport theory. However, these off-diagonal elements are generated in the complex device and represent entanglement of the various modes, such as that occurring in figure 12 [218].

Generally, the device is embedded within an environment, such as shown in figure 21, and this environment can include other devices. So, the total Hamiltonian can be written as [31]

$$\begin{align}H = {H_{en}} + {H_{dev}} + {H_{e - d}},\end{align} \tag{ 27 }$$

where the last term represents interactions between the device and its environment. Since measurements are always made in the environment [59], the coupling between the device and the environment is crucially important. The reduction to the dashed box in figure 21 can be achieved through the use of projection operators [219], in which the goal is to find the equation of motion for the reduced density matrix

$$\begin{align}{\rho _{dev}} = T{r_{en}}\left\{ \rho \right\}.\end{align} \tag{ 28 }$$

The problem is that the equation of motion arises from the normal Liouville equation in which

$$\begin{align}i\hbar \frac{{\partial \rho }}{{\partial t}} = H\rho - \rho H = \left[ {H,\rho } \right],\end{align} \tag{ 29 }$$

where the last term defines the commutator relationship. It is to this equation that the projection operators are applied [31, 53, 57]. The details will not be presented here, as they are complicated and the result agrees with intuition. The end result is that the equation of motion for the device density matrix is found to be

$$\begin{align}i\hbar \frac{{\partial {\rho _{dev}}}}{{\partial t}}& = \left[ {\left( {{H_{dev}} + \overline {{H_{e - d}}} } \right),{\rho _{dev}}} \right]\nonumber\\ & \quad + \mathop \int \limits_0^t {\text{d}}t^{^{\prime}} \left[ {{{\Sigma }}\left( {t^{^{\prime}}} \right),{\rho _{dev}}\left( {t - t^{^{\prime}} } \right)} \right].\end{align} \tag{ 30 }$$

where in the last term, ${{\Sigma }}\left( {t^{^{\prime}} } \right)$ is the self-energy mentioned above and represents scattering within the device and between the device and its environment. The reduced interaction (represented by the overhead bar) in the first term represents effects of the environment on the device, such as boundary conditions affected by the environment.

There are two common approaches followed from here to obtain solution techniques of use for the solving of the properties of a device. In the first approach, a set of balanced equations for the density, momentum (current), and energy are obtained [220, 221], just as are normally done with the Boltzmann equation [222]. These equations are referred to as the (quantum) hydrodynamic equations. These equations were used to study tunneling in single and double-barrier structures [64, 223]. Often, these approaches bring the explicit quantization in a confined structure into the problem via an effective potential, and this approach was used to study the carrier density distributions in a GaAs HEMT [224].

It is possible to directly use the Liouville equation (28), and try to solve this equation directly. However, it is not obvious that one can do this for an application in a small semiconductor device. One problem is that we cannot just do the diagonal terms in the density matrix, but must also properly account for the off-diagonal terms [225]. But there is a long history of this approach [226]. The presence of the reservoirs, or contacts, generally lead to interferences that appear within the off-diagonal terms [227–230]. Indeed, this approach has been used to study the transport through a small ${n^ + } - n - {n^ + }$ structure in Si [231, 232]. More recently, this approach has been applied to a double-gate ultra-thin-body Si MOSFET [233]. The channel region was 5 nm thick and 10 nm long, and lightly doped. The source and drain regions were more heavily doped. In figure 22, the potential throughout this device is shown. The effect of phonon scattering and interface roughness scattering are shown in the distribution functions in figure 23, where the diagonal terms of the density matrix are plotted as a function of energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second approach develops a proper quantum kinetic equation and then solves it directly for the density matrix, an approach that is amenable to using ensemble Monte Carlo particle methods [234]. The density matrix is more complicated than a classical approach as there are two positions and the time along either of these must be considered. But a path integral equation can be constructed from the equation of motion for the density matrix, and this serves as the basis for developing a Monte Carlo procedure. More importantly, this procedure is not limited to just the lowest-order scattering processes [235].

General Monte Carlo techniques coupled with the density matrix have been discussed by Jacoboni [236]. This approach was used to follow time-dependent wave packet trajectories, and the Bohm potential, in order to build up the density matrix [237]. Monte Carlo was also used with the density matrix to study interband tunneling of holes at high electric fields in GaN [238].

With the classical Boltzmann equation, one has a nearly equilibrium situation in which the distribution function is a Maxwellian or a Fermi–Dirac at a temperature possibly above that of the lattice. One often couples this with an assumption that the contacts are in near-equilibrium, but this is not always true, as was seen in figure 13, and this assumption is no more likely to be true in the quantum case. Then in the hot carrier system, and that means most devices, the distribution is unknown and has to to be found via one of several possible methods to solve this version of the Boltzmann equation, as the currents (and therefore the device characteristics) are found from integrals over this distribution.

In quantum mechanics, the system is quite the same, giving rise to excitations that propagate under the action of the potential. The NEGF include new functions that must be found to describe this very non-equilibrium distribution [239, 240]. With these Green's functions, the wave functions $\psi$ and ${\psi ^\dagger }$ go beyond the simple forms (25), as the coefficients become expressed as fermion creation and annihilation operators. Thus, the wave functions are referred to as field operators [57], and the Green's functions determine the transport of an excitation over a certain distance (given by the two positions and times in the arguments) to where it is annihilated (or destroyed). In this sense the Green's functions describe transport through the system. In terms of these two field operators, one uses the anti-commutators for fermions as

$$\begin{align}\left\{ {A,B} \right\} = AB + BA,\end{align} \tag{ 31 }$$

and the four Green's functions in NEGF are described by

$$\begin{align}{G_r}\left( {r,r^{^{\prime}};t,t^{^{\prime}} } \right)& = - i{{\Theta }}\left( {t - t^{^{\prime}} } \right)\left\{ {\hat \psi \left( {r,t} \right),{{\hat \psi }^\dagger }\left( {r^{^{\prime}},t^{^{\prime}} } \right)} \right\} \nonumber\\ {G_a}\left( {r,r^{^{\prime}};t,t^{^{\prime}} } \right)& = i{{\Theta }}\left( {{t_0} - t} \right)\left\{ {{{\hat \psi }^\dagger }\left( {r^{^{\prime}},t^{^{\prime}} } \right),\hat \psi \left( {r,t} \right)} \right\} \nonumber\\ {G^ < }\left( {r,r^{^{\prime}};t,t^{^{\prime}} } \right)& = i{{\hat \psi }^\dagger }\left( {r^{^{\prime}},t^{^{\prime}} } \right)\hat \psi \left( {r,t} \right)\nonumber\\ {G^ > }\left( {r,r^{^{\prime}};t,t^{^{\prime}} } \right)& = - i\hat \psi \left( {r,t} \right){{\hat \psi }^\dagger }\left( {r^{^{\prime}},t^{^{\prime}} } \right).\end{align} \tag{ 32 }$$

Here, ${G_r}$ and ${G_a}$ are the retarded and advanced Green's functions respectively, and derive from the near equilibrium case where the distribution function is a Fermi–Dirac function. These define the spectral density as

$$\begin{align}A\left( {k,\omega } \right) = - Im\left\{ {G_r^{}\left( {k,\omega } \right) - G_a^{}\left( {k,\omega } \right)} \right\},\end{align} \tag{ 33 }$$

where k and ω are Fourier transforms of $r - r^{\prime}$ and $t - t^{^{\prime}}$ respectively. The spectral density relates the energy to the momentum of the quantum state. Classically, this is a delta function asserting that, for example $E = {p^2}/2{m^*}$, where the momentum $p = \hbar k$ and ${m^*}$ is the effective mass in the semiconductor. In quantum mechanics, the presence of state broadening leads to a breakdown of this equality, thus leading to the need for the spectral density. These two Green's functions, ${G_r}$ and ${G_a}$, are operator expressions, hence the anti-commutation (curly) brackets common for fermions.

The other two Green's functions in (32) are not operators, but are correlation functions related to the distribution function that must be found for the non-equilibrium device. The connection between these two functions is not accurately known in the far from equilibrium case, but has been hypothesized to maintain the same form as found near equilibrium [240]. This ansatz thus gives the less-than function as [241]

$$\begin{align}{G^ < }\left( {k,\omega } \right) = f\left( \omega \right)A\left( {k,\omega } \right),\end{align} \tag{ 34 }$$

and the greater-than function as

$$\begin{align}{G^ > }\left( {k,\omega } \right) = \left[ \,{f\left( \omega \right) + 1} \right]A\left( {k,\omega } \right),\end{align} \tag{ 35 }$$

where $f\left( \omega \right)$ is the nonequilibrium distribution function. As one can infer from the multitude of functions, using the Green's functions is far more difficult than using either wave functions or the density matrix.

The correlation functions (34) and (35) require a pair of equations of motion for these quantities. With the added complications of having four functions, these equations will be more complicated than the simple Boltzmann equation. Keldysh [242] introduced a general method with a single matrix Green's function, and hence a single matrix equation to be solved. This does not reduce the overall effort required, but it simplifies the equations in which these functions are described. The Keldysh matrix involves a contour that defines the evolution of time between $t$ and $t^{^{\prime}}$, either forward or reverse in time. This contour is drawn from the initial time (t₀ or 0) to the time of interest (the lower part of the overall contour) and then back to the initial time (the upper part of the contour). Each of these field operators can be on the upper line or the lower line of the overall matrix. The matrix is a 2 × 2 matrix. The wave functions enter the Green's functions through (32). The rows of the Keldysh matrix are defined by the operator $\psi \left( {r,t} \right)$. That is, row one of the matrix corresponds to when t is on the lower line, while row two is defined when t is on the upper line of the trajectory. Similarly, the operator ${\psi ^\dagger }\left( {r^{^{\prime}},t^{^{\prime}} } \right)$ defines the columns of the matrix. Column one is defined when $t^{^{\prime}}$ is on the lower part of the trajectory, while column two is defined when $t^{^{\prime}}$ is on the upper part of the trajectory. This is discussed fully in more detailed treatments [57], but is not all that germane to the discussion here. This contour ordered matrix may be written as:

$$\begin{align}{G_C} = \left[ {\begin{array}{*{20}{c}} {{G_r}}&{{G_K}} \\ 0&{{G_a}} \end{array}} \right],\end{align} \tag{ 36 }$$

where the Keldysh function is

$$\begin{align}{G_K} = {G^ < } + {G^ > },\end{align} \tag{ 37 }$$

and ${G^ < }$ and ${G^ > }\,$are given by the third and fourth lines of (32), respectively. From the Liouville equation (28), the equation of motion for the Green's functions can be written as [242]

$$\begin{align}& \left[ {i\hbar \frac{\partial }{{\partial t}} - {H_0}\left( r \right) - V\left( r \right)} \right]{G_K} = \hbar I - {{\Sigma }}{G_K} \nonumber\\ & \left[ { - i\hbar \frac{\partial }{{\partial t^{^{\prime}} }} - {H_0}\left( {r^{\prime} } \right) - V\left( {r^{\prime} } \right)} \right]{G_K} = \hbar I - {G_K}{{\Sigma }}.\end{align} \tag{ 38 }$$

Here, the bold-face characters are the matrix forms following (36). The Keldysh approach and its failures will be discussed further in the appendix.

Most modern approaches to the use of NEGF for devices base their work on the Keldysh [242] formulation. The scattering self-energies ${{\Sigma }}$ depend on the correlation functions ${G^ < }$ and ${G^ > }$ and on the greater and lesser Green's functions of the phonons, which account for the occupancy of the phonon states and depend on the phonon energies. Keldysh himself uses the S-matrix from equilibrium theory that works only upon the assumption that the system is close to equilibrium, so that the normal techniques can be used. Keldysh also assumes the interaction representation, in which the scattering-derived ${{\Sigma }}$ is assumed to follow from a unitary time-translation operator. But this assumption fails when the energies contain self-energies so that the Hamiltonian is no longer Hermitian. More importantly, Danielewicz has pointed out that correlated initial conditions also make the S-matrix questionable at best [243].

The entire concept of the contour mentioned above is also debatable. When electric fields or forces are applied to the device, the entire complex system undergoes a phase transition that breaks time-reversal symmetry. The device, within its environment, then seeks a steady-state if it exists; a far-from-equilibrium stable state that balances the driving forces and the dissipative forces [53]. During the transient evolution to this stable state, the system may evolve though a number of intermediate phases: e.g. homogeneous, inhomogeneous, linear, nonlinear, etc. When these driving forces are removed, the system response does not reverse its course through these different phases, but seeks a relaxation toward the equilibrium steady-state from which it initially deviated. This excitation/relaxation cycle may involve significant hysteresis. The excitation process generates entropy [220]. The relaxation process does not remove this entropy from the system, but generates even more as dissipation still occurs during this relaxation. The adoption of a time-ordered contour that smoothly retraces itself to its initial state is contrary to this physics (discussed further in the appendix).

Many approaches to using the NEGF use this perturbation theory expansion to arrive at Σ, even in the case of impurity scattering. It has already been pointed out that this approach will not yield the interferences seen in figure 12. More properly, the terms ${{\Sigma }}{G_K}$ and ${G_K}{{\Sigma }}$ should be two-particle Green's functions, and their partition into a pair of single particle Green's functions is not clearly evident, as it ignores all correlations. The evaluation of a two-particle Green's function is usually much more complicated and cannot be reduced to a simple Dyson's equation. Rather, it needs the full evaluation of the Bethe-Salpeter polarization [57]

$$\begin{align} \Pi \left( {k,\omega } \right) = {G_r}\left( {k,\omega } \right){G_a}\left( {k,\omega } \right)& \Bigg\{ {{k^{\prime}} } \cdot k\delta \left( {{k^{\prime}} - k} \right) + \int {\frac{{{d^3}{k^{\prime}} }}{{{{\left( {2\pi } \right)}^3}}}}\nonumber\\ & \;\;\; \times \frac{{{k^{\prime}} \cdot k}}{{{k^2}}}\Lambda \left( {{k^{\prime}} - k} \right)\Pi \left( {{k^{\prime}} ,\omega^{\prime} } \right) \Bigg\}, \end{align} \tag{ 39 }$$

which then enters the conductivity as

$$\begin{align}\sigma \left( {k,\omega } \right) = - \frac{{{e^2}\hbar }}{{{m^{*2}}}}\int {\frac{{{d^3}{k^{\prime}} }}{{{{\left( {2\pi } \right)}^3}}}} \int {\frac{{d\omega^{\prime} }}{{2\pi }}} \Pi \left( {{k^{\prime}} ,\omega^{\prime} } \right)\frac{{\partial f\left( {\omega^{\prime} } \right)}}{{\partial \omega^{\prime} }}.\end{align} \tag{ 40 }$$

The distribution function $f\left( \omega \right)$ comes from (34). Failure to incorporate the Bethe–Salpeter equation should cause a questioning of results found with the use of NEGF.

Most approaches for actually solving for device response with NEGF rely on a version of the CBR mentioned above. With this approach, the Green's functions can be reduced to those of the device itself (the region in the dashed box in figure 21). The CBR leads to equations of the form [244]

$$\begin{align}{{{{G}}}^R}\left( E \right)& = {{{{A}}}^{ - 1}}\left( E \right){{{{G}}}_0}\left( E \right)\nonumber\\ & = \left[ {\begin{array}{*{20}{c}} {{{{A}}}_l^{ - 1}{{{G}}}_l^0}& {{{{A}}}_l^{ - 1}{{{G}}}_l^0} \\ { - {{{{A}}}_{dl}}{{{A}}}_l^{ - 1}{{{G}}}_l^0 + {{{G}}}_{dl}^0} & { - {{{{A}}}_{dl}}{{{A}}}_l^{ - 1}{{{G}}}_{dr}^0 + {{{G}}}_d^0} \end{array}} \right],\end{align} \tag{ 41 }$$

with

$$\begin{align}A\left( E \right) = I - {G_0}\left( E \right){{\Sigma }}\left( E \right).\end{align} \tag{ 42 }$$

The subscripts are written according to figure 21. The matrix A can easily be evaluated using the properties of the self-energy, which is presumed to be non-zero only where the system makes contact with the external leads. This obviously means that ballistic transport is being assumed throughout the device. But in the case of ballistic transport, NEGF is an unnecessary afterthought. For ballistic transport, the coupling to leads gives the self-energies that define the transmission from one lead through the device to the next lead as [245]

$$\begin{align}T = {{{\Gamma }}_l}{G_r}{{{\Gamma }}_r}{G_a},\end{align} \tag{ 43 }$$

and the trace over the matrix T is used in the Landauer equation (16). The Γ are the imaginarly parts of the self-energies that arise from T_dl and T_dr . Since the two Green's functions in this equation are the equilibrium Green's functions, use of NEGF seems to be overkill and is used more for hype than for need.

It is clear that the interference from impurities can be treated with NEGF from the earliest applications to real devices [103]. More recently, the role of fixed charges in Si nanowire FETs was examined, taking full account of the atomistic positions of the fixed charge within Poisson's equation (and not with perturbation theory) [246]. These authors clearly demonstrated how a few random dopants can greatly affect the nanowire FETs. The effect of random dopants in an InAs drain region for a Si tunnel FET was also examined [247]. Then, using DFT to calculate both the band structure and the full asymmetric wave functions of donor electrons in Si nanowire FETs, it was clearly established how this anisotropy could affect the transport properties [104].

NEGF was used early to study the properties of double-gate MOSFETs and showed that channel length modulation and drain-induced barrier lowering became significant for channel lengths below 15 nm [248]. It was also used to study a double-gate SOI Si MOSFET [249], and to study the role of phonon scattering in multi-material nanowire FETs and FinFETs [250]. In figure 24, the current spectra found in this latter work, overlaid on the potentials, are shown for two Si nanowire FETs of 2.2 × 2.2 N m² and 3.6 × 3.6 N m². It may be seen that the larger cross-section device has less change in the current. The hot electrons have not relaxed in this larger device even after 65 nm of drain region. It may be concluded from this that the denser density of states in the smaller wire has more dissipative scattering than the slightly larger wire, and that relaxation of the carriers extends well into the drain regions for quite a distance.

Zoom In Zoom Out Reset image size

Others have also studied nanowire FETs of various types [251] and FinFETs [252], as well as various tunneling structures [253]. Strong simulation packages using full band structure and NEGF for transport have appeared for FETs over the years [254, 255]. There have been several studies of alternate material FETs, typically graphene and TMDCs, that have used NEGF as well.

In classical physics one uses the Boltzmann equation, which yields a distribution function that is defined in position, momentum, and time. Hence, there has always been a desire to have a phase space formulation in quantum mechanics, and one was developed quite early by Wigner [68]. To approach this, and to clarify the two positions and time that appear in the density matrix, as well as the NEGF, the center-of-mass and difference coordinates are defined through

$$\begin{align}\begin{array}{*{20}{c}} {R = \frac{1}{2}\left( {r + r^{\prime} } \right)} \\ {s = r - r^{\prime} } \end{array}.\end{align} \tag{ 44 }$$

Like the density matrix, the Wigner function has only a single time, so that the times in NEGF are set equal as $t = t^{^{\prime}} \equiv t$. Then, the density matrix may be transformed on the difference coordinate as

$$\begin{align}{f_W}\left( {R,t} \right) = \frac{1}{{{h^3}}}\int {\text{d}} s\rho \left( {R + \frac{s}{2},R - \frac{s}{2},t} \right){{\text{e}}^{ip \cdot s/\hbar }}.\end{align} \tag{ 45 }$$

Using the Liouville equation (28), this leads to an equation of motion for the Wigner function that is expressed as (in one dimension)

$$\begin{align} \frac{{\partial {f_W}}}{{\partial t}} + \frac{p}{{{m^*}}}\frac{{\partial {f_W}}}{{\partial R}} - \frac{1}{h}\int {\text{d}} p^{\prime} W\left( {R,p^{\prime} } \right){f_W}\left( {p + p^{\prime} } \right) = 0 , \end{align} \tag{ 46 }$$

where

$$\begin{align} W\left( {R,p^{\prime} } \right) = \int {\text{d}} s\;{\mkern 1mu} {\text{sin}}\left( {\frac{{sp^{\prime} }}{\hbar }} \right) \times \left[ {V\left( {R + \frac{s}{2}} \right) - V\left( {R - \frac{s}{2}} \right)} \right]\!. \end{align} \tag{ 47 }$$

The term $W\left( {R,p^{^{\prime}}} \right)$ is often called the Wigner potential. Scattering has not been included, but it is easily incorporated by a gain-loss term added to the right-hand side of the equation. This is usual for treatment of devices, and we will return to this below.

The Wigner function is not a positive definite function under far-from-equilibrium conditions [256]. This is a consequence of the uncertainty relationship. If the Wigner function is integrated over position or momentum, a positive definite function gives the probabilities (square magnitude of the wave functions) in momentum or position respectively. The negative excursions of the function exist over phase space regions whose volume is of order ${\hbar ^3}$ (in three dimensions), so that smoothing the function over a volume corresponding to the uncertainty principle will produce a positive definite function. The Wigner function is always positive definite in the ground state appropriate to equilibrium. But the onset of negative excursions of the Wigner function are viewed as the appearance of correlation and/or entanglement [58]. Such correlation or entanglement is connected to off-diagonal elements of the density matrix [218], and these lead to oscillatory behavior in the Wigner function and the negative excursions. Because of the phase-space nature of the Wigner distribution, it is possible to identify where quantum corrections enter a problem by comparing it with the classical version. Wigner himself derived an effective quantum potential that can be used as a correction term for the potentials in thermodynamics [68]. It is worth remarking that many start their study of devices with NEGF, but finally resort to using the Wigner function [257].

As with Green's functions, early use of Wigner functions for transport considered resonant tunneling diodes [207, 258]. Such tunneling can be seen in figure 25, where a Gaussian wave packet is in the process of tunneling through a 3 nm AlGaAs barrier, with GaAs regions on either side [259]. The total wave packet still has components that represent the initial wave packet (on the right below the barrier, represented by the two black lines), the reflected wave (on the left below the barrier), and the transmitted wave (above the barrier). It is worth noting that a wave packet is composed of a great many waves, and the figure shows how the reflected components actually slow down and then accelerate into the backward direction. Reflection and transmission is not instantaneous in quantum mechanics, but occurs over a period of time, which is characteristic of the motion of the wave packet. The wave here is incident on the barrier at a slight angle, which is observed in the 'tilts' of the partial waves, and there is a small electric field applied which is apparent in the increasing momentum as the wave moves forward (upward).

Zoom In Zoom Out Reset image size

The Wigner equation of motion (45) is quite similar to the Boltzmann equation, except for the complication of the Wigner potential, which arises from the nonlocal properties of quantum mechanics. Nevertheless, it is possible to transform the equation of motion into a path integral form, from which a Monte Carlo procedure can be developed. As has been mentioned in the Introduction, the idea of using particles in quantum mechanics was discussed quite early by Kennard [2]. The concept of a Wigner path is thus useful, as it provides a representation of the quantum evolution, and leads to numerical simulation of the Wigner equation of motion. Such a path is followed by a small sample of the Wigner function as it evolves through phase space (these are represented by the particles used in the simulation). Moreover, this path evolves quite like a classical particle, except for scattering and for rapidly varying quantum potentials, as described by Kennard [2] and later by Bohm [51]. Using paths guided by the Bohm potential has been used extensively by Oriols [260]. In using particle approaches with the Wigner function, there are a few quantum issues that have to be addressed. First is the problem of including interference effects leading to negative excursions for the Wigner function. Second will be the nonlocal effects that have to be treated in scattering processes, especially those that appear in treating FETs.

The problem of phase is dealt with first. An early approach used a weighted Monte Carlo in which each particle was assigned a particular weight [261]. A variant to include paths with negative weights handles the negative parts of the Wigner function [262]. Because some of the Monte Carlo particles carry negative weights, a problem arises as to how to keep the density at its nominal value, as there may well be more simulation particles than real particles. This can be handled easily within the simulation, and the weights are termed affinities. The Wigner function at a point in phase space can be described by a sum over the particles as

$$\begin{align}{f_W}\left( {r,k} \right) = \mathop \sum \limits_{i = 1}^N {A_i}\delta \left( {r - {r_i}} \right)\delta \left( {k - {k_i}} \right),\end{align} \tag{ 48 }$$

where N is the number of simulation particles and ${A_i}$ is the weight, or affinity, attached to each of these particles. Absorbing boundary conditions on the exit boundaries from the device may be instituted in a very easy manner. An important aspect of the affinity is its sign and magnitude change during the simulation. The nonlocal Wigner potential determines the change in the affinity that each particle possesses. That is, the wave nature of the particles is maintained through the variation of each particle's affinity. This affinity was used for figure 3.

The weight and the affinity are really artificial numerical quantities whose purpose is to simulate the quantum phase interference that occurs during real propagation within a semiconductor device. A melding of these two concepts was pursued in another approach [263], which led to the adoption of a pure sign convention [264, 265]. This creates a model in which the Wigner function is considered to include a Boltzmann-like scattering term, and a generation term. The quantum information is carried by the sign of the quasi-particles. This approach is most useful when treating the interaction with the non-local potential as a scattering event. When a scattering event from the potential occurs, two new particles are created, and thus two new wave packets, one with the momentum increased by q and one where it is reduced by q, with q determined randomly from the probability distribution of the potential's spatial Fourier transform. The sign on one of these new particles is taken to be the same as the incident particle, while the sign on the other is the opposite. These signs are taken into account in each averaging process that is used to find average values. Equivalent particles with opposite signs annihilate one another when they meet in phase space, i.e. when their wave packets overlap. This approach was used for figure 12. The Wigner Monte Carlo approach was used to simulate double-gate MOSFETs [266], especially on the nano-scale [267]. In figure 26, the schematic of a nano-scale double-gate MOSFET is shown for these simulations. The source and drain regions are heavily doped and are degenerate, while the channel region is undoped and has a gate length of 6 nm. The Si layer thickness is taken to be 3 nm. In figure 27, the quantum calculation is compared to a classical calculation for the lowest sub-band occupancy in the sub-threshold region. Interference fringes may be seen in the quantum calculation, but not in the classical simulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above simulation for the sub-threshold behavior gives results similar to those for the transistor when it is turned on. In figure 28, the Wigner function is shown for the lowest propagating sub-band in a gate-all-around cylindrical nanowire Si MOSFET with diameter of 50 nm diameter and channel length of 60 nm [268]. It may be seen that the lowest sub-band, which is pictured, moves smoothly through the device, accelerating in momentum as it follows the electric field. Both ballistic transport and scattering due to impurities and acoustic phonons are shown. In either case, as well as the sub-threshold behavior of figure 27, the current-carrying carriers are not relaxed well in the drain. One can clearly see the separation of these carriers from the main body of carriers in the drain in figure 28. Similar behavior has been reported by Croitoru et al [269, 270].

Zoom In Zoom Out Reset image size

As remarked above, many approximations are used to include collision effects and scattering within the Wigner approach. These are typically the normal semi-classical scattering functions, modified slightly to account for off-shell effects and collisional broadening. Levinson has developed a formal approach to Wigner scattering [271], and his results support this approach. In quantum transport, scattering from a phonon, for example, does not have to be local; the event may spread over a small spatial distance. The collision is not instantaneous, but takes a few femtoseconds to be completed [272]. With the intra-collisional field effect, in which the particle continues to be accelerated by the field during this short time, the particle also moves a short distance, which makes it nonlocal. In addition, the loss of a one-to-one correspondence between energy and momentum leads to broadening of the energy conservation that arises from replacement of the delta function by the spectral density function.

While it is not something that is often discussed with the Wigner function, the momentum arises from the difference in the two positions of the wave functions, as may be seen in (44) and (45). The momentum is therefore normal to the main diagonal in the density matrix. The time duration is the time required to transit between the two positions at the corresponding momentum. These different variables are thus tied together in a self-consistent manner. Many methods have been suggested to include the collisional displacement as a modification to the use of semi-classical scattering within the Monte Carlo method [273–275]. To properly include the role of the finite collision duration, we need to account for two effects: (a) the collision duration itself and the resulting energy shift, and (b) the position change during the collision. During the collision, the position changes approximately as [276, 277]

$$\begin{align}\delta x = \left[ {v\left( t \right) + \frac{{eF \cdot v\left( t \right)}}{{{m^*}v\left( t \right)}}\frac{{{\tau _c}}}{2}} \right]{\tau _c},\end{align} \tag{ 49 }$$

where F is the electric field and ${\tau _c}$ is the collision duration. The first term is just simple displacement due to the particle's motion. The second term represents the intra-collisonal field effect. During the interaction, there is an overlap between the two wave packets at the two positions separated by $\delta x$. The two wave functions are Wigner packets which may be taken to be Gaussian in shape, which leads to the probability for a given shift

$$\begin{align}P\left( {\delta x} \right) = \frac{1}{{4\pi {\sigma ^2}}}{e^{ - {{\left( {\delta x} \right)}^2}/2{\sigma ^2}}}{\text{ }}{\text{}}.\end{align} \tag{ 50 }$$

Here, σ is the rms width of the wave packets. The direction of the shift is given by the angle between the velocity and the electric field in the dot product of (49). This provides a method of incorporating the displacement in time and space when using the semi-classical (instantaneous) approximation for the scattering process.

Generally, the collision duration may be considered to be exponentially distributed, and then may be estimated by

$$\begin{align}{\tau _c} = - {\tau^{\prime} _c}ln\left( r \right){\text{ ,}}\end{align} \tag{ 51 }$$

where ${\tau^{\prime} _c}$ is the maximum collision duration determined from direct calculations [272] and r is a random number between 0 and 1. In figure 29, the average displacement during the collision is plotted for a simulation of a resonant tunneling diode [277], with the tunneling structure centered around x = 0. The RTD structure contains 3 nm barriers of 0.3 eV height and a 5 nm well. It may be seen that, away from the tunneling region, the rms displacement is about 0.6 nm. But in the tunneling region, the rms displacement can be much larger. This is because the field is considerably larger around the tunneling structure (due to the barriers), where most of the potential drop occurs. Here, the displacement is peaked in the direction of the electric field, and clearly shows the increase due to the intra-collisional field effect.

Zoom In Zoom Out Reset image size