Self-Consistent Solution of Schrödinger-Poisson Equations For Modeling Realistic Silicon Transistors

1. Introduction

With continued scaling of semiconductor devices to meet the aggressive goals set forth by the International Technology Roadmap for Semiconductors (ITRS) [1], cost of fabrication are exponentially increasing, adding equivalent value to the power of predictive analysis through accurate modeling. Moreover, due to aggressive scaling of device dimensions, quantum transport is gradually becoming the dominant transport mechanism for carriers. The current generation of transistors [2,3,4,5,6,7,8,9] have active regions of length comparable to the wavelength of an electron. At such small dimensions, the basic assumption of the semi-classical transport physics [10,11], that electrons are point-like objects with defined position and momentum, falls short. A full quantum study of electron transport is needed to simulate and understand the behavior of these devices, including explicit quantum mechanical effects like tunneling, quantization and resonance [12,13,14,15].

The three well-known methods to study quantum transport are the Wigner function method [16,17,18], the Non-equilibrium Green’s function method [19,20,21,22,23] and the self-consistent solution of the Schrödinger, Pauli master [24,25,26] and Poisson equations. In our research we will be using the latter to model ballistic and dissipative electron transport. The model we employ uses the single-electron envelope approximation and the conduction bands are approximated with parabolic ellipsoidal or spherical valleys (six ellipsoidal valleys in the case of silicon). In particular, we are extending the popular Quantum Transmitting Boundary method (QTBM) for solving the two-dimensional Schrödinger equation and therefore validate its application towards behavior analysis of more realistic devices.

2. Simulation Methodology

Our main objective is to solve a system that can communicate with its environment through its contacts, namely an open system [27]. This requires discretizing the continuum of energies in the contacts. The optimal discretization is obtained by sampling the density-of-states using the solution of the system in equilibrium, namely a closed system [28,29]. Section 2.1 details our theoretical scheme of obtaining these solutions of the closed system. We solve the Schrödinger equation with closed boundary conditions self-consistently [30,31] with the Poisson equation in the two-dimensional (2-D) plane of the device to achieve this. The device is assumed wide enough so that a 3-D simulation is not required. In Section 2.3, our simulation methodology is extended to model an open system driven far from equilibrium under the influence of a drain-to-source bias. The Quantum Transmitting Boundary method (QTBM) [32,33] is used to enforce open boundary conditions to model particle exchange, through current-carrying leads, between the device and environment. Furthermore, ways to derive various device characteristics of practical importance are highlighted here. Finally we provide an application for our ballistic simulation tool.

We start with the equilibrium condition of the system. This is obtained by solving the two-dimensional (2-D) Schrödinger and Poisson equations self consistently [30] for different gate biases and at different temperatures, with no drain-to-source voltage applied. For very wide devices, we can assume translational invariance of the device profile along the out-of-plane y direction and therefore our 2-D approach is a good approximation. We define this closed state of the system as one in which the system does not exchange particles (and energy, thermodynamically) with the rest of the `universe’ through its contacts. Therefore, in its quantum description, either the wavefunctions must vanish outside the system or there should be zero net probability current entering or exiting the system. In order to obtain the correct closed system solution, it is better to solve the Schrödinger equation with both Dirichlet and Neumann boundary conditions, so that one gets the right density near the contacts and, so, the correct potential. The Dirichlet and Neumann solutions behave `sine-like’ and `cosine-like’ at the device-contact interface, respectively, as Figure 1 shows.

Thereby, a combined set of these two forms a complete set of physical solutions and any other solution to the closed system can be represented as a linear combination of the two. In all the cases, zero-value Dirichlet boundary conditions will be applied at the oxide-semiconductor interface as well as at the substrate of the device to account for the idealized yet physically consistent situation of having zero charges beyond the said interfaces.

Next, we describe how to obtain a self-consistent 2-D solution of the Schrödinger equation, considering both Dirichlet (Sec. 2.1) and QTBM boundary conditions (Sec. 2.3), and Poisson equations. When solving the Poisson equation, in order to account for electrostatic control through gates using externally applied bias, Dirichlet boundary conditions are applied at the gate-dielectric interfaces, while Neumann conditions are imposed elsewhere to ensure confinement of the electric field to the simulation domain.

2.1. Two-Dimensional Schrödinger Equation Solution for Quasi-static Charge Distribution

The time-independent single electron `effective mass’ Luttinger-Kohn (Schrödinger) equation [34] is solved using the parabolic-band approximation to account for the Silicon band structure. The numerical method of choice is the finite differences method. The Luttinger-Kohn equation can be written as:

$$[E_m(-i\nabla )+V(x,y,z)]\psi (x,y,z)=E\psi (x,y,z),$$

(1)

where, using the effective mass approximation, $E_m\left(k\right)=\frac{{\hslash }^2{k}^2}{2m^{*}}$ and $m^{*}$ is the effective mass of band m. For Silicon channels oriented along the $<100>$ direction, the two-dimensional Schrödinger equation can be written as [30,35]:

$$-\frac{{\hslash }^2}{2}\left[\frac{1}{m_{\mathrm{x}}^v}\frac{{\partial }^2{\xi }^v(x,z)}{\partial x^2}+\frac{1}{m_{\mathrm{z}}^v}\frac{{\partial }^2{\xi }^v(x,z)}{\partial z^2}\right]+V(x,z){\xi }^v(x,z)=E_{xz}^v{\xi }^v(x,z),$$

(2)

where the envelope wavefunction is $\psi (x,y,z)=e^{\mathrm{i}{k}_{\mathrm{y}}y}\xi (x,z)$ and $m_{\mathrm{x}}$, $m_{\mathrm{y}}$ and $m_{\mathrm{z}}$ are the transport mass (along transport x direction), the out-of-plane mass (along the y-direction into the plane of the device), and the quantization mass (along z, perpendicular to the carrier transport). V is the potential energy on the $(x,z)$ plane, $E_{xz}$ is the energy of the two-dimensional wavefunction $\xi (x,z)$ and $k_{\mathrm{y}}$ represents the wavevector in the y direction. The superscript v is used to represent the different values taken by the annonated quantities for the six ellipsoidal valleys of the Si Brillouin zone All computations are repeated once for each pair of inequivalent valley-orientations.

