S2DS: Physics-based compact model for circuit simulation of two-dimensional semiconductor devices including non-idealities

There has been growing interest in two-dimensional (2D) semiconductor devices and circuits after recent developments in the fabrication and growth of 2D materials beyond graphene.¹ Layered semiconductors include transition metal dichalcogenides (TMDs) such as MoS₂ and the so-called “X-enes,” for example, phosphorene. The thickness of one monolayer (1L) of such 2D materials is defined as the layer spacing from X-ray diffraction of the bulk material, and it is typically less than 1 nm (for MoS₂, t_2D = 6.15 Å).² While 1 L of some materials (e.g., phosphorene) has been difficult to isolate and measure due to the lack of stability in ambient conditions, most sulfur-based TMDs such as MoS₂ are stable and can be grown as large-area monolayers with promising electrical properties.^3,4 Sub-nanometer channel thickness is expected to enable good gate control and to minimize short channel effects in field-effect transistors (FETs) with sub-10 nm gate lengths.^5,6 It also appears that the mobility (and therefore carrier mean free paths) of 2D materials is better preserved than that of bulk materials (e.g., Si) when the channel thickness is scaled below ∼3 nm.⁷ In addition, several simulation studies have recently suggested that 2D FETs could benefit from the larger band gap of 1L TMDs, reducing the off-state leakage current.^8,9

Although there remains much room for improvement in 2D devices, nascent efforts are already underway to integrate them into simple circuits with few transistors.^10–12 With recent reports of wafer-scale uniform growth and improved mobilities,^3,4 as well as smaller contact resistance,⁷ larger-scale circuits will soon become feasible. This push towards the applications of 2D FETs requires a platform to evaluate such devices at the circuit and system level. With many 2D materials being tested, there is also a growing requirement to analyze a larger material parameter space. In this context, Jiménez introduced a model that captured the essential physics of ideal 2D FETs, but did not include non-idealities and parasitics present in realistic devices.¹³ Other 2D FET models accounted for device electrostatics based on Poisson's equation with varying degrees of complexity but without the inclusion of high-field velocity saturation,¹⁴ and fringing fields, or self-heating and temperature-dependent effects.^14,15

In this paper, we describe the Stanford 2D Semiconductor (S2DS) transistor model, a physics-based and data-calibrated model to simulate circuits and systems made from 2D materials. S2DS has been implemented in Verilog-A (Ref. 16) and is freely available online.¹⁷ The model description provided in this paper focuses on monolayer 2D materials because they are most promising from a scaling point of view, but S2DS can also treat few-layer FETs if properly calibrated. The paper is organized as follows. First, we derive an analytic expression for the drain current and calibrate the model with existing experimental data. Then, we describe additional effects, including high-field velocity saturation, self-heating, contact resistance, and fringing fields that are essential to understand the realistic behavior of 2D FETs. Importantly, self-heating is found to play a strong role, particularly as the quality of 2D materials improves and the current densities achieved in experiments increase. The ultimate purpose of S2DS is to help researchers understand transistor measurements and to assist in designing further experiments. The SPICE-compatible model also enables design and benchmarking of large-scale circuits and systems composed of 2D devices with varying substrates of differing thermal properties (e.g., silicon versus glass or plastics).

Figure 1 shows the schematic of a FET based on a monolayer 2D semiconductor, with its main parasitic components. For the sake of concreteness, a typical device structure considered in this work includes the 1L material on an insulator (e.g., the bottom oxide (BOX) shown in Fig. 1), a top-gate oxide and gate metal stack on top. A finite underlap distance (L_U) appears between the edge of the top-gate and the source and drain contacts, respectively. The contacts themselves have a finite length L_C, which can play a role when this becomes comparable to the current transfer length.⁷ Beneath the bottom oxide, the substrate can be conductive as a doped Si wafer (which is often used as a back-gate in experiments), or insulating such as glass, quartz, or a flexible polymer.

