Large-signal model of 2DFETs: compact modeling of terminal charges and intrinsic capacitances

We present a physics-based circuit-compatible model for double-gated two-dimensional semiconductor based field effect transistors, which provides explicit expressions for the drain current, terminal charges and intrinsic capacitances. The drain current model is based on the drift-diffusion mechanism for the carrier transport and considers Fermi-Dirac statistics coupled with an appropriate field-effect approach. The terminal charge and intrinsic capacitance models are calculated adopting a Ward-Dutton linear charge partition scheme that guarantees charge-conservation. It has been implemented in Verilog-A to make it compatible with standard circuit simulators. In order to benchmark the proposed modeling framework we also present experimental DC and high-frequency measurements of a purposely fabricated monolayer MoS2 FET showing excellent agreement between the model and the experiment and thus demonstrating the capabilities of the combined approach to predict the performance of 2DFETs.


INTRODUCTION
Since the emergence of graphene, over a surprisingly short period of time, an entire new family of two-dimensional materials (2DMs) has been discovered. 1 Some of them are already used as channel materials in FETs, being promising candidates to augment silicon and III-V compound semiconductors via heterogeneous integration for advanced applications. In this context, the development of electrical models for 2D semiconductor-based FETs (2DFETs) is essential to: (1) interpret experimental results; (2) evaluate their expected digital/RF performance; (3) benchmark against existing technologies; and (4) eventually provide guidance for circuit design and circuit-level simulations.
Several compact models for three-terminal FETs based on graphene and related 2D materials (GRMs) have been recently published encompassing both static and dynamic regimes. [2][3][4][5][6][7][8][9][10] In particular, Suryavanshi et al. 2 proposed a semi-classical transport approach for the drain current combining the intrinsic FET behavior with models of the contact resistance, traps and impurities, quantum capacitance, fringing fields, high-field velocity saturation and self-heating. However, the dynamic regime is described using a charge model derived for bulk MOSFETs 11 that fails to capture the specific physics of the 2D channel, since it considers that the channel is always in a weak-inversion regime and therefore the channel charge can be assumed to be linearly distributed along the device. Wang et al. 3 reported on a compact model restricted to the static regime of graphene FETs and based on the "virtual-source" approach, valid for both the saturation and linear regions of device operation. Also based on the "virtualsource" approach, Rakheja et al. 4 studied quasi-ballistic graphene FETs. Nevertheless, they also employed the approximation derived for bulk MOSFETs 11 for the dynamic regime description even though the channel material considered is graphene. Jiménez early developed a model describing the drain current for single layer transition metal dichalcogenide (TMD) FETs based on surface potential and drift-diffusion approaches but the dynamic behavior of the TMD-FETs was not treated. 5 Also based on surface potential and drift-diffusion arguments, Parrish et al. 6 presented a compact model for graphene-based FETs for linear and non-linear circuits, but the dynamic regime description was roughly estimated by equally splitting the total gate capacitance between gate-drain and gate-source capacitances. A charge-based compact model for TMD-based FETs only valid for the static regime but simultaneously including interface traps, ambipolar current behavior, and negative capacitance was proposed by Yadav et al. 7 Taur et al. 8 presented an analytic drain current model for short-channel 2DFETs combining a subthreshold current model based on the solution to 2D Poisson's equation and a drift-diffusion approach to get an all-region static description. Rahman et al. 9 discussed a physics-based compact drain current model for monolayer TMD-FETs considering drift-diffusion transport description and the gradual channel approximation to analytically solve the Poisson's equation. Unfortunately, none of them faced the dynamic regime operation. Finally, Gholipour et al. 10 proposed a compact model that combines both the drift-diffusion and the Landauer-Buttiker approaches to properly describe the long-and short-channel FETs in flexible electronics, respectively. As a key feature, the mechanical bending is considered by linearly decreasing the channel material bandgap with respect to the applied strain. However, the capacitance model was roughly approximated by considering that the top (back) gate-drain and top (back) gatesource capacitances are just the geometrical top (back) gate capacitance. In summary, up to date the dynamic description of 2DFETs has been either roughly approximated or not specifically addressed for such devices. Our aim is to develop an accurate and physics-based compact model of the terminal charges and corresponding intrinsic capacitances of four-terminal (4T) 2DFETs valid for all operation regimes. In doing so, we present such a model together with the drain current equation, developed by some of us, 12 based on drift-diffusion theory where the surface potential is the key magnitude of the model. So that, the combination of both approaches gives rise to a compact largesignal model suitable for commercial TCAD tools employed in the simulation of circuits based on 2DFETs. The model has been shown to reproduce with good agreement DC and high-frequency measurements of a fabricated n-type MoS 2 -based FET, and has been used to predict the output transient of a three-stage ring oscillator.