The expression $E_{xz}^v+{\hslash }^2{k}_{\mathrm{y}}^2/2m_{\mathrm{y}}^v$ denotes the total energy of an electron state described by the wavefunction $\psi (x,y,z)$. Numerically, Eq. (2) can be recast as:

$$\begin{array}{c}{a}_1\left[\frac{{\xi }^v(x_{i+1},z_i)}{\Delta x_i^2}+\frac{{\xi }^v(x_{i-1},z_i)}{\Delta x_i\Delta x_{i-1}}\right]+a_2\left[\frac{{\xi }^v(x_i,z_{i+1})}{\Delta z_i^2}+\frac{{\xi }^v(x_i,z_{i-1})}{\Delta z_i\Delta z_{i-1}}\right]-[\frac{a_1}{\Delta x_i^2}+\frac{a_1}{\Delta x_i\Delta x_{i-1}}\hfill \\ +\frac{a_2}{\Delta z_i^2}+\frac{a_2}{\Delta z_i\Delta z_{i-1}}+V(x_i,z_i)]{\xi }^v(x_i,z_i)=E_{xz}^v{\xi }^v(x_i,z_i),\end{array}$$

(3)

where $a_1=-{\hslash }^2/2m_{\mathrm{x}}^v$, $a_2=-{\hslash }^2/2m_{\mathrm{z}}^v$ and, i is the index for the N discretized mesh points. We deploy a 2-D mesh obtained by the tensor product of two 1-D meshes consisting of $N_{\mathrm{x}}$ and $N_{\mathrm{z}}$ points, so that $N=N_{\mathrm{x}}\times N_{\mathrm{z}}$. The $N\times N$ Hamiltonian matrix (${\mathbf{H}}^{\mathrm{D}}$) is created using the terms on the left-hand-side of Eq. 3. Since we are modeling a closed system, the wavefunctions should vanish outside the device and therefore Dirichlet boundary conditions are assumed at the simulation boundaries. With these assumptions, Eq. (3) takes the form of an eigenvalue problem of rank N:

$$\left[{\mathbf{H}}^{\mathrm{D}}+\mathbf{V}\right]·{\mathbf{\xi }}_{\mu }=E_{\mu }{\mathbf{\xi }}_{\mu },$$

(4)

where ${\mathbf{\xi }}_{\mu }$ and $E_{\mu }$ are the ${\mu }^{\mathrm{th}}$ eigenfunction and eigenvalue, respectively, and, $\mu$ ranges from 1 to N. Incorporating non uniformity of the mesh, the Hamiltonian matrix (${\mathbf{H}}^{\mathrm{D}}$) is given by:

$${\mathbf{H}}^{\mathrm{D}}={\mathbf{H}}_{\mathrm{x}}^{\mathrm{D}}+{\mathbf{H}}_{\mathrm{z}}^{\mathrm{D}},$$

(5)

where ${\mathbf{H}}_{\mathrm{x}}^{\mathrm{D}}=$

$$\left[\begin{array}{cccccccc}-2{\mathrm{A}}_{1,1}^{1,1}& {\mathrm{A}}_{1,2}^{1,1}& \cdots & 0& \cdots & \cdots & \cdots & 0\\ {\mathrm{A}}_{2,1}^{2,1}& \begin{array}{c}-{\mathrm{A}}_{2,2}^{2,2}\\ -{\mathrm{A}}_{2,1}^{2,2}\end{array}& {\mathrm{A}}_{2,2}^{2,3}& 0& & \cdots & \cdots & 0\\ 0& \ddots & \ddots & \ddots & \cdots & & & 0\\ ⋮& & \ddots & \ddots & \ddots & \cdots & & ⋮\\ 0& \cdots & {\mathrm{A}}_{\mathrm{i},\mathrm{i}-1}^{\mathrm{i},\mathrm{i}-1}& \begin{array}{c}-{\mathrm{A}}_{\mathrm{i},\mathrm{i}}^{\mathrm{i},\mathrm{i}}\\ -{\mathrm{A}}_{\mathrm{i},\mathrm{i}-1}^{\mathrm{i},\mathrm{i}}\end{array}& {\mathrm{A}}_{\mathrm{i},\mathrm{i}}^{\mathrm{i},\mathrm{i}+1}& 0& \cdots & 0\\ 0& \cdots & \cdots & \ddots & \ddots & \ddots & \cdots & 0\\ 0& \cdots & \cdots & & & & {\mathrm{A}}_{\mathrm{N},\mathrm{N}-1}^{\mathrm{N},\mathrm{N}-1}& -2{\mathrm{A}}_{\mathrm{N},\mathrm{N}}^{\mathrm{N},\mathrm{N}}\end{array}\right],$$

and, ${\mathbf{H}}_{\mathrm{z}}^{\mathrm{D}}=$

$$\left[\begin{array}{ccccccc}-2{\mathrm{B}}_{1,1}^{1,1}& 0& \cdots & {\mathrm{B}}_{1,1}^{1,1+{\mathrm{N}}_{\mathrm{x}}}& 0& \cdots & \\ 0& -2{\mathrm{B}}_{2,2}^{2,2}& 0& \cdots & {\mathrm{B}}_{2,2}^{2,2+{\mathrm{N}}_{\mathrm{x}}}& 0& \cdots \\ 0& \cdots & \ddots & & & \cdots \\ ⋮& & \cdots & \ddots & & \cdots \\ & & {\mathrm{B}}_{\mathrm{i},\mathrm{i}-{\mathrm{N}}_{\mathrm{x}}}^{\mathrm{i},\mathrm{i}-{\mathrm{N}}_{\mathrm{x}}}& \cdots & \begin{array}{c}-{\mathrm{B}}_{\mathrm{i},\mathrm{i}}^{\mathrm{i},\mathrm{i}}\\ -{\mathrm{B}}_{\mathrm{i},\mathrm{i}-{\mathrm{N}}_{\mathrm{x}}}\end{array}& 0& \cdots & {\mathrm{B}}_{\mathrm{i},\mathrm{i}}^{\mathrm{i},\mathrm{i}+{\mathrm{N}}_{\mathrm{x}}}& \cdots \\ ⋮& \cdots & & & & \ddots & \cdots & \\ 0& & \cdots & & & & & {\mathrm{B}}_{\mathrm{N},\mathrm{N}-{\mathrm{N}}_{\mathrm{x}}}^{\mathrm{N},\mathrm{N}-{\mathrm{N}}_{\mathrm{x}}}& {\mathrm{B}}_{\mathrm{N},\mathrm{N}}^{\mathrm{N},\mathrm{N}}\end{array}\right].$$