To simplify the analysis, we present equations for n-type transistors (n-FET), considering electron transport in the channel and positive voltages (V_GS, V_DS > 0). Of course, these can be easily modified to simulate p-type transistors (p-FET), as is done for inverter simulations in the latter part of this study. We build on the model developed by Jiménez to derive the current-voltage (I–V) characteristics, with several extensions described below, particularly with respect to “extrinsic” transistor aspects, including self-heating, velocity saturation, and fringing fields.^13,18

First, we include the effect of the band structure for charge calculation. For example, in most TMDs (including MoS₂) two conduction band valleys may participate in transport, one at the K point and the other halfway along the K and Γ points of the Brillouin zone, sometimes labeled the Q valley.^19–21 Figure 2(a) shows the schematic of this band structure, where ΔE_KQ is the energy separation between the K and Q conduction band valleys. We calculate the charge density as

where f(E) is the Fermi-Dirac distribution and DOS_2D(E) is the 2D density of states corresponding to the lowest band. The Fermi energy E_F = qV_C, where V_C is the voltage across the quantum capacitance for the 2D channel and q is the elementary charge. Simplifying the charge density expression, we obtain n_2D = N_2D ln(1 + α) where α = exp[(qV_C − E₀)/(k_BT)], E₀ = E_G/2, and

where k_B is the Boltzmann constant and T is the average device temperature. We use the semiconductor mid-gap as our energy reference (E = 0) such that the conduction band energy is E₀ and the valence band energy is −E₀. g_K and g_Q are the degeneracy of the K and Q conduction valleys, respectively, and m_effK and m_effQ are their respective DOS effective masses. For MoS₂g_K = 2, g_Q = 6,^14,m_effK = 0.48m₀, and m_effQ = 0.57m₀.²² ΔE_KQ is the energy separation between K and Q conduction valleys (∼0.11 eV for monolayer MoS₂).^19–21 For the sake of simplicity, we take the band gap (E_G) to be the same as the photoluminescence (PL) gap, which is ∼1.85 eV for MoS₂ and ∼1.65 eV for WSe₂. However, we note that the true electronic band gap can be affected by the dielectric screening environment, strain, and proximity to grain boundaries.^20,23

Along with the free charge carriers, impurities (N_Dop) and traps (N_it) also contribute to the total channel charge (Q_ch) as

As shown in Fig. 2(a), we model the interface traps as acceptors, situated at an effective energy E_it below the conduction band, with an effective trap density D_it. To simplify the model, D_it is assumed here to be a delta function in energy, but this approach could be generalized. At a particular bias, the number of trapped carriers (N_it) is given by

The device electrostatics are guided by the distributed capacitive circuit model shown in Fig. 2(b).^13,24 The top and the back oxide capacitance are C_t (= ε_OX/t_OX) and C_b (=ε_BOX/t_BOX), respectively. ε_OX and ε_BOX are the dielectric constants, and t_OX and t_BOX are the oxide thicknesses of the top and bottom oxide, respectively. C_q is the quantum capacitance and C_it is the capacitance due to traps at the oxide-semiconductor interface, taken as a combination from both interfaces of the ultra-thin 2D channel. The quantum capacitance C_q and the trap capacitance C_it are given by

V_GS–V_GS0 and V_BS–V_BS0 are internal voltages from the top- and the back-gate, respectively. V_GS0 and V_BS0 are flatband voltages of the top- and back-gate, treated as fit parameters to the experimental data. (For example, if the top-gate metal work function is increased, V_GS0 will be higher.) The total charge in the 2D channel (Q_ch from Eq. (3)) and the quantum potentials at the source (V_Cs) and the drain (V_Cd) are calculated iteratively as discussed in  Appendix A, including doping and trapped charges. We neglect the channel depletion capacitance because the channel thickness is less than 1 nm. We note that the effect of the fringing field from the drain through the BOX can also be incorporated in our model by including an additional capacitance between the drain node and the channel in the circuit shown in Fig. 2(b).²⁵ In thicker multi-layer channels, the depletion capacitance should be accounted for, in a similar manner as it is done for silicon-on-insulator (SOI) transistors.²⁶ 

