We present a physics-based compact model for two-dimensional (2D) field-effect transistors (FETs) based on monolayer semiconductors such as MoS_{2}. 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.

## I. INTRODUCTION

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 MoS_{2} 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_{2}, *t*_{2D} = 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 MoS_{2} 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 (*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.^{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 (*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_{2}) 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*_{0})/(*k*_{B}*T*)], *E*_{0} = *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*_{0} and the valence band energy is −*E*_{0}. *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_{2}*g*_{K} = 2, *g*_{Q} = 6,^{14}^{,}*m*_{effK} = 0.48*m*_{0}, and *m*_{effQ} = 0.57*m*_{0}.^{22} Δ*E*_{KQ} is the energy separation between K and Q conduction valleys (∼0.11 eV for monolayer MoS_{2}).^{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_{2} and ∼1.65 eV for WSe_{2}. 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

*β*= exp[−

*E*

_{it}/(

*k*

_{B}

*T*)].

*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).^{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 (*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_{2} 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_{2} 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 *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.

### B. Mobility

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 Si^{29} 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 (*T*_{0}), *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_{2} (Ref. 27), and the fit provides *η* = 6.8 and *F*_{C} = 305 V/*μ*m. Previous studies on 1L MoS_{2} have observed *η* = 1.45 and *F*_{C} = 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 MoS_{2}.^{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*_{0}*/*(1 + *N*_{OP}) where the OP (optical phonon) occupation is *N*_{OP} = 1/[exp(ℏω_{OP}/(*k*_{B}*T*)) −1]. Here, ℏω_{OP} is the OP energy and *v*_{0} can be interpreted as the saturation velocity extrapolated to zero Kelvin. For MoS_{2}, the best fit against experimental high-field data (for monolayer MoS_{2} grown by chemical vapor deposition on SiO_{2}) is obtained with ℏω_{OP} ≈ 30 to 40 meV and *v*_{0} ≈ 2 to 3 × 10^{6} cm/s.^{33} However, when modeling *I–V* curves of different devices in the literature, we must treat *v*_{0} 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 $(RTH)$ such that the average device temperature rise is Δ*T* = *T* − *T*_{0} = *P*$RTH$, where *P* = *I*_{D}(*V*_{DS} − 2*I*_{D}*R*_{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^{−7} m^{2} KW^{−1} for monolayer MoS_{2} on SiO_{2},^{39} and the “thermal area” of the device is *A* ≈ *W*(*L*_{G} + 2*L*_{U}). Due to heat spreading effects in the SiO_{2}, the effective thermal width at the SiO_{2}/Si interface can be approximated as *W*_{eff} ≈ *W* + 2*t*_{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_{2} on 90 nm SiO_{2} on Si substrate, *L*_{H} = [*k*_{2D}*t*_{2D}(*W*/*g** **+** *$RTBR$)]^{1/2} ≈ 100 nm, if we take *k*_{2D} ≈ 80 Wm^{−1} K^{−1} as the in-plane thermal conductivity of MoS_{2} at room temperature (*k*_{2D} could become a function of length in shorter devices.^{40}). For “longer” devices (with *L*_{G} + 2*L*_{U} > 3*L*_{H}), the thermal resistance is given simply by $Rth$ ≈ 1/[*g*(*L*_{G} + 2*L*_{U})]. For “shorter” channel length devices (*L*_{G} + 2*L*_{U}* **<** *3*L*_{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.^{41}

To quantify the impact of device self-heating, we simulate a typical 1L MoS_{2} 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_{2}/Si substrate (dashed), with self-heating on 300 nm SiO_{2}/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 *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.^{43} 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.^{43} 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

*ε*

_{sp}is the dielectric constant of the surrounding spacer region, and the other quantities are defined in Figs. 1 and 5(a).

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^{44}

*Q*

_{S}) and the drain (

*Q*

_{D}) are written in terms of the position dependent channel charge

*Q*

_{n}(

*x*) = −

*qn*

_{2D}(

*x*). We note, as described in Section II A, that

*n*

_{2D}(

*x*) =

*N*

_{2D}ln(1 +

*α*) where

*α*= exp[(

*qV*

_{C}(

*x*) −

*E*

_{0})/(

*k*

_{B}

*T*)], and additional details are discussed in Appendix B. Using nodal charges, we calculate the internal field capacitances between node

*m*and node

*j*as

*C*

_{if}= −∂

*Q*(

_{m}/∂V_{j}*m*

*≠*

*j*) where

*m*and

*j*are the transistor nodes (gate, source, or drain).

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^{43}

We note that when *t*_{OX} ≪ *t*_{G} ≈ *L*_{U}, then *ϕ* ≈ 1.

### E. Extrinsic resistance

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.^{7} 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*_{U}*C*_{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}.

## III. MODEL CALIBRATION

We calibrate the model to monolayer 2D semiconductors including MoS_{2} (n-type) and WSe_{2} (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^{26} and band structure. Here, we use monolayer MoS_{2} 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^{2} 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 MoS_{2} devices. We also fit the model to long-channel 1 L WSe_{2} (Ref. 27) and extract effective *μ*_{0} = 245 cm^{2} V^{−1} s^{−1}, device doping (*N*_{Dop} ∼ 10^{13} cm^{−2}), and trap density (*D*_{it} ∼ 10^{12} cm^{−2}). The model comparison to *I*_{D}–*V*_{DS} and *I*_{D}–*V*_{GS} experimental data for 1L WSe_{2} 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}, *β*, *ξ*, *V*_{GS0}, as described above. Other material parameters such as carrier effective masses,^{22} thermal conductivities,^{40} 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).

## 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 (*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_{2} and p-type devices to 1L WSe_{2}. For MoS_{2}, we use *μ*_{0} = 81 cm^{2} V^{−1} s^{−1} (Ref. 48), and for WSe_{2}, we use *μ*_{0} = 245 cm^{2} V^{−1} s^{−1} (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_{2} to MoS_{2} 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_{2} and MoS_{2}, 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.

## V. CONCLUSION

In summary, we described S2DS, a physics-based, predictive 2D transistor model capable of simulating devices and circuits from 2D semiconductors such as MoS_{2} and WSe_{2}. 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.

## ACKNOWLEDGMENTS

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

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* = −d*V*_{n}(*x*)/d*x* 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

### 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 FETs^{52} as well as graphene FETs.^{53} 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*_{G}*C*_{t}(1 + *ω*)]. *ω* is a fit parameters ranging from 0 to 1.^{44}

## References

_{2}