${\mathrm{A}}_{i,j}^{k,l}=\frac{a_1}{\Delta x_i\Delta x_j}$ and ${\mathrm{B}}_{i,j}^{k,l}=\frac{a_2}{\Delta z_i\Delta z_j}$ are the matrix elements in row k and column l. ${\mathbf{H}}^{\mathrm{D}}$ actually gives wavefunctions that vanish just outside the device domain. All elements, except those on the diagonal, of the rows representing the left and right edges of the device, i.e., rows having indexes $j{N}_{\mathrm{x}}$ or $j{N}_{\mathrm{x}}$ (j being an integer within range $(0,N_{\mathrm{y}})$), are set to zero. This ensures application of zero-value Dirichlet boundary condition exactly at the left and right edge of the device. $\mathbf{V}$ is a diagonal matrix containing the potential distribution, as shown below,

$$\mathbf{V}=\left[\begin{array}{c}{V}_1\\ & V_2\\ & & \ddots \\ & & & \ddots \\ & & & & V_{N-1}\\ & & & & & V_N\end{array}\right],$$

(6)

where the index i represents each mesh point.

Typically the <110> channel orientation is preferred for current-generation MOSFET devices due to added advantage/favorability towards hole mobility, while having relatively small negative impact on electron mobility [36]. In this case, the effective mass tensor becomes a $3\times 3$ matrix, instead of a relatively simpler diagonal matrix, resulting in the presence of mixed second-order derivatives in the Schrödinger equation.

Let $(x^{\prime },y^{\prime },z^{\prime })$ and $(x,y,z)$ be the co-ordinate systems based on the $\left[110\right]$ and $\left[100\right]$ orientations, respectively. Based on simple geometry, we can write

$$x^{\prime }=\frac{1}{\sqrt{2}}(x+y)y^{\prime }=\frac{1}{\sqrt{2}}(x-y)z^{\prime }=z.$$

(7)

$$\mathrm{Re}-\mathrm{writing},x=\frac{1}{\sqrt{2}}(x^{\prime }+y^{\prime })y=\frac{1}{\sqrt{2}}(x^{\prime }-y^{\prime })z=z^{\prime }.$$

$$\mathrm{Now},\frac{\partial \xi }{\partial x}=\frac{\partial \xi }{\partial x^{\prime }}\frac{\partial x^{\prime }}{\partial x}+\frac{\partial \xi }{\partial y^{\prime }}\frac{\partial y^{\prime }}{\partial x}=\frac{1}{\sqrt{2}}(\frac{\partial \xi }{\partial x^{\prime }}+\frac{\partial \xi }{\partial y^{\prime }}).$$

(8)

Rotating the co-ordinate system from $(x,y,z)$ to $(x^{\prime },y^{\prime },z^{\prime })$ by incorporating the above relations in Eq. (1) we obtain

$$\begin{array}{cc}-\frac{{\hslash }^2}{2}\left\{\frac{1}{2}\left[\frac{1}{m_{\mathrm{x}}^v}+\frac{1}{m_{\mathrm{y}}^v}\right]\left[\frac{{\partial }^2}{\partial x^{\prime 2}}+\frac{{\partial }^2}{\partial y^{\prime 2}}\right]+\frac{1}{m_{\mathrm{z}}^v}\frac{{\partial }^2}{\partial z^{\prime 2}}\right.\hfill & \\ \hfill +\left.\left[\frac{1}{m_{\mathrm{x}}}-\frac{1}{m_{\mathrm{y}}}\right]\frac{{\partial }^2}{\partial x^{\prime }\partial y^{\prime }}+V(x^{\prime },y^{\prime })\right\}& {\psi }^v(x^{\prime },y^{\prime },z)=E^v{\psi }^v(x^{\prime },y^{\prime },z^{\prime }).\hfill \end{array}$$

(9)

We can express the solution as ${\psi }^v(x^{\prime },y^{\prime },z^{\prime })=e^i{k}_{y^{\prime }}(\alpha x^{\prime }+y^{\prime }){\xi }^v(x^{\prime },z^{\prime })$ to remove the mixed 2nd order derivatives. Note here that we take into account the translational invariance of the potential in the $y^{\prime }$ direction by assuming a free electron solution in that direction. Substituting the guess solution in Eq. (9) we obtain,

$$\begin{array}{c}-\frac{{\hslash }^2}{2}[\frac{1}{2}(\frac{1}{m_{\mathrm{x}}^v}+\frac{1}{m_{\mathrm{y}}^v})(-{\alpha }^2{k}_{y^{\prime }}^2\xi (x^{\prime },z^{\prime })+2i\alpha k_{y^{\prime }}\frac{\partial {\xi }^v(x^{\prime },z^{\prime })}{\partial x^{\prime }}+\frac{{\partial }^2{\xi }^v(x^{\prime },z^{\prime })}{\partial x^{\prime 2}}-k_{y^{\prime }}^2)\hfill \\ +\frac{1}{m_{\mathrm{z}}^v}\frac{{\partial }^2{\xi }^v(x^{\prime },z^{\prime })}{\partial z^{\prime 2}}+(\frac{1}{m_{\mathrm{x}}^v}-\frac{1}{m_{\mathrm{y}}^v})(-\alpha k_{y^{\prime }}^2{\xi }^v(x^{\prime },z^{\prime })+i{k}_{y^{\prime }}\frac{\partial {\xi }^v(x^{\prime },z^{\prime })}{\partial x^{\prime }})\\ \hfill +V(x^{\prime },z^{\prime }){\xi }^v(x^{\prime },z^{\prime })]=E^v{\xi }^v(x^{\prime },z^{\prime }).\end{array}$$

