We present a physics-based compact model for two-dimensional (2D) field-effect transistors (FETs) based on monolayer semiconductors such as MoS2. A semi-classical transport approach is appropriate for the 2D channel, enabling simplified analytical expressions for the drain current. In addition to intrinsic FET behavior, the model includes contact resistance, traps and impurities, quantum capacitance, fringing fields, high-field velocity saturation, and self-heating, the latter being found to play an important role. The model is calibrated with state-of-the-art experimental data for n- and p-type 2D-FETs, and it can be used to analyze device properties for sub-100 nm gate lengths. Using the experimental fit, we demonstrate the feasibility of circuit simulations using properly scaled devices. The complete model is implemented in SPICE-compatible Verilog-A, and a downloadable version is freely available at the nanoHUB.org.
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.1 Layered semiconductors include transition metal dichalcogenides (TMDs) such as MoS2 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 MoS2, t2D = 6.15 Å).2 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 MoS2 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.7 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,7 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.13 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,14 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.17 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).
II. COMPACT MODEL
A. Intrinsic I-V characteristics
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 (LU) appears between the edge of the top-gate and the source and drain contacts, respectively. The contacts themselves have a finite length LC, which can play a role when this becomes comparable to the current transfer length.7 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 (VGS, VDS > 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 MoS2) 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 ΔEKQ 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 DOS2D(E) is the 2D density of states corresponding to the lowest band. The Fermi energy EF = qVC, where VC is the voltage across the quantum capacitance for the 2D channel and q is the elementary charge. Simplifying the charge density expression, we obtain n2D = N2D ln(1 + α) where α = exp[(qVC − E0)/(kBT)], E0 = EG/2, and
where kB 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 E0 and the valence band energy is −E0. gK and gQ are the degeneracy of the K and Q conduction valleys, respectively, and meffK and meffQ are their respective DOS effective masses. For MoS2gK = 2, gQ = 6,14,meffK = 0.48m0, and meffQ = 0.57m0.22 ΔEKQ is the energy separation between K and Q conduction valleys (∼0.11 eV for monolayer MoS2).19–21 For the sake of simplicity, we take the band gap (EG) to be the same as the photoluminescence (PL) gap, which is ∼1.85 eV for MoS2 and ∼1.65 eV for WSe2. 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 (NDop) and traps (Nit) also contribute to the total channel charge (Qch) as
As shown in Fig. 2(a), we model the interface traps as acceptors, situated at an effective energy Eit below the conduction band, with an effective trap density Dit. To simplify the model, Dit 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 (Nit) 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 Ct (= εOX/tOX) and Cb (=εBOX/tBOX), respectively. εOX and εBOX are the dielectric constants, and tOX and tBOX are the oxide thicknesses of the top and bottom oxide, respectively. Cq is the quantum capacitance and Cit 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 Cq and the trap capacitance Cit are given by
VGS–VGS0 and VBS–VBS0 are internal voltages from the top- and the back-gate, respectively. VGS0 and VBS0 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, VGS0 will be higher.) The total charge in the 2D channel (Qch from Eq. (3)) and the quantum potentials at the source (VCs) and the drain (VCd) 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).25 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.26
We solve for semi-classical drift transport to obtain an expression of the drain current (ID = −IS) 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 MoS2 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, Ct and Cb are the top and the bottom oxide capacitances per unit area, and other variables are defined earlier. We recall that α is a function of VC, and thus ID is calculated as the difference of Eq. (6) evaluated at VCs and VCd. (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 ID,eff = ID(1 + λVDS). However, sometimes λ is not needed, especially when fitting the model against long channel back-gated MoS2 FETs.3 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 ID vs. VGS 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 VGS at the kink in the ID vs. VGS curve is the voltage at which the Fermi energy is closest to the trap energy and charges a significant amount of traps. For large VGS, 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,28 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 Si29 should not, in principle, affect 2D semiconductors without dangling bonds.28 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 μ0 is the effective mobility at zero field and room temperature (T0), T is the average device temperature, FV is the vertical electric field, γ is a positive constant that depends on dominant phonons, and fitting parameters η and FC depend on the material and the quality of the interface. Figure 4(b) compares this model with experimental data for 1 L WSe2 (Ref. 27), and the fit provides η = 6.8 and FC = 305 V/μm. Previous studies on 1L MoS2 have observed η = 1.45 and FC = 90 V/μm.14 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 MoS2.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 vsat 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 vsat = v0/(1 + NOP) where the OP (optical phonon) occupation is NOP = 1/[exp(ℏωOP/(kBT)) −1]. Here, ℏωOP is the OP energy and v0 can be interpreted as the saturation velocity extrapolated to zero Kelvin. For MoS2, the best fit against experimental high-field data (for monolayer MoS2 grown by chemical vapor deposition on SiO2) is obtained with ℏωOP ≈ 30 to 40 meV and v0 ≈ 2 to 3 × 106 cm/s.33 However, when modeling I–V curves of different devices in the literature, we must treat v0 and ℏωOP as fitting parameters.
C. Self-heating model
Considering that 2D devices can carry high current densities,3,34,35 unlike their organic counterparts in flexible and transparent electronics,36 such 2D transistors can generate significant heat. We can model the FET self-heating by including a thermal resistance such that the average device temperature rise is ΔT = T − T0 = P, where P = ID(VDS − 2IDRC) 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 ( = /A), the spreading resistance of the bottom oxide , and the spreading thermal resistance into the substrate .37,38 The thermal resistance per unit length is given as
Here, kBOX and tBOX are the thermal conductivity and thickness of the bottom oxide (BOX), respectively, and ksub is the thermal conductivity of the substrate. The TBR per unit area is ≈ 10−7 m2 KW−1 for monolayer MoS2 on SiO2,39 and the “thermal area” of the device is A ≈ W(LG + 2LU). Due to heat spreading effects in the SiO2, the effective thermal width at the SiO2/Si interface can be approximated as Weff ≈ W + 2tBOX.
The thermal expression in Eq. (9) strictly only applies to “longer” channel devices, at least three times longer than the lateral thermal healing length (LH) along the channel. For 1L MoS2 on 90 nm SiO2 on Si substrate, LH = [k2Dt2D(W/g + )]1/2 ≈ 100 nm, if we take k2D ≈ 80 Wm−1 K−1 as the in-plane thermal conductivity of MoS2 at room temperature (k2D could become a function of length in shorter devices.40). For “longer” devices (with LG + 2LU > 3LH), the thermal resistance is given simply by ≈ 1/[g(LG + 2LU)]. For “shorter” channel length devices (LG + 2LU < 3LH), 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.41
To quantify the impact of device self-heating, we simulate a typical 1L MoS2 transistor (LG = 1 μm) in Fig. 5(b) under four circumstances: without self-heating (solid line, top curve), with self-heating on 90 nm SiO2/Si substrate (dashed), with self-heating on 300 nm SiO2/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.
D. Capacitance modeling
We now discuss the key parasitic capacitances that contribute to CGS (gate to source) and CGD (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 (Cif), outer fringing fields through the surrounding region (Cof) and normal fringing fields between the gate, and the source or the drain (Cnf). 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 Cnf 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.43 We modify Cnf from Ref. 43 to only include the part of the gate sidewall (tG + tOX − tC), which is higher than the contact metal as shown in the Fig. 5(a) inset. We assume that the contact length (LC) is larger than the underlap length (LU). Cof includes the fringing fields from the horizontal edges of the gate metal to the horizontal edges of the contact metal.43 In addition, Cof includes the parallel capacitance between vertical sidewalls of the gate and the source or drain, approximated with an average distance ()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 Cif, we separate the contribution of the channel charge between the source and the drain terminals by using the Ward-Dutton charge partition scheme44
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 (Ctf) from the gate metal sidewall to the underlap region. The analytical expression for Ctf is obtained similarly to Eq. (10) using conformal mapping43
We note that when tOX ≪ tG ≈ LU, then ϕ ≈ 1.
E. Extrinsic resistance
Generally speaking, the total extrinsic resistance Rext = RU + RC of the intrinsic device includes resistance due to the underlap region (RU) and the contact resistance between metal and semiconductor (RC). The underlap RU 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 (RC) 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.7 For simplicity, here we assume that RC is optimized in the fabrication and is Ohmic, which is a good approximation at higher lateral field and high VDS. 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 CF (=Cof + Cnf) for different dielectric constants of the spacer region. As expected, the capacitance is largest for small LU and drops off to a small value for larger underlap lengths. For this particular device geometry (see Fig. 6 caption), we plot RU and CF in Fig. 6(d). To minimize the RUCF product, the optimal LU ≈ 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 LU is to decrease the dielectric constant of surrounding material, εsp.
III. MODEL CALIBRATION
We calibrate the model to monolayer 2D semiconductors including MoS2 (n-type) and WSe2 (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 capacitance26 and band structure. Here, we use monolayer MoS2 devices of various lengths ranging from LG ≈ 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 cm2 V−1 s−1) are independently measured through experimental transfer line method (TLM), reassuring our fit.3
The model provides a good fit for long-channel back-gated [Fig. 7(a)] as well as top-gated [Fig. 7(b)] monolayer MoS2 devices. We also fit the model to long-channel 1 L WSe2 (Ref. 27) and extract effective μ0 = 245 cm2 V−1 s−1, device doping (NDop ∼ 1013 cm−2), and trap density (Dit ∼ 1012 cm−2). The model comparison to ID–VDS and ID–VGS experimental data for 1L WSe2 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 μ0, β, ξ, VGS0, as described above. Other material parameters such as carrier effective masses,22 thermal conductivities,40 and band structure19–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).
IV. VERILOG-A AND CIRCUIT SIMULATIONS
The S2DS compact model is implemented in Verilog-A (Ref. 16) to perform circuit simulations in SPICE, and the code is freely available online.17 Standard modeling guidelines were followed to capture the device behavior in Verilog-A.47 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.41 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 (VCs and VCd).
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 MoS2 and p-type devices to 1L WSe2. For MoS2, we use μ0 = 81 cm2 V−1 s−1 (Ref. 48), and for WSe2, we use μ0 = 245 cm2 V−1 s−1 (Ref. 27). For this test case, we consider channel lengths LG = 100 nm, effective oxide thickness (EOT) = 2 nm, and contact resistance RC = 200 Ω⋅μm. Over the range of gate voltages applied (supply voltage VDD = 1 V), the mobility ratio of WSe2 to MoS2 is 3:1. Thus, we used a width ratio Wp/Wn = 1:3 to balance the inverter, where Wp = 1 μm and Wn = 3 μm are widths of FETs from WSe2 and MoS2, 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 (CL) to emulate the interconnect capacitance. The period of the oscillator is close to 0.5 ns for VDD = 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 (VGS0 and VBS0) using different gate metals, and optimization of device dimensions (e.g., LU) 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 MoS2 and WSe2. 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.17 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.
APPENDIX A: DERIVATION OF I–V CHARACTERISTICS
The above equation solved iteratively with Eq. (3) provides the quantum potential VC at the drain [calculated at Vn(x) = VDS] and at the source [calculated at Vn(x) = 0]. Ct and Cb are the capacitances of the top and bottom oxide, respectively. VGS and VBS are the voltages at the top-gate and the back-gate with respect to the source terminal. VGS0 and VBS0 are the flat band voltages of the top- and back-gate, respectively, for an undoped channel, and are treated as fitting parameters. Vn(x) and VC(x) are the channel potential and the quantum potential, respectively. Differentiating both sides of Eq. (A1) with respect to VC(x), we obtain
where Cq and Cit are provided in Eq. (5). We assume drift transport in the channel,
By substituting F = −dVn(x)/dx and integrating the equation over x from [0, LG], we obtain
where we have assumed the mobilty to be constant along the channel. After changing the variable from x to channel potential Vn, the above equation is re-written as
Note that Qn(Vn) [=−qn2D(Vn)] is the charge due to the free carriers (electrons in this case). VCd and VCs 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
APPENDIX B: CALCULATION OF THE CHANNEL CHARGE PARTITION
The above expressions (B1 and B2), although derived for bulk MOSFETs,44 have been previously used to calculate charge at the drain and the source nodes in carbon nanotube FETs52 as well as graphene FETs.53 The symbol θ is an empirical function of VDS, 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 − Fq such that Fq = (VDS/VDSsat)/[(VDS/VDSsat)a + 1]1/a, with a typical value a = 3 to 4, to make θ continuous. VDSsat is the saturation voltage at which the channel charge near the drain becomes zero or n2D (at x = LG) = 0. It is approximated as VDSsat = Qn0/[WLGCt(1 + ω)]. ω is a fit parameters ranging from 0 to 1.44