Large-signal model of four-terminal 2DFETs
We first deal with the development of the large-signal model of 4T 2DFETs by the sequential formulation of the electrostatics, the drain current equation, the terminal charge model, and the corresponding intrinsic capacitance description.
The cross-section of the considered dual-gated 2DFET is depicted in Fig. 1a. For the sake of brevity, we derive here the expressions for an n-type device (the extension to a p-type transistor can be obtained straightforwardly). Provided that the semiconductor bandgap is not too small, we can assume that p(x) ≪ n(x) in all regimes of operation, where p(x) and n(x) are the hole and electron carrier densities, respectively. Upon application of 1D Gauss' law to the double-gate stack shown in Fig. 1b, the electrostatics can be described using the equivalent circuit depicted in Fig. 1c: where Q net (x) = q(p(x)n(x)) is the overall net mobile sheet charge density and q is the elementary charge.
is the top (back) oxide capacitance per unit area where ε t (ε b ) and t t (t b ) are the top (back) gate oxide relative permittivity and thickness, respectively; and V g -V g0 (V b -V b0 ) is the overdrive top (back) gate voltage. These latter quantities comprise workfunction differences between the gates and the 2D channel and any possible additional charges due to impurities or doping. The energy qV c (x) = E C (x) − E F (x) represents the shift of the quasi-Fermi level with respect to the conduction band edge and -qV(x) = E F (x) is the quasi-Fermi level along the channel. This latter quantity must fulfill the boundary conditions: (1) V(x = 0) = V s (source voltage) at the source end; (2) V(x = L) = V d (drain voltage) at the drain end. 2D crystals are far from perfect and suffer from diverse nonidealities. In particular, interface traps are unavoidable in FETs, even for those processed by state-of-the-art CMOS technology, 13 impacting negatively the performance of 2DFETs. 14 Therefore, it is mandatory to include them in order to achieve an accurate predictive TCAD tool. For n-type devices, acceptor-like traps (which are negatively charged when occupied by electrons and are energetically located in the upper half of the bandgap 15 ) contribute the most to the device electrical characteristics. 16 Assuming the traps are situated at an effective energy E it = -qV it below the conduction band, the occupied trap density can be written as n it = N it /(1 + exp(V c − V it /V th )), where N it is the effective trap density considered to be a delta function in energy. The trap charge density is then Q it = -qn it and the interface trap capacitance C it at the oxide-semiconductor interface (taken as a combination from both interfaces of the ultra-thin 2D channel) can be computed as: where V th = k B T/q is the thermal voltage, k B is the Boltzmann constant, and T is the temperature. The term C it adds to the quantum capacitance C q , as shown in Fig. 1c. We can, thus, find an expression to calculate the overall net mobile sheet density assuming a parabolic dispersion relationship modeled in the effective mass approximation and using Fermi-Dirac statistics 12 : where C dq = q 2 D 0 is defined as the degenerated-quantum capacitance, which corresponds to the upper-limit achievable when the 2D carrier density becomes heavily degenerated. 17 is the 2D density of states, with ħ being the reduced Planck's constant, g K (g Q ) the degeneracy factor and m K (m Q ) the conduction band effective mass at the K (Q) band valley. In most TMDs, the second conduction valley, labeled as Q valley, is non-negligible since the energy separation between the K and Q conduction valleys, ΔE 2 , is only around 2k B T. 18,19 Thus, two conduction band valleys may participate in the transport process. The rest of valleys are, on the contrary, far away in energy to contribute to the electrical conduction, 20 and hence can be neglected for practical purposes. The quantum capacitance of the 2DM can be computed by evaluating C q = dQ net /dV c resulting in: Due to the complexity of Eq. (1) and considering the relation between Q net and V c given by Eq. (3), it is not possible to get an explicit expression for V c as a function of the applied bias. However, an iterative Verilog-A algorithm can be implemented to evaluate the chemical potential at the source and drain edges,  which are the relevant quantities determining the drain current. In doing so, V cs = V c | V=Vs and V cd = V c | V=Vd are obtained, respectively, from Eqs. (1) and (3) by a construct 21,22 which lets the circuit simulator to iteratively solve both equations during run-time.
In the diffusive regime, the drain current of a 2DFET can be accurately computed by the following compact expression 12 : where u s = u(V cs ) and u d = u(V cd ); C tb = C t + C b ; W and L are the device width and length, respectively; and µ is the carrier mobility which has been assumed to be dependent of the applied electric field, carrier density, and temperature as in ref. 2 Once we have a suitable expression for describing the static behavior of the device, the dynamic response must be analyzed to build a largesignal model. In doing so, the terminal currents in the time domain can be computed as: where k stands for g, d, s, and b, i.e. top gate, drain, source and back-gate, respectively. The terminal currents, thus, can be obtained by either computing the time derivative of the corresponding terminal charge or by the weighted sum of the time derivative of the terminal voltages, where the weights are given by the intrinsic capacitances. The intrinsic capacitances of FETs are modeled in terms of the terminal charges. Specifically, from the electrostatics given in Eq.
(1) the following relations are derived 23 : where The charge controlled by both the drain and source terminals is computed based on the Ward-Dutton's linear charge partition scheme, 24 which guarantees charge conservation: The above expressions can be conveniently written using u as the integration variable, as it was done to model the drain current in Eq. (5). Based on the fact that the drain current is the same at any point x along the channel (i.e. we are under the quasi-static approximation framework), we get from the transport model the following equations needed to evaluate the charges in Eqs. (7)-(11): Replacing Eq. (12) in Eqs. (7) and (10), the following compact expressions are obtained for the terminal charges Q gb and Q d : where Q 2D = WLC dq V th . Equation (13) can be replaced in Eqs. (8), (9) and (11) to get the terminal charges Q g , Q b , and Q s , respectively. Q d in Eq. (13) has been simplified by assuming that the term e (2ud+2us) is prevalent over other low-order terms, e.g., e ud and e us . Once the terminal charges are evaluated, a four-terminal FET can be modeled with 4 self-capacitances and 12 intrinsic transcapacitances given by: where C ij describes the dependence of the charge at terminal i with respect to a varying voltage applied to terminal j assuming that the voltage at any other terminal remains constant. Due to charge conservation and considering a reference-independent model, only nine independent capacitances are needed to describe the 4T device. In addition, taking advantage of the relations between the top-and back-gate capacitances, 25 the computation of only four capacitances is enough; for instance, C gg , C gd , C dg , and C dd can be chosen as the independent set.