We solve for semi-classical drift transport to obtain an expression of the drain current (I_D = −I_S) for all transistor operating regions. The semi-classical approach is appropriate even for 2D FETs near 10 nm channel length, as the present-day experimental mean free path in monolayer 2D semiconductors such as MoS₂ is ∼2 to 3 nm (see Fig. S10 in the supplementary material of Ref. 3). Similarly, our approach should hold for channel widths greater than 10–20 nm, for which edge scattering effects can be safely ignored. (And all experimental data for 2D semiconductors are typically taken on micron-width devices, to obtain larger current drives.) Thus

Here μ is the carrier mobility, C_t and C_b are the top and the bottom oxide capacitances per unit area, and other variables are defined earlier. We recall that α is a function of V_C, and thus I_D is calculated as the difference of Eq. (6) evaluated at V_Cs and V_Cd. (The complete derivation is given in  Appendix A.) Gate and diffusion currents are not included in the present version of the S2DS model, and thus, leakage power will be underestimated. Nevertheless, this could be a reduced component in TMD FETs, which have larger band gaps than semiconductors such as Si and Ge.

When fitting to some (but not all) experimental data, we need to introduce a finite output resistance modeled by a fit parameter λ as I_D,eff = I_D(1 + λV_DS). However, sometimes λ is not needed, especially when fitting the model against long channel back-gated MoS₂ FETs.³ For fitting the model with experimental data on top-gate transistors, we used a finite value 0 < λ < 0.1.^11,27 We note that the current saturation region is also influenced by device self-heating, which is taken into account self-consistently, as we will discuss below.

With these considerations, Fig. 3 displays the I_D vs. V_GS curve for a few trap densities and trap energies. For large trap densities, a part of the gate voltage is used to charge the traps. As a result, less voltage is available to induce mobile charges in the channel, leading to smaller drain currents [Fig. 3(a)]. In Fig. 3(b), the V_GS at the kink in the I_D vs. V_GS curve is the voltage at which the Fermi energy is closest to the trap energy and charges a significant amount of traps. For large V_GS, all traps are charged, and the drain current remains constant for different trap energies. Similarly, Fig. 4(a) shows the impact of doping the channel material for an n-type 2D FET. Large doping shifts the flatband voltage in the negative direction, increasing the current.

The electron mobility in 2D materials depends on the vertical and the lateral electric field, as well as on the temperature. The dependence on vertical (gate) field comes in through the dependence on carrier density. Higher carrier density can partially screen scattering with ionized impurities and remote polar phonons,²⁸ and higher carrier density also raises the Fermi level, which changes the effective density of states for scattering. Classical “sixth-power law of thickness” surface roughness scattering presented in ultra-thin (<3 nm) bulk semiconductors such as Si²⁹ should not, in principle, affect 2D semiconductors without dangling bonds.²⁸ However, the vertical field could affect scattering with microscopic roughness of the gate dielectric, including the remote phonons mentioned above. The mobility dependence on lateral field mostly comprises high-field effects, i.e., drift velocity saturation. The temperature dependence of mobility comes in through scattering with intrinsic phonons (of the 2D material) and remote dielectric phonons.

Keeping the above considerations in mind, we fit the mobility behavior (at low lateral field) with the following semi-empirical relationship:

where μ₀ is the effective mobility at zero field and room temperature (T₀), T is the average device temperature, F_V is the vertical electric field, γ is a positive constant that depends on dominant phonons, and fitting parameters η and F_C depend on the material and the quality of the interface. Figure 4(b) compares this model with experimental data for 1 L WSe₂ (Ref. 27), and the fit provides η = 6.8 and F_C = 305 V/μm. Previous studies on 1L MoS₂ have observed η = 1.45 and F_C = 90 V/μm.¹⁴ The value of γ is also obtained by fitting the model to experimental data, yielding γ (bulk) = 2.6 (Ref. 30) and γ (1L) = 1 to 1.6 for electron mobility in MoS₂.^3,4