(10)

Equating the coefficients of the first order derivatives w.r.t. x to 0,

$$i\alpha k_{y^{\prime }}(\frac{1}{m_{\mathrm{x}}^v}+\frac{1}{m_{\mathrm{y}}^v})+i{k}_{y^{\prime }}(\frac{1}{m_{\mathrm{x}}^v}-\frac{1}{m_{\mathrm{y}}^v})=0,$$

(11)

giving $\alpha =-\frac{m_{\mathrm{c}}^v}{2m_{\mathrm{xy}}^v}$, with $\frac{1}{m_{\mathrm{c}}^v}=\frac{1}{2}(\frac{1}{m_{\mathrm{x}}^v}+\frac{1}{m_{\mathrm{y}}^v})$ and $\frac{1}{m_{\mathrm{xy}}^v}=\frac{1}{m_{\mathrm{x}}^v}-\frac{1}{m_{\mathrm{y}}^v}$. Removing the `prime’ from the co-ordinate axes for convenience, the updated Schrödinger equation for representing transport in the $\left[110\right]$ direction then takes the form

$$\begin{array}{cc}-\frac{{\hslash }^2}{2}\left[\frac{1}{m_{\mathrm{c}}^v}\frac{{\partial }^2{\xi }^v(x,z)}{\partial x^2}+\frac{1}{m_{\mathrm{z}}^v}\frac{{\partial }^2{\xi }^v(x,z)}{\partial z^2}\right]+V(x,z){\xi }^v(x,z)\hfill & \\ \hfill =E^v{\xi }^v(x,z)-& \frac{{\hslash }^2{k}_{\mathrm{y}}^2}{2}\left[\frac{1}{m_{\mathrm{c}}^v}-\frac{m_{\mathrm{c}}^v}{4{\left(m_{\mathrm{xy}}^v\right)}^2}\right]{\xi }^v(x,z),\hfill \end{array}$$

(12)

where the full (envelope) wavefunction ${\psi }^v(x,y,z)$ is given as

$${\psi }^v(x,y,z)=e^{\mathrm{i}{k}_{\mathrm{y}}y}{e}^{-\mathrm{i}\frac{m_{\mathrm{c}}^v}{2m_{\mathrm{xy}}^v}{k}_{\mathrm{y}}x}{\xi }^v(x,z).$$

(13)

The electron wavefunctions ${\xi }_{\mu }(x,z)$ (labeled by the index $\mu$) and the corresponding eigenvalues $E_{xz}^{\mu }$ are determined by solving this eigenvalue problem.

The Schrödinger equation for <110> direction, Eq. (12), must be utilized to treat only four of the six conduction band Si valleys, since the in-plane rotation from the $<100>$ to the $<110>$ direction does not affect the other two valleys with $m_{\mathrm{z}}=m_{\mathrm{L}}$. Here $m_{\mathrm{L}}=0.91m_0$ is the longitudinal effective mass and $m_{\mathrm{T}}=0.19m_0$ (used later) represents the transverse effective mass ($m_{(}0)$ is the electron mass). These valleys can be treated using the relatively simpler expression given by Eq. (3). Note that the calculations, shown here for Si, can be extended for any other material with similar orientations of their conduction band minima (provided their band structure can considerably approximated by effective masses).

2.2. Determination of Carrier Concentration And Use in Poisson Equation

The electron wavefunctions obtained as eigen-solutions to the closed-system Hamiltonian are used to compute the electron charge distribution in the device. The information regarding the DoS along the x and z directions is already contained in these quantized 2-D wavefunctions. We have only to use the 1-D DoS representing the continuous energy spectrum along the ’homogenous’ out-of-plane y direction given by the expression

$$D_{1D}^v\left(E_{\mathrm{y}}\right)=\frac{1}{\pi \hslash }\sqrt{\frac{m_{\mathrm{y}}^v}{2E_{\mathrm{y}}}}.$$

(14)

The expression for the electron charge distribution is

$$n(x,z)=\sum _v\sum _{\mu }{\int }_{E_{xz}^{\mu ,v}}^{\infty }\mathrm{d}{E}_{\mathrm{y}}{D}_{1D}^v\left(E_{\mathrm{y}}\right)f(E_{\mathrm{y}},E_{\mathrm{F}}){\left|{\xi }_{\mu }^v(x,z)\right|}^2,$$

(15)

where ${\xi }_{\mu }$ and $E_{\mu }$ represent the wavefunction and energy of the electronic state $\mu$, respectively. Using Eq. (14), the above expression becomes

$$n(x,z)=\sum _v\sum _{\mu }\frac{1}{\pi \hslash }\sqrt{\frac{m_{\mathrm{y}}^v{k}_{\mathrm{B}}T}{2}}{F}_{-\frac{1}{2}}\left(\frac{E_{\mathrm{F}}-E_{xz}^{\mu ,v}}{k_{\mathrm{B}}T}\right){\left|{\xi }_{\mu }^v(x,z)\right|}^2,$$

(16)

Here $F_{-\frac{1}{2}}$ is the Fermi-Dirac integral of order -1/2, computed numerically using the Gauss-Legendre quadrature method [37].

The hole charge distribution, $p(x,y,z)$, is calculated semiclassically using the following expression based on 3-D density-of-states,

$$\begin{array}{c}\hfill p(x,y,z)=\frac{1}{2\sqrt{\pi }}{\left(\frac{2m_{\mathrm{h}}{k}_{\mathrm{B}}T}{\pi {\hslash }^2}\right)}^{\frac{3}{2}}{F}_{\frac{1}{2}}\left(\frac{V(x,z)-(E_{\mathrm{F}}+E_{\mathrm{g}})}{k_{\mathrm{B}}T}\right),\end{array}$$

(17)

where $m_{\mathrm{h}}=0.8m_0$ is the hole effective mass, the $E_{\mathrm{g}}$ is the band gap energy of silicon, and, $F_{\frac{1}{2}}$ is the Fermi-Dirac integral of order 1/2 [37].

Note that the inclusion of temperature as an input parameter in the device simulation is achieved through the temperature dependence of Eqs. (16) and (17).

Here we briefly describe the numerical solution of Poisson equation which constitutes an important part of most transport simulations. We solve the Poisson equation in the 2-D plane $(x,z)$ of the device. To incorporate correctly the electrostatic effects caused by transitions in the electric permittivity along the semiconductor-dielectric interface, we solve the generalized Poisson equation:

$$\nabla .\left[ϵ(x,z)\nabla V(x,z)\right]={\mathrm{e}}^2\left[p(x,z)-n(x,z)+N_{\mathrm{A}}(x,z)-N_{\mathrm{D}}(x,z)\right],$$

(18)

where e is the electron charge and ${ϵ}_{\mathrm{Si}}(x,z)$ is the permittivity of the material present at $(x,z)$. The centered finite-differences method is used to solve this linear system problem :

$$\mathbf{P}·\mathbf{V}=\mathbf{D},$$

(19)

where $\mathbf{P}$ is an $N\times N$ matrix representing the differential operators on the left hand side of Eq. (19), $\mathbf{V}$ and $\mathbf{D}$ are $N\times 1$ matrices expressing the potential energy distribution and the charge terms on the right hand side of Eq. (18), respectively, for each of the N mesh points. Neumann boundary conditions are employed along all the domain boundaries. Here, the gate potential $V_{\mathrm{GS}}$ is applied by setting Dirichlet boundary conditions. Mathematically this is done by assigning zeros to all the elements of the rows representing the said portion of the oxide-metal interface, except for the diagonal elements where unity is assigned. The gate potential in electron volts is added to the corresponding rows of the charge matrix $\mathbf{D}$. This assignment might lead to singularity issues while solving the linear system Eq. (19) computationally. The solution consists in multiplying the unity diagonal terms on $\mathbf{P}$ and the corresponding gate potential terms on $\mathbf{D}$ by a factor that has the same order of magnitude as the other non-zero terms on $\mathbf{P}$. The Schrodinger and Poisson equations are solved selfconsistently following the flow highlighted in Figure 2.

2.3. 2-D Schrödinger Equation with Quantum Transmitting Boundary Method: Device Under Bias

Our main scenario of interest is the `open’ system in which electrons can flow into and out of the active region of the semiconductor device. This essentially describes the behavior of a device in the presence of an applied drain-to-source bias $V_{\mathrm{DS}}$. The method we follow is the QTBM. [32]. The source and drain contacts are hypothesized as infinitely long leads culminating into the device, as shown schematically in Figure 3.