Model validation
In order to validate the presented model, we perform the DC and RF measurements of an experimental monolayer MoS 2 -FET. The fabricated device consists of an n-type channel of a chemical vapor deposited (CVD) MoS 2 monolayer transferred onto 285 nm intrinsic Si/SiO 2 wafer. The 150-nm-long MoS 2 channel is contacted with a stack of 2/70 nm Cr/Au and electrostatically controlled by an embedded gate formed by a 10-nm barrier of atomic layer deposited Al 2 O 3 and a gate metal stack consisting of 2/23 nm Ti/Au. More details on the fabrication and characterization process can be found in the Methods section. Table 1 summarizes the input parameters used for simulating the device under test. The extracted 2D semiconductor-metal contact resistance for the device is 3.5 kΩ µm and has been included in the simulation by connecting lumped resistors to the source and drain terminals.
We first test our model by comparing the DC transfer characteristics (TCs) obtained at two different drain voltages. A good agreement between both, measurements and simulations, has been achieved as shown in Fig. 2. The model predicts the DC behavior accurately in all regimes of operation with current values ranging between 10 −3 up to 10 2 µA/µm. However, the experimental device departs from the exponential trend for values in the limit of the IRDS specification for low-power applications, being the OFF current actually limited by the onset of hole current at negative V gs (see Supplementary Notes and Supplementary Fig. S1 for a detailed discussion). This small mismatch, however, does not impact the dynamic operation prediction. Indeed, we have assessed the expected RF performance of such device in Fig. 3 by benchmarking the simulated small-signal current gain (h21) and unilateral power gain (U) against the results extracted from  Fig. 4 (solid lines). The relative error between the numerical solution of Eq. (14) by using the expressions in Eqs. (8)- (11) and the analytical computation of Eq. (14) by using the compact expressions in Eq. (13) is less than 2% for the bias window considered in Fig. 4. The results are presented together with the intrinsic capacitances relying on the well-known conventional charge model developed for bulk MOSFETs 11 and elsewhere used to model 2DFETs 2 . As can be seen, the latter is accurate only either when a low drain bias is applied or if the device is operated in the subthreshold regime (see Supplementary Notes and Supplementary Fig. S2 for more details). The trends in the intrinsic capacitances can be better understood looking at the overdrive gate and drain bias dependences of both gate (Q g ) and drain (Q d ) charges, which are depicted in Fig.  5a, b, respectively. Figure 5c, d shows the same bias dependence of a set of four independent intrinsic capacitances of the 4T 2DFET under test. As the device is an n-type transistor, we have also plotted the shift of the quasi-Fermi level with respect to the conduction band edge at both the drain/source sides (corresponding to the chemical potentials V cd and V cs ) in Fig. 5e, f. For the sake of clarity, we provide a separate discussion of the intrinsic capacitances according to the transistor operation regime.
Subthreshold regime (OFF operation): If both V cd and V cs > 10V th , then Q g , Q d ∼ 0 (see Fig. 5a). Therefore the channel is empty of carriers and the device operates in the OFF state. This situation happens for overdrive gate biases lower than the threshold voltage, V A , according to Fig. 5e. Given the aforementioned conditions then C gg ≈ C t C b WL/C tb and C gd ≈ C dg ≈ C dd ≈ 0, see to the left of A1 point in Fig. 5c.
Saturation regime (ON operation): On the other hand, when V cd > 10V th while V cs < 10V th , the pinch-off is originated at the drain side. This situation is produced for overdrive gate biases between V A and V B in Fig. 5e ("A-B" section). Both C gg ("A1-B3" section) and C dg ("A1-B2" section) jump because the channel charge at the source side becomes very sensitive to the Fermi level location when it is close to the conduction band edge. However, C gd and C dd are negligible ("A1-B1" section) because the depletion region close to the drain prevents this terminal to control the channel charge. The saturation regime is also observed at the right column of Fig. 5 for drain biases higher than V C (drain saturation voltage) in Fig. 5f (right of C point). The intrinsic capacitances C gg (right of C3 point) and C dg (right of C2 point), shown in Fig. 5d, decrease with respect to the linear regime (left of C point, explained below) because the channel is depleted at the drain side, so both Q g and Q d become less sensitive to V gs as compared to the linear regime. On the other hand, C gd and C dd are negligible (right of C1 point), as Q g and Q d cannot be controlled anymore by the drain terminal after the depletion of the channel drain side, which is confirmed in Fig. 5b.
Linear (ohmic) regime (ON operation): This regime takes place when both V cd and V cs < 10V th . This situation is produced for gate overdrive voltages higher than V B in Fig. 5e (right of B point), where the channel starts to leave the depletion scenario at the drain side. Then a jump is observed in C gg (right of B3 point), C dg (right of B2 point), C gd and C dd (right of B1 point) as a consequence of the recovered electrical connection between the channel and the drain terminal (see Fig. 5a, c). The linear regime can be also observed at the right column of Fig. 5 for drain biases lower than V C (see Fig. 5f). The negative value of V cs means that the channel is degenerated at the source side for the specific bias considered in this study. In Fig. 5d, it can be seen that C gg (left of C3 point), C dg (left of C2 point), C gd and C dd (left of C1 point) decrease with V ds as the channel get closer to the drain depletion condition reached at V ds = V C .
Finally, with the aim of showing the potential of the developed TCAD tool, a relevant circuit commonly used to evaluate the performance of a digital technology, namely, a ring oscillator (RO) is simulated. The design of a three-stage based RO encompasses both DC and transient simulations. Each device has been described by the parameters shown in Table 1 that have been demonstrated to accurately reproduce the fabricated MoS 2 -based   Figs. 2 and 3. Specifically, we have designed the three-stage RO to oscillate at a specific frequency f osc = 3 1/2 /(2πR dd C) ≈ 1 GHz by using a resistor connected to each drain of R dd = 20.4 kΩ and capacitors connected to the output of each stage with C = 13.5 fF. In doing so, each single stage has been designed to show a voltage gain of 3.7 at f osc , upon biasing it at V gs = 1.25 V and V ds = 3.5 V. Figure 6 shows how the technology assessed in this work provides a suitable switching characteristic around 1 GHz. The simulation predicts a distorted sinusoidal as the AC component goes over a large range of bias with imperfect linearity. It is worth noticing that a small-signal model (in contrast with a large-signal model) could not predict that feature. Supplementary Fig. S3 provides a comparison of the compact model outcome against experimental measurements of bilayer MoS 2 FETs 26 comprising the assessment of transfer characteristics, output characteristics, the performance of an inverter, and a five-stage ring oscillator based on such devices, the latter validating the large-signal prediction capability of the model.