At high lateral field, the carrier drift velocity begins to saturate and the effective mobility decreases. We include this effect using a semi-empirical relation

where ξ is a fitting parameter with typical value around 2–4, F is the lateral electric field, and v_sat is the saturation velocity. Note that it is this μ which is then used when calculating the current in Eq. (6). The temperature dependence of the saturation velocity is incorporated similarly to the models for graphene and Si,^31,32 as v_sat = v₀/(1 + N_OP) where the OP (optical phonon) occupation is N_OP = 1/[exp(ℏω_OP/(k_BT)) −1]. Here, ℏω_OP is the OP energy and v₀ can be interpreted as the saturation velocity extrapolated to zero Kelvin. For MoS₂, the best fit against experimental high-field data (for monolayer MoS₂ grown by chemical vapor deposition on SiO₂) is obtained with ℏω_OP ≈ 30 to 40 meV and v₀ ≈ 2 to 3 × 10⁶ cm/s.³³ However, when modeling I–V curves of different devices in the literature, we must treat v₀ and ℏω_OP as fitting parameters.

Considering that 2D devices can carry high current densities,^3,34,35 unlike their organic counterparts in flexible and transparent electronics,³⁶ such 2D transistors can generate significant heat. We can model the FET self-heating by including a thermal resistance $(RTH)$ such that the average device temperature rise is ΔT = T − T₀ = P$RTH$⁠, where P = I_D(V_DS − 2I_DR_C) is the power input of the device without the contacts. As illustrated in the inset of Fig. 5(a), the total thermal resistance has three components: the thermal boundary resistance (TBR) between the channel and the bottom oxide (⁠$RB$ = $RTBR$/A), the spreading resistance of the bottom oxide $(RBOX)$⁠, and the spreading thermal resistance into the substrate $(Rsi)$⁠.^37,38 The thermal resistance per unit length is given as

Here, k_BOX and t_BOX are the thermal conductivity and thickness of the bottom oxide (BOX), respectively, and k_sub is the thermal conductivity of the substrate. The TBR per unit area is $RTBR$ ≈ 10⁻⁷ m² KW⁻¹ for monolayer MoS₂ on SiO₂,³⁹ and the “thermal area” of the device is A ≈ W(L_G + 2L_U). Due to heat spreading effects in the SiO₂, the effective thermal width at the SiO₂/Si interface can be approximated as W_eff ≈ W + 2t_BOX.

The thermal expression in Eq. (9) strictly only applies to “longer” channel devices, at least three times longer than the lateral thermal healing length (L_H) along the channel. For 1L MoS₂ on 90 nm SiO₂ on Si substrate, L_H = [k_2Dt_2D(W/g + $RTBR$⁠)]^1/2 ≈ 100 nm, if we take k_2D ≈ 80 Wm⁻¹ K⁻¹ as the in-plane thermal conductivity of MoS₂ at room temperature (k_2D could become a function of length in shorter devices.⁴⁰). For “longer” devices (with L_G + 2L_U > 3L_H), the thermal resistance is given simply by $Rth$ ≈ 1/[g(L_G + 2L_U)]. For “shorter” channel length devices (L_G + 2L_U < 3L_H), heat flow into the metal contacts becomes non-negligible and can be taken into account as described in Ref. 38. We note that heat flow into a top metal gate can, in general, be neglected, partly due to TBR at the two top oxide interfaces, but mainly due to the presence of the larger heat sink (i.e., Si substrate) at the back-side.⁴¹ 

To quantify the impact of device self-heating, we simulate a typical 1L MoS₂ transistor (L_G = 1 μm) in Fig. 5(b) under four circumstances: without self-heating (solid line, top curve), with self-heating on 90 nm SiO₂/Si substrate (dashed), with self-heating on 300 nm SiO₂/Si substrate (dotted), and finally on a poor thermal substrate where the Si was replaced by a plastic like polyethylene naphthalate (PEN) or acrylic (dash-dotted). We observe a reduction in saturation current of approximately 20%, 26%, and 70% from the “ideal” scenario without any self-heating. Thus, self-heating becomes crucial for devices on substrates with poor thermal conductivity, especially when the FET channel is a high-quality 2D material which has large current-carrying capability.