A distinct coordinate system $({\omega }_s,{\kappa }_s)$ is defined for each lead s, as shown in the figure. This is just a matter of convenience as the same derivations/expressions become applicable for all the leads. The potential is assumed to be constant along the direction ${\omega }_s$ while the potential profile along ${\kappa }_s$, is set equal to that along the lead-device interface. The wavefunction outside the system of leads and device is considered zero. Keeping these assumptions in mind, the wavefunctions within the leads can be broken into two independent components - traveling waves along ${\omega }_s$ and wavefunctions with a discretized energy spectrum along ${\kappa }_s$ due to the quantum confinement in that direction. The latter part is computed by solving the 1-D Schrödinger equation along ${\kappa }_s$,

$$-\frac{{\hslash }^2}{2m_{\mathrm{z}}^v}\frac{{\partial }^2{\phi }_m^{r,v}\left({\kappa }_s\right)}{\partial {\kappa }_s^2}+V_s\left({\kappa }_s\right){\phi }_m^{s,v}\left({\kappa }_s\right)=E_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right),$$

(20)

where $V_s\left({\kappa }_s\right)$ is the electron potential energy along ${\kappa }_s$ in lead s and $m_z$ is the quantization mass. Eq. (20) becomes an eigenvalue problem as a consequence of the zero-value Dirichlet boundary conditions imposed by the lead edges. ${\phi }_m^s$ and $E_m^s$ represent the $m_{\mathrm{th}}$ eigenstate and eigen-energy, respectively, in lead s of the one-dimensional (1-D) eigenvalue problem which can be written as,

$${\mathbf{H}}_{1\mathrm{D}}·{\mathbf{\phi }}_m^{s,v}=E_m^{s,v}{\mathbf{\phi }}_m^{s,v},$$

(21)

where the 1-D Hamiltonian matrix ${\mathbf{H}}_{1\mathrm{D}}$ is a $N_{\mathrm{x}}\times N_{\mathrm{x}}$ tridiagonal matrix, built as follows:

$${\mathbf{H}}_{1\mathrm{D}}=\left[\begin{array}{ccccccc}-2A_1& B_1& 0& \cdots & \cdots & \cdots & 0\\ B_2& -2A_2& B_2& 0& \cdots & \cdots & 0\\ 0& \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0& \cdots & \cdots & \cdots & B_{n-1}& -2A_{n-1}& B_{n-1}\\ 0& \cdots & \cdots & \cdots & 0& B_n& -2A_{n-1}\end{array}\right].$$

(22)

The computed eigenstates are normalized along the ${\kappa }_s$ direction and the total wavefunction inside lead s is then given by:

$$\begin{array}{c}{\Psi }_{\beta ,s}^v({\omega }_s,{\kappa }_s)=\sum _{m=1}^{N_s^v}\left[a_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)e^{\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}+b_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)e^{-\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}\right]\hfill \\ \hfill +\sum _{m=N_s^v+1}^{\infty }{b}_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)e^{k_{m,\beta }^{s,v}{\omega }_s}.\end{array}$$

(23)