DISCUSSION
A physics-based large-signal compact model of four-terminal 2D semiconductor-based FETs, implemented in Verilog-A, has been presented. The model captures the terminal charges and capacitances covering all the operation regimes, so accurate time domain simulations and frequency response studies at the circuit level are feasible within the validity of the quasi-static approximation. The model can be incorporated to existing commercial circuit simulators enabling the simulation of digital and RF applications based on emerging 2D technologies. We have checked that the   6 Three-stage ring oscillator switching at a frequency around 1 GHz based on the 2DFETs described in Table 1. The supply bias is 3.5 V. model outcome is consistent with our experimental measurements of MoS 2 -FET devices targeting RF applications. Finally, a design of a three-stage ring oscillator has been carried out to exhibit the potential of the TCAD tool presented and showing the feasibility of using this technology in switching applications in the gigahertz range.

METHODS Fabrication
The MoS 2 FET reported in this paper belongs to the same batch as the devices published in ref. 27 by some of the authors. The fabrication process begins with patterning two embedded gate fingers on intrinsic Si/SiO 2 (>20 kΩ cm). Then, the embedded 150-nm gate metal stack consisting of 2/23 nm Ti/Au was defined and deposited by using electron beam lithography (EBL) and e-beam evaporation. Large area atomic single layer CVD MoS 2 , grown by a standard vapor transfer process as described in ref. , 27 was then transferred by poly(methyl methacrylate)-assisted wet transfer. Phosphoric acid etching was used to connect the embedded gate fingers and the gate pad. The active MoS 2 channel was etched using Cl 2 /O 2 plasma. Finally, source and drain (S/D) contacts consisting of 2/70 nm Cr/ Au were patterned through a final EBL step.