We now discuss the key parasitic capacitances that contribute to C_GS (gate to source) and C_GD (gate to drain). The total parasitic capacitance between the gate and the channel nodes (source and drain) is due to internal fields through the channel (C_if), outer fringing fields through the surrounding region (C_of) and normal fringing fields between the gate, and the source or the drain (C_nf). We display these fields in the schematic shown in the inset of Fig. 5(a). The S2DS model does not include the capacitance between the gate and metal plugs at the drain or the source, but such capacitance can easily be included following Ref. 42.

The capacitance C_nf is obtained by mapping the perpendicular surfaces of the gate sidewall and the top surface of contact metal to equivalent parallel surfaces using conformal mapping.⁴³ We modify C_nf from Ref. 43 to only include the part of the gate sidewall (t_G + t_OX − t_C), which is higher than the contact metal as shown in the Fig. 5(a) inset. We assume that the contact length (L_C) is larger than the underlap length (L_U). C_of includes the fringing fields from the horizontal edges of the gate metal to the horizontal edges of the contact metal.⁴³ In addition, C_of includes the parallel capacitance between vertical sidewalls of the gate and the source or drain, approximated with an average distance (⁠$LU2+tOX2$⁠)^1/2 between sidewalls. By solving for the specific geometry in Fig. 5(a), we obtain the analytical form of the total fringing capacitance as

To obtain C_if, we separate the contribution of the channel charge between the source and the drain terminals by using the Ward-Dutton charge partition scheme⁴⁴ 

We also consider the impact of the fringing field from the top-gate on the carriers in the underlap region by including an effective fringing capacitor (C_tf) from the gate metal sidewall to the underlap region. The analytical expression for C_tf is obtained similarly to Eq. (10) using conformal mapping⁴³ 

We note that when t_OX ≪ t_G ≈ L_U, then ϕ ≈ 1.

Generally speaking, the total extrinsic resistance R_ext = R_U + R_C of the intrinsic device includes resistance due to the underlap region (R_U) and the contact resistance between metal and semiconductor (R_C). The underlap R_U can be reduced by adjusting the back-gate voltage (i.e., “electrostatic doping” using a back-gate plane under the entire device), or by chemical doping, the latter being preferred in realistic devices. The impact of both adjustments is shown in Fig. 6(a) and 6(b) for different underlap lengths. In addition, the underlap resistance can also be reduced by increasing the mobility of the 2D channel material.

The contact resistance (R_C) for 2D devices can display non-linear behavior with respect to drain and gate voltages due to the Schottky barrier at the metal-semiconductor interface.⁷ For simplicity, here we assume that R_C is optimized in the fabrication and is Ohmic, which is a good approximation at higher lateral field and high V_DS. Nonetheless, following work on organic transistors,^45,46 it is possible to model each contact with a pair of anti-parallel Schottky diodes.

In Fig. 6(c), we plot the fringe capacitance C_F (=C_of + C_nf) for different dielectric constants of the spacer region. As expected, the capacitance is largest for small L_U and drops off to a small value for larger underlap lengths. For this particular device geometry (see Fig. 6 caption), we plot R_U and C_F in Fig. 6(d). To minimize the R_UC_F product, the optimal L_U ≈ 25 nm in this case. If the underlap resistance is reduced (e.g., by higher doping or higher mobility in this region), the optimized underlap length will increase. Another solution to reduce the optimized L_U is to decrease the dielectric constant of surrounding material, ε_sp.

We calibrate the model to monolayer 2D semiconductors including MoS₂ (n-type) and WSe₂ (p-type) using the existing experimental data. We note that this calibration is not universal, i.e., it can be improved or adjusted as more data become available from improved nanofabrication of such devices. The model is also not restricted to monolayer devices, although thicker devices could require additional treatment of depletion capacitance²⁶ and band structure. Here, we use monolayer MoS₂ devices of various lengths ranging from L_G ≈ 80 nm to a few μm in order to calibrate the model.^3,11 The device contact resistance (e.g., ∼1 to 6 kΩ·μm depending on doping and back-gate voltage) and the mobility (20–40 cm² V⁻¹ s⁻¹) are independently measured through experimental transfer line method (TLM), reassuring our fit.³ 