Here, the expression inside the first summation represent $N_s^v$ traveling waves with energy (along the x and z directions) $E_{\beta }^v>E_m^{s,v}$, going into and reflecting out of the device, respectively, through the lead s. The index $\beta$ denotes the wavefunctions ${\Psi }_{\beta ,s}^v$ and the corresponding energies $E_{\beta }^v$ in the conduction band valley v. The expression within the second summation represents the gradually decaying modes with energy $E_{\beta }^v<E_m^{s,v}$. The coefficients $a_m^{s,v}$’s are chosen as inputs for the different waves traveling into the device, while the coefficients $b_{m,\beta }^{s,v}$ need to be determined. For all other leads $j\ne s$, the injection amplitudes $a_m^{j,v}=0$. The wavevectors $k_{m,\beta }^{s,v}$ for the traveling modes are given by: ${\left[2m_{\mathrm{x}}^v(E_{\beta }^v-E_m^{s,v})\right]}^{1/2}/\hslash$ and, for the evanescent modes, by: ${\left[2m_{\mathrm{x}}^v(E_m^{s,v}-E_{\beta }^v)\right]}^{1/2}/\hslash$. The energy $E_{\beta }^v$ will be referred to as the `injection energy’. Note that ${\omega }_s=-x\mathrm{or}x$ depending on the left or right contact, respectively.

The boundary conditions at the interface ${\Gamma }_s$ between the device and lead s dictate the continuity of both the wavefunction at the interface, ${\varphi }_{\beta ,s}^v{|}_{{\Gamma }_s}={\Psi }_{\beta ,s}^v({\omega }_s=0,{\kappa }_s)$, and, the normal derivative, $\nabla {\varphi }_{\beta }^{s,v}.\widehat{{\kappa }_s}{|}_{{\Gamma }_s}=\nabla {\Psi }_{\beta ,s}^v({\omega }_s=0,{\kappa }_s).\widehat{{\kappa }_s}$ for all $m\le N_s^v$. Here ${\varphi }_{\beta ,s}^v(x,z)$ denotes the wavefunction inside the device. Using Eq. (23) and combining the aforementioned boundary conditions together, we obtain:

$$\begin{array}{c}\nabla {\varphi }_{\beta }^{s,v}.\widehat{{\kappa }_s}{|}_{{\Gamma }_s}=\sum _{m=1}^{N_s^v}[a_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{\partial \left(e^{\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}\right)}{\partial {\omega }_s}{|}_{{\omega }_s=0}\hfill \\ \hfill +b_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{\partial \left(e^{-\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}\right)}{\partial {\omega }_s}{|}_{{\omega }_s=0}]+\sum _{m=N_s^v+1}^{\infty }{b}_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{\partial \left(e^{k_{m,\beta }^{s,v}}\right)}{\partial {\omega }_s}{|}_{{\omega }_s=0}.\end{array}$$

(24)

Using finite-differences to discretize the right hand side of the above equation, we obtain

$$\begin{array}{c}\frac{\partial {\varphi }_{\beta }^{s,v}}{\partial x}{|}_{{\omega }_s=0}=\sum _{m=1}^{N_s^v}[a_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{\mathrm{i}{k}_{m,\beta }^{s,v}({\omega }_s=0)}\left(e^{\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1\right)}{\Delta {\omega }_s}\hfill \\ +b_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{-\mathrm{i}{k}_{m,\beta }^{s,v}({\omega }_s=0)}\left(e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1\right)}{\Delta {\omega }_s}]\\ \hfill +\sum _{m=N_s^v+1}^{\infty }{b}_{m,\beta }^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{k_{m,\beta }^{s,v}({\omega }_s=0)}\left(e^{k_{m,\beta }^{s,v}\Delta {\omega }_s}-1\right)}{\Delta {\omega }_s}.\end{array}$$

(25)

Now, multiplying both sides of Eq. (23) by ${\phi }_{m^{\prime }}^{s,v*}\left({\kappa }_s\right)$ and integrating over ${\kappa }_s$ in the range $(0,d_s)$, where $d_s$ is the vertical height of lead s, we obtain for ${\omega }_s=0$,

$$\begin{array}{c}{\int }_0^{d_s}\mathrm{d}{\kappa }_s{\phi }_{m^{\prime }}^{s,v*}\left({\kappa }_s\right){\Psi }_{\beta ,s}^v({\omega }_s=0,{\kappa }_s)=\sum _{m=1}^{N_s^v}\left[a_m^{s,v}+b_{m,\beta }^{s,v}\right]{\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{{\mathrm{m}}^{\prime }}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}}\left({\kappa }_{\mathrm{s}}\right)\hfill \\ \hfill +\sum _{m=N_s^v+1}^{\infty }{b}_{m,\beta }^{s,v}{\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{{\mathrm{m}}^{\prime }}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}}\left({\kappa }_{\mathrm{s}}\right).\end{array}$$

(26)

The wavefunctions $\phi$ are eigenstates of a Hermitian operator (${\mathbf{H}}_{1\mathrm{D}}$) and so are mutually orthogonal, ${\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{{\mathrm{m}}^{\prime }}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}}\left({\kappa }_{\mathrm{s}}\right)={\delta }_{{\mathrm{mm}}^{\prime }}$. Using this relation and replacing ${\Psi }_{\beta }$ with ${\varphi }_{\beta }$, observing continuity of the wavefunction at the interface, Eq. (26) gives us

$$b_{m,\beta }^{s,v}={\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\varphi }_{\beta }^{\mathrm{s},\mathrm{v}}({\omega }_{\mathrm{s}}=0,{\kappa }_{\mathrm{s}})-{\mathrm{a}}_{\mathrm{m}}^{\mathrm{s},\mathrm{v}}.$$

(27)

Note that since evanescent waves are not injected, the above relation will not contain the $a_m^{s,v}$ term for $m>N_s^v$. Eq. (25) can be re-written as:

$$\begin{array}{cc}\hfill \frac{\partial {\varphi }_{\beta }^{s,v}}{\partial x}{|}_{{\omega }_s=0}& =\sum _{m=1}^{N_s^v}[a_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)\left(\frac{e^{\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}-\frac{e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}\right)\hfill \\ \hfill & +\left({\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\varphi }_{\beta }^{\mathrm{s},\mathrm{v}}({\omega }_{\mathrm{s}}=0,{\kappa }_{\mathrm{s}})\right){\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}]\hfill \\ \hfill & +\sum _{m=N_s^v+1}^{\infty }\left({\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\varphi }_{\beta }^{\mathrm{s},\mathrm{v}}({\omega }_{\mathrm{s}}=0,{\kappa }_{\mathrm{s}})\right){\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{k_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}\hfill \\ \hfill & =\sum _{m=1}^{N_s^v}[\frac{2\mathrm{i}{a}_m^{s,v}{\phi }_m^{s,v}\left({\kappa }_s\right)sin\left(k_{m,\beta }^{s,v}\Delta {\omega }_s\right)}{\Delta {\omega }_s}\hfill \\ \hfill & +\left({\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\varphi }_{\beta }^{\mathrm{s},\mathrm{v}}({\omega }_{\mathrm{s}}=0,{\kappa }_{\mathrm{s}})\right){\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}]\hfill \\ \hfill & +\sum _{m=N_s^v+1}^{\infty }\left({\int }_0^{d_s}\mathrm{d}{\kappa }_{\mathrm{s}}{\phi }_{\mathrm{m}}^{\mathrm{s},\mathrm{v}*}\left({\kappa }_{\mathrm{s}}\right){\varphi }_{\beta }^{\mathrm{s},\mathrm{v}}({\omega }_{\mathrm{s}}=0,{\kappa }_{\mathrm{s}})\right){\phi }_m^{s,v}\left({\kappa }_s\right)\frac{e^{k_{m,\beta }^{s,v}\Delta {\omega }_s}-1}{\Delta {\omega }_s}.\hfill \end{array}$$

(28)

In practice the wavefunctions $\varphi$ are calculated separately for each injected traveling wave and thus they depend on m as well (a superscript m will be used to denote this dependence henceforth). In our problem it has been seen that the extent of numerical errors is reduced with the introduction of more discretization into the system. For this reason, we also use a discretized version of the wave vectors $k_{m,\beta }^{s,v}$. The deduction of the same is given below in brief. Our starting point is the Schrödinger equation for the forward traveling wave $e^{i\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}$,

$$-\frac{{\hslash }^2}{2m_{\mathrm{x}}^v}\frac{{\partial }^2\left(e^{\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_r}\right)}{\partial x^2}=(E_{\beta }^v-E_m^{s,v})e^{\mathrm{i}{k}_{m,\beta }^{s,v}{\omega }_s}.$$

Using centered differences method, the second order derivative can be written as:

$$\begin{array}{c}\hfill -\frac{{\hslash }^2}{2m_{\mathrm{x}}^v}\frac{e^{\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}+e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-2}{\Delta {\omega }_s^2}=E_{\beta }^v-E_m^{s,v}\\ \hfill ,\mathrm{or},cos\left(k_{m,\beta }^{s,v}\Delta {\omega }_s\right)-1=-m_{\mathrm{x}}^v\Delta {\omega }_s^2(E_{\beta }^v-E_m^{s,v})/{\hslash }^2.\end{array}$$

Using basic trigonometric identities, the discretized wave vector can be retrieved as,

$$k_{m,\beta }^{s,v}=\frac{2}{\Delta {\omega }_s}arcsin\left(\sqrt{m_{\mathrm{x}}^v\Delta {\omega }_s^2(E_{\beta }^v-E_m^{s,v})/\left(2{\hslash }^2\right)}\right).$$

(29)

Similarly for the evanescent waves $(m>N_s^v)$, the wave vector can be written as:

$$k_{m,\beta }^{s,v}=\frac{2}{\Delta {\omega }_s}arcsinh\left(\sqrt{m_{\mathrm{x}}^v\Delta {\omega }_s^2(E_m^{s,v}-E_{\beta }^v)/\left(2{\hslash }^2\right)}\right).$$

(30)

The QTBM Hamiltonian used to calculate the wavefunctions ${\varphi }_{m,\beta }^{s,v}$ is built using the closed system Hamiltonian ${\mathbf{H}}^{\mathrm{D}}$, described in Sec. 2.1. The matrix elements corresponding to drain and source ends of the device are modified to capture the QTBM boundary conditions. Let $l^1,l^2,...,l^{N_{\mathrm{L}}}$ and $r^1,r^2,...,r^{N_{\mathrm{R}}}$ be mesh indices of the $N_{\mathrm{L}}$ ($N_{\mathrm{L}}\le N_{\mathrm{x}}$) and $N_{\mathrm{R}}$ ($N_{\mathrm{R}}\le N_{\mathrm{x}}$) points in the range $(1,N_{\mathrm{x}}\times N_{\mathrm{y}})$, representing the left and right contacts, respectively. Two $N\times N$ matrices ${\Sigma }_{\mathrm{L}}$ and ${\Sigma }_{\mathrm{R}}$ are built which include the reflected and transmitted (both traveling and evanescent) waves traveling into and out of the device, respectively, shown below,

$$\begin{array}{c}{\Sigma }_{\mathrm{L}}=\hfill \\ \hfill \left[\begin{array}{ccccccccc}0& \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0& \cdots & {\mathrm{A}}^{l^1,l^1}-{\mathrm{L}}^{l^1,l^1}& \cdots & -{\mathrm{L}}^{l^1,l^2}& \cdots & -{\mathrm{L}}^{l^1,l^j}& \cdots & 0\\ ⋮& \cdots & \cdots & \ddots & \cdots & \cdots & \cdots & \cdots & 0\\ ⋮& \cdots & -{\mathrm{L}}^{l^2,l^1}& \cdots & {\mathrm{A}}^{l^2,l^2}-{\mathrm{L}}^{l^2,l^2}& \cdots & -{\mathrm{L}}^{l^2,l^j}& \cdots & 0\\ 0& \cdots & \cdots & \cdots & \cdots & \ddots & \cdots & \cdots & 0\\ 0& \cdots & -{\mathrm{L}}^{l^j,l^1}& \cdots & -{\mathrm{L}}^{l^2,l^2}& \cdots & {\mathrm{A}}^{l^j,l^j}-{\mathrm{L}}^{l^j,l^j}& \cdots & 0\\ 0& \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \ddots & \end{array}\right]\end{array}$$

(31)

$$\begin{array}{c}{\Sigma }_{\mathrm{R}}=\hfill \\ \hfill \left[\begin{array}{ccccccccc}0& \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0& \cdots & {\mathrm{A}}^{r^1,r^1}+{\mathrm{L}}^{r^1,r^1}& \cdots & {\mathrm{L}}^{r^1,r^2}& \cdots & {\mathrm{L}}^{r^1,r^i}& \cdots & 0\\ ⋮& \cdots & \cdots & \ddots & \cdots & \cdots & \cdots & \cdots & 0\\ ⋮& \cdots & {\mathrm{L}}^{r^2,r^1}& \cdots & {\mathrm{A}}^{r^2,r^2}+{\mathrm{L}}^{r^2,r^2}& \cdots & {\mathrm{L}}^{r^2,r^i}& \cdots & 0\\ 0& \cdots & \cdots & \cdots & \cdots & \ddots & \cdots & \cdots & 0\\ 0& \cdots & {\mathrm{L}}^{r^i,r^1}& \cdots & {\mathrm{L}}^{r^2,r^2}& \cdots & {\mathrm{A}}^{r^i,r^i}+{\mathrm{L}}^{r^i,r^i}& \cdots & 0\\ 0& \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \ddots & \end{array}\right]\end{array}$$

(32)

where the A’s are defined in the same way as for Eq. (4) and

$$\begin{array}{c}{\mathrm{L}}^{p,q}=\frac{a_1}{\Delta x}[\sum _{m=1}^{N_s^v}{\phi }_m^{s,v*}({\kappa }_s=q){\varphi }_{m,\beta }^{s,v}({\omega }_s=0,{\kappa }_s=q)\Delta {\omega }_r{\phi }_m^{s,v}({\kappa }_s=p)\hfill \\ \times (e^{-\mathrm{i}{k}_{m,\beta }^{s,v}\Delta {\omega }_s}-1)+\sum _{m=N_s^v+1}^{\infty }{\phi }_m^{s,v*}({\kappa }_s=q){\varphi }_{m,\beta }^{s,v}({\omega }_s=0,{\kappa }_s=q)\Delta {\omega }_s\\ \hfill \times {\phi }_m^{s,v}({\kappa }_s=p)(e^{k_{m,\beta }^{s,v}\Delta {\omega }_s}-1)].\end{array}$$

(33)

Here, p and q represent the mesh indices $l^i$ (or $r^j$), i (or j) being an integer in the range $(1,N_{\mathrm{L}})$ (or $(1,N_{\mathrm{R}})$, respectively). We also define a $N\times 1$ matrix ${\mathbf{B}}^m$, expressing the component of the wave (traveling) injected into the device:

$${\mathbf{B}}^m=\left[\begin{array}{c}0\\ ⋮\\ {\mathrm{G}}^{l^1}\\ 0\\ ⋮\\ {\mathrm{G}}^{l^j}\\ 0\\ ⋮\end{array}\right],$$

(34)

where ${\mathrm{G}}^{l^j}=-a_1\left[2\mathrm{i}{a}_m^{s,v}{\phi }_m^{s,v}({\kappa }_s=l^j)sin\left(k_{m,\beta }^{s,v}\Delta {\omega }_s\right)\right]/\Delta x^2$. The resulting linear system can be written as:

$$\left[{\mathbf{H}}^{\mathrm{D}}+{\Sigma }_{\mathrm{L}}+{\Sigma }_{\mathrm{R}}+\mathbf{V}\right]{\mathbf{\varphi }}_{m,\beta }^{s,v}={\mathbf{B}}^m.$$

(35)

Here ${\mathbf{\varphi }}_{m,\beta }^{s,v}$ is a $N\times 1$ vector expressing the value of electron wavefunction ${\varphi }_{m,\beta ,s}$ at the $N^{\mathrm{th}}$ mesh point. Eq. (35) is solved to obtain ${\varphi }_{m,\beta ,s}$ separately for the leads s, the injection energies $E_{\beta }^v$ and the traveling modes m ( and valley v). Based on the computed wavefunctions, the transmission coefficient, local density-of-states (LDOS) and, most importantly, the current can be determined [26,38].

3. Application Towards Realistic Device Modeling

We simulate the electrical behavior of a Si ultra-thin body double-gate (DG) nMOS with channel length of 10 nm (Figure 4). The active region (simulation region) is 4 nm thick (${\tau }_{\mathrm{Si}}$) and 50 nm long with 1 nm (EOT ≈ 0.3 nm) oxide at each top and bottom gate. The significance of this study lies in the fact that our modeled device is essentially a cross-section of the current generation of nanosheet-based gate-all-around transistors [2,7,39]. In other words, our 2-D simulation based approach can help analyze the electrical behavior of 3-D devices of industrial relevance. The channel is intrinsically doped with highly doped n-type source and drain regions. The doping profile within the source/drain is assumed as a dual Gaussian profile with peaks located at the oxide-semiconductor interfaces. The Schrodinger and Poisson equations are solved self-consistently in the 2-D simulation space with our QTBM-based platform.

Figure 5 shows the simulated $I_{\mathrm{DS}}$-$V_{\mathrm{GS}}$ characteristics of the double-gate nMOS at 10 K and 300 K with equal $V_{\mathrm{GS}}$ applied at both gates. Characteristic CMOS behavior is observed with fast switching action represented by a low subthreshold slope (SS) of 64 mV/dec at 300 K, while at 10K the SS drops below this mark as expected. For device operation deep inside saturation, the charge-distribution plot in Figure 6(a) shows the phenomenon of channel inversion (effective reduction of the inversion-layer depth), whereas volume inversion is observed in the linear region of operation (Figure 6(b)). Due to this, the current density takes the shape of a single centered beam close to injection at the source and splits into two beams when flowing through the dual inverted channels at high gate bias, before it finally merges at the drain.

Thus, we show that our QTBM-based quantum transport modeling platform can model 2-D cross-sections of realistic-sized devices and allow in-depth understanding of their electrical behavior by taking into account explicit quantum transport.

4. Summary

In summary, we have built a novel simulation methodology by solving the Schrödinger and Poission equations self-consistently with QTBM boundary conditions on a 2-D device plane. This quantum-transport based methodology can allow more accurate electrical and predictive analyses of realistic semiconductor devices (like FinFET, Gate-all-around and CFET architectures) by taking into account explicit quantum transport while still having manageable computational cost requirement. As an application, we demonstrate the capability of our tool by modeling the electrical behavior of a double-gate MOSFET which essentially acts as a simple investigative study into the behavior of nanosheet-based transistors of current-generation technology nodes.