Electrical measurements
The electrical DC characterization was done on a Cascade Microtech Summit 11000B-AP probe-station using an Agilent B1500A parameter analyzer. Microwave performance was characterized using an Agilent twoport vector network analyzer (VNA-E8361C). All measurements were taken at room temperature, in ambient atmosphere, and in the dark. The intrinsic microwave performance is obtained after de-embedding the measured data using OPEN and SHORT structures. The OPEN de-embedding was performed on the as-measured device-under-test by etching away the MoS 2 in the active regions. The SHORT de-embedding is subsequently carried out by depositing a strip of metal across the channel region, shorting out all pads.

Circuit simulations
The developed large-signal model of 4T 2DFETs is implemented in Verilog-A and included as a separate module in Keysight© Advanced Design System (ADS), a popular electronic design automation software for RF and microwave applications. It calculates the electrostatics by iteratively solving Eqs. (1) and (3) during run-time using a construct, 21,22 and computes the static and dynamic response of the device by evaluating Eqs. (5) and (6), respectively, using the compact expressions derived in this work.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.  Fig. S1 shows the transfer characteristics (TCs) over an extended range of gate voltages. According to our simulations (solid lines), the threshold voltage of the analyzed device is Vth ~ 0.42 V. It can be observed that from the near-threshold down to the subthreshold region, the model departs from the experimental DC measurements (symbols). This is likely to happen because of the channel ambipolarity [1], [2], which has not been considered in the model. In spite of the n-type contact design, the fabricated device was not optimized to suppress hole conduction. Consequently, a residual hole current at negative Vgs appears as can be expected in intrinsic or low doped channels. It is worth noticing that the mismatch between the model and the experimental measurements is observed for current values below the IRDS limit for the OFF current in both low-power (10 -4 µA/µm) and high-performance applications (10 -1 µA/µm) (horizontal dashed lines in Fig. S1) what limits its impact in practical circuit applications and has a negligible effect in the circuit-level predictions of the model. Nevertheless, our model works properly for other devices that do not exhibit the residual hole current, as demonstrated in Refs. [3], [4].