The model provides a good fit for long-channel back-gated [Fig. 7(a)] as well as top-gated [Fig. 7(b)] monolayer MoS₂ devices. We also fit the model to long-channel 1 L WSe₂ (Ref. 27) and extract effective μ₀ = 245 cm² V⁻¹ s⁻¹, device doping (N_Dop ∼ 10¹³ cm⁻²), and trap density (D_it ∼ 10¹² cm⁻²). The model comparison to I_D–V_DS and I_D–V_GS experimental data for 1L WSe₂ p-type FET are shown in Figs. 7(c) and 7(d).

We note that the S2DS model can reproduce all experimental (I–V) data sets we have examined (from our lab and the published literature) with approximately ten fitting parameters. These include μ₀, β, ξ, V_GS0, as described above. Other material parameters such as carrier effective masses,²² thermal conductivities,⁴⁰ and band structure^19–21 are generally imported from first principles simulations or published experimental data and are not treated as fitting parameters. Nevertheless, these are available in the model and can be adjusted as future data become available (e.g., influence of strain on effective masses).

The S2DS compact model is implemented in Verilog-A (Ref. 16) to perform circuit simulations in SPICE, and the code is freely available online.¹⁷ Standard modeling guidelines were followed to capture the device behavior in Verilog-A.⁴⁷ Device self-heating is currently modeled in DC mode only, but time-dependent heating could also be modeled by including the device thermal capacitance, which will be dominated by the materials surrounding the 2D channel rather than the 2D channel itself.⁴¹ We utilize the multi-physics support in Verilog-A to include the lumped self-heating model using a thermal node. We use additional internal nodes to calculate the quantum potentials near the source and the drain (V_Cs and V_Cd).

Before concluding, we wish to illustrate the use of the S2DS model in simple circuit simulations, such as the inverters and ring oscillators in Fig. 8 (Additional system-level simulations, including the role of variability, will be the scope of future work.). For these circuit demonstrations, we calibrate n-type devices to 1L MoS₂ and p-type devices to 1L WSe₂. For MoS₂, we use μ₀ = 81 cm² V⁻¹ s⁻¹ (Ref. 48), and for WSe₂, we use μ₀ = 245 cm² V⁻¹ s⁻¹ (Ref. 27). For this test case, we consider channel lengths L_G = 100 nm, effective oxide thickness (EOT) = 2 nm, and contact resistance R_C = 200 Ω⋅μm. Over the range of gate voltages applied (supply voltage V_DD = 1 V), the mobility ratio of WSe₂ to MoS₂ is 3:1. Thus, we used a width ratio W_p/W_n = 1:3 to balance the inverter, where W_p = 1 μm and W_n = 3 μm are widths of FETs from WSe₂ and MoS₂, respectively. The saturation drive currents of the optimized devices are approximately ∼350 μA.

From our data-calibrated model, we predict a DC inverter gain of ∼130, as seen in Fig. 8(a). We also simulate a 3-stage ring oscillator [inset of Fig. 8(b)] in SPICE using the above-mentioned inverter designs. We load the individual inverters with a 30 fF load (C_L) to emulate the interconnect capacitance. The period of the oscillator is close to 0.5 ns for V_DD = 1 V as seen in Fig. 8(b). These results are not fundamental, but reliant on the particular test case chosen here, to illustrate the practicality of the S2DS model. Further improvements in mobility, drive currents, tailoring of flat-band voltages (V_GS0 and V_BS0) using different gate metals, and optimization of device dimensions (e.g., L_U) will further help to improve the delays and DC gain.