Simplified capacitance model used to compare the proposed model
The well-known simplified capacitance model presented for bulk MOSFETs [5] and later used to model 2DFETs [4] is obtained by evaluating the partial derivatives in (S.2) using the following simplified expressions for the terminal charges:  [6] going from the linear regime η = 1 to the saturation regime η = 0. It also assumes that the net charge density is linear along the channel (Qn0), being this approach accurate only when either a low drain bias is applied or if the device is operated in the subthreshold regime [5], which in turn is in agreement with the comparison shown in FIG. 4 of the main text. In addition, Fig. S2 depicts the charge density profile along the 2D channel for different drain biases showing that it departs from linearity for higher drain biases, a fact that is captured by the model presented in this work (S.1).  Table I   To assess the validity of the proposed compact model, we have performed a variety of simulations comprising dc transfer (TCs) and output characteristics (OCs), the inverter characteristics and the performance of a five-stage ring oscillator based on depletion (D-) and enhancement (E-) mode bilayer MoS2 FETs reported in Ref. [7]. Both D-and E-mode 2DFETs effectively shift their threshold voltage by about ~0.76V due to a difference in the work function of each metal gate employed in the fabrication process. Specifically, the metals were Al and Pd [7]. Table S1 shows the input parameters used for simulating both devices, where a difference of 0.8V has been applied to the parameter Vg0 to split the FET operation into D-and E-modes. This parameter has been defined in the main text as a quantity that comprises work-function differences between the gate materials and the 2D channel and any possible additional charges due to impurities or doping. Fig. S3A shows the comparison between measurements and simulations of the TCs of both devices. Similarly to the TCs of our device (see Fig. S1), residual holes dominate the current at negative gate voltages. It is worth noticing a shift in the threshold voltage likely induced by the drain bias, which is not captured by our model. However, that effect is not apparent in the device analyzed in the main manuscript. Fig. S3B shows the

A B C D
comparison of the OCs of the E-mode bilayer MoS2 FET. In the saturation region, a mismatch between the model prediction and the dc measurements is detected which we attribute to a channel-length modulation effect [5], not considered in our model. In spite of these small mismatches the model captures to a very good degree large-signal operation. In particular, Fig. S3C shows the output voltage as a function of the input voltage in a bilayer MoS2 based inverter formed by the series combination of a D-mode FET in resistor configuration and an E-mode FET (see inset of Fig. S3C). A voltage gain close to 5 is achieved and used to build a five-stage ring oscillator. In doing so, a loop is constructed by connecting five inverters in series. Fig. S3D shows the output voltage switching at a fundamental frequency of ~1.6 MHz corresponding to a propagation delay per stage of 62.5 ns.