In summary, we described S2DS, a physics-based, predictive 2D transistor model capable of simulating devices and circuits from 2D semiconductors such as MoS₂ and WSe₂. S2DS is exclusively geared toward 2D materials, which cannot be readily treated with the existing models for ultra-thin body (UTB) silicon on insulator (SOI).^49–51 S2DS includes the quantum capacitance of 2D monolayers, anisotropic 2D material properties (e.g., thermal anisotropy), as well as the 2D band structure and mobility calibrated with the existing experimental data. S2DS also includes a thermal model, which is important for self-consistent high-field simulations, illustrating that 2D-FET saturation currents are thermally limited, especially on thermally insulating (e.g., flexible, plastic) substrates. The model is presented with sufficient generality that it can be easily updated as additional data become available for 2D materials and devices.

We also demonstrate the capability of S2DS to extract and quantify unknown variables such as doping density and trap density from electrical data. The model can be used to optimize device dimensions such as the gate-source underlap length, as well. Importantly, S2DS is implemented in Verilog-A, and the code is freely available on the nanoHUB.org.¹⁷ All the equations are analytical and can be comfortably used in the SPICE environment for multi-transistor circuit simulations. Such a compact model is a major step towards developing systems from 2D materials, and it could also be used to benchmark and select various 2D materials with future systems in mind.

The authors would like to thank Kirby Smithe, Chris English, Shaloo Rakheja, and Chi-Shuen Lee for sharing experimental data and useful discussion. This work was supported in part by the NCN-NEEDS program, which is funded by the National Science Foundation, contract 1227020-EEC, and by the Semiconductor Research Corporation (SRC). This study was also partly supported by the NSF EFRI 2-DARE grant 1542883, and by Systems on Nanoscale Information fabriCs (SONIC), one of the six SRC STARnet Centers, sponsored by MARCO and DARPA.

We solve for the dependence of channel charge (Q_ch, Eq. (3)), bias potentials, and oxide capacitances as follows:^13,24

The above equation solved iteratively with Eq. (3) provides the quantum potential V_C at the drain [calculated at V_n(x) = V_DS] and at the source [calculated at V_n(x) = 0]. C_t and C_b are the capacitances of the top and bottom oxide, respectively. V_GS and V_BS are the voltages at the top-gate and the back-gate with respect to the source terminal. V_GS0 and V_BS0 are the flat band voltages of the top- and back-gate, respectively, for an undoped channel, and are treated as fitting parameters. V_n(x) and V_C(x) are the channel potential and the quantum potential, respectively. Differentiating both sides of Eq. (A1) with respect to V_C(x), we obtain

where C_q and C_it are provided in Eq. (5). We assume drift transport in the channel,

By substituting F = −dV_n(x)/dx and integrating the equation over x from [0, L_G], we obtain

where we have assumed the mobilty to be constant along the channel. After changing the variable from x to channel potential V_n, the above equation is re-written as

Note that Q_n(V_n) [=−qn_2D(V_n)] is the charge due to the free carriers (electrons in this case). V_Cd and V_Cs are the quantum potentials at the drain and source, respectively. The final equation is derived by substituting capacitance expressions from Eqs. (A2) and (5). After simplifying, we obtain the following analytical form of the drain current

We use the Ward-Dutton partition⁴⁴ to simplify the expressions in Eq. (11)

The above expressions (B1 and B2), although derived for bulk MOSFETs,⁴⁴ have been previously used to calculate charge at the drain and the source nodes in carbon nanotube FETs⁵² as well as graphene FETs.⁵³ The symbol θ is an empirical function of V_DS, which is used to provide a continuum between the different regions of transistor operation. From the linear regime θ = 1 to the saturation regime θ = 0.^53,θ = 1 − F_q such that F_q = (V_DS/V_DSsat)/[(V_DS/V_DSsat)^a + 1]^1/a, with a typical value a = 3 to 4, to make θ continuous. V_DSsat is the saturation voltage at which the channel charge near the drain becomes zero or n_2D (at x = L_G) = 0. It is approximated as V_DSsat = Q_n0/[WL_GC_t(1 + ω)]. ω is a fit parameters ranging from 0 to 1.⁴⁴