In this paper, we develop an analytic model based on the theory of virtual-source emission-diffusion (VS-ED) to describe ambipolar current conduction in ultrathin black phosphorus (BP) field-effect transistors (FETs). Unlike the VS model which is strictly applicable to quasiballistic devices, the VS-ED model can be applied to long-channel devices with drift-diffusive transport. The model comprehends the in-plane band structure anisotropy in BP, as well as the asymmetry in electron and hole current conduction characteristics. The model also includes the effect of Schottky-type source/drain contact resistances, which are voltage-dependent and can significantly limit current conduction in the on-state in BP FETs. Model parameters are extracted using measured data of back-gated BP transistors with gate lengths of 1000 nm and 300 nm with BP thicknesses of 7.3 nm and 8.1 nm, and for the temperature range 180–298 K. Compared to previous BP models that are validated only for room temperature and near-equilibrium bias conditions (low drain-source voltage), we demonstrate an excellent agreement between the model and data over a broad range of bias and temperature values. The model is also validated against numerical technology computer-aided design data of back- and top-gated BP transistors with a channel length of 300 nm and a thickness of 8.1 nm. The model is implemented in Verilog-A, and the capability of the model to handle both dc and transient circuit simulations is demonstrated using SPECTRE. The model not only provides physical insight into technology-device interaction in BP transistors but can also be used to design and optimize BP-based circuits using a standard hierarchical circuit simulator.

## I. INTRODUCTION

Over the past decade, two dimensional (2D) materials, including graphene,^{1–3} transition metal dichalcogenides (TMDs),^{4,5} silicene,^{6,7} and germanane,^{8} have emerged as promising candidates for future generations of nanoelectronic devices due to their excellent electrostatic integrity, vertical scalability, and electronic properties that are considerably different from those in their bulk parental materials.^{9–12} In graphene, carriers have a limited phase space for scattering, which can allow for quasiballistic transport at room temperature.^{12–14} Several high-frequency devices using graphene have been experimentally demonstrated such as radio frequency field-effect transistors (FETs), optical modulators, and photodetectors.^{9} However, the use of graphene in digital switching applications is challenging due to its low on–off current ratio, which results from its zero bandgap.^{15} On the other hand, TMD-based FETs, with the bandgap in the range of 1–2 eV, possess high on–off current ratios;^{12} however, the carrier mobility in TMDs is much lower than that in graphene.^{16}

More recently, black phosphorus (BP) has emerged as one of the most interesting 2D materials for high-performance transistor applications.^{17,18} In its bulk form, BP exhibits a bandgap of $\u2248$0.3 eV, while in its monolayer form, the bandgap increases to $\u2248$2 eV.^{9,18,19} Additionally, the hole mobility in BP is reported to be as high as 1000 cm$2$/V s at room temperature.^{18,20} Because of its puckered crystal structure, there exists a significant difference in the effective mass of carriers along the zigzag and armchair directions, with the armchair direction featuring light mass of carriers.^{9,18,19,21} With its tunable bandgap, band structure anisotropy, and high carrier mobility, BP could be an important material platform to realize digital transistors for both high-performance and low-power electronic applications.^{19,22}

Current state-of-the-art BP FETs are thick channel devices ($>$5 nm) and exhibit ambipolar current conduction with on–off current ratio around $104$ at low $Vds$, which degrades to $\u224810$ at higher $Vds$. The on–off current ratio is strongly correlated to the channel thickness, which controls the bandgap, density of states (DOSs), and gate efficiency—all of which are key parameters in determining the device performance. Calculations in Ref. 23 show that the DOS in BP is nearly independent of its thickness for energy levels close to the bottom of the conduction band, i.e., $E\u2212Ec<0.3$ eV. Here, $E$ is the electron energy and $Ec$ is the bottom of the conduction band. While the DOS increases with channel thickness at higher energy levels, its impact on maximum on-current is negligible as confirmed by numerical simulations using nonequilibrium Green’s function of a 15 nm ballistic BP FET.^{23} Thicker channel devices may also suffer from lower on–off current ratio due to ambipolarity in the device arising from significant electron injection from the drain Schottky contact into the channel. The electron injection increases with $Vds$, consequently, degrading the on–off current ratio at high $Vds$. Additionally, the gate efficiency is reduced in thicker channel devices, which further degrades the off-state characteristics of the device. For a fixed off-current, more on-current is delivered in BP FETs employing contacts that have a lower Schottky barrier with the channel. With technology optimization, including channel (thickness and length) and contact engineering, the on–off current ratio in BP FETs can be improved over a broad-bias range.

To understand, design, and simulate electronic circuits built with BP transistors, a physics-based compact model of current conduction is needed.^{24} So far, only a few models that describe current conduction in BP transistors have been reported in the literature.^{25,26} Prior BP FET models are mainly suited for near-equilibrium transport conditions, i.e., when the drain-source bias, $Vds$, is comparable to the thermal voltage, $\varphi t=kBT/q$ ($kB$ is Boltzmann’s constant, $T$ is the operating temperature, and $q$ is the elementary charge). Moreover, these models also do not include the effect of interface traps and nonlinear contact resistances, which are important to interpret experimental data. The 2D FET models reported elsewhere either introduce significant empiricism in their approach due to thermally activated and hopping-based transport^{27–29} or suffer from a limitation of not being able to reproduce the ambipolar characteristics. For example, the current–voltage (I-V) model for short-channel 2D transistors presented in Ref. 30 focuses mainly on intrinsic current conduction, while extrinsic effects due to contacts and traps are not included. The drift-diffusive I-V model presented in Ref. 31 is based on the calculation of the surface potential throughout the channel. However, the desired error-tolerance in the calculation of the surface potential must be in the range of sub-nV, which makes this model susceptible to convergence issues. Additionally, the broad-bias validity of the model in Ref. 31 requires numerical integration, which increases the computational complexity of the proposed model. As such, the model fit to experimental data has been demonstrated only for a limited bias range (gate-source bias greater than $\u22121.5V$ for transfer characteristics and drain-source bias of $\u22121.2V$ for output characteristics). Models based on the Landauer transport theory^{32,33} have also been reported in the literature.^{25,26} In Ref. 25, a Schottky barrier (SB) model to describe current conduction in off-state is developed. This model is extended in Ref. 26 to cover the on-state regime by including the channel transmission. However, this model is not suited for broad-bias circuit simulations as its validity is restricted to $Vds\u2264$ 10 mV at 300 K. We also note that none of the existing compact models for BP transistors has been implemented in Verilog-A to enable device-circuit codesign and optimization.

In this paper, we present an analytic I-V model of BP transistors to capture the ambipolar nature of current conduction over a broad range of bias voltages and temperatures. The model uses the theory of virtual-source emission-diffusion (VS-ED) to model transport in both long- and short-channel devices.^{34} Unlike the VS model, introduced previously for graphene,^{24} the VS-ED model predicts the injection velocity of carriers at high-bias based on the channel transmission coefficient and the channel potential profile. In the case of drift-diffusive transport, the channel potential profile is calculated using the gradual-channel approximation (GCA), which is valid throughout the channel length for both linear and saturation regimes. In the VS-ED model for BP FETs, charges are calculated in terms of threshold voltage, defined as a function of applied bias and ambipolar point for both electron and hole branches. To handle the nonlinear behavior of Schottky source/drain contacts, we develop a new contact resistance model, whose physics is verified through technology computer-aided design (TCAD) simulations. The model can also capture the low-temperature behavior of BP FETs by introducing temperature dependence in key model parameters, namely, the low-field mobility and threshold voltage. The model is validated by applying it to study BP FETs with gate lengths of 300 nm and 1000 nm and BP thicknesses of 7.3 nm and 8.1 nm. The model has been implemented in Verilog-A and is used to simulate the circuit behavior of BP-based inverters and ring oscillators.

The remainder of the paper is organized as follows. A description of the VS-ED model for BP FETs is presented in Sec. II. The model includes a description of charges, VS velocity, and contact resistances. Section III introduces temperature-dependent effects in BP FETs. Model validation against measured data is presented in Sec. IV. The paper concludes in Sec. V. The paper also includes three appendixes focusing on (i) parameter extraction methodology, (ii) effect of carrier concentration on mobility and I-V characteristics, and (iii) implementation of the model in Verilog-A and its use in circuit simulations.

## II. MODEL DESCRIPTION

In transistors with ambipolar current conduction, the net drain-source current, $Ids$, for a given device width, $W$, is given as^{24}

To capture the effects of carrier scattering in drift-diffusive BP FET devices with gate lengths much longer than the mean free path, we use the VS-ED model in which the empirical drain-current saturation function of the VS model is replaced with a physically derived saturation function. Moreover, the VS-ED model predicts the high-bias injection velocity rather than empirically fitting it to the experimental data.^{34–36} Per the VS-ED model for ambipolar transport, holes and electrons are in thermal equilibrium with different contacts and are, therefore, injected into the channel from different VS points. Accordingly,

Note all quantities calculated at the top-of-the-barrier or the virtual source (VS) are denoted with the subscript x$0$. In the above equation, $Qx0,e$ and $Qx0,h$ are the electron and hole charges, respectively, at the VS point in the channel (see Sec. II A for details). For both types of carriers, the injection velocity, $vinj$, at their respective VS points is expressed as a product of the thermal velocity, $vT$, and a physically derived saturation function, $Fsat$, based on emission-diffusion theory^{37} to capture the transition of current between linear and saturation regimes (see Sec. II B for the details of injection velocity). Unlike graphene FETs in which electron and hole typically have similar mobilities and velocities, BP FETs feature different physical properties for electrons and holes. As a result of this difference, $vT,e$ ($Fsat,e$) and $vT,h$ ($Fsat,h$) are treated separately in this paper.^{38}

All quantities in Eq. (2) vary with $Vgsi$ and $Vdsi$, which are the intrinsic gate-source and drain-source voltages, respectively. Intrinsic voltages are given as^{24}

where $Vgs$ and $Vds$ denote the external gate-source and drain-source voltage drops, respectively. Contact resistances corresponding to electron and hole branches are denoted as $Relec$ and $Rhole$, respectively. Section II C presents the analytic model of contact resistances.

### A. Model of VS charges

Using the threshold voltage-based model of VS charges

where $Cg$ is the gate-channel capacitance in strong inversion, $Vte$ ($Vth$) is the threshold voltage of the electron (hole) branch, $ne$ ($nh$) is the nonideality factor of the electron (hole) branch, $Vdgi=Vdsi\u2212Vgsi$. The gate-channel capacitance $Cg=\u03f5ox/CET$, where $\u03f5ox$ is the oxide permittivity and $CET$ is the capacitance equivalent thickness of the oxide given as $CET=tox+\u27e8z\u27e9\u03f5ox/\u03f5ch$. Here, $tox$ is the physical oxide thickness, $\u27e8z\u27e9$ is the finite depth of the wave-function of carriers in the channel (measured from insulator-channel interface), and $\u03f5ch$ is the dielectric permittivity of the channel. As the gate-source bias increases, the wave-function in the channel moves closer to the channel-insulator interface leading to a slight increase of $Cg$ with $Vgs$. However, for devices with thick oxides, a bias-independent $CET$, approximately equal to the physical oxide thickness, suffices.

The threshold voltage in FETs changes by a few $\varphi t$ between weak and strong inversion regimes.^{24} To accommodate this shift, the threshold voltage of electron and hole branches is modeled as

The factors $\alpha e$ and $\alpha h$ control the magnitude of threshold voltage shift between weak and strong inversion, while the functions $FFe$ and $FFh$ are Fermi functions that guarantee smooth transition between the two threshold voltage values,

The model of VS charges is validated against numerical simulations conducted in Sentaurus (for Sentaurus setup, see Sec. IV C). Figure 1 shows the net gate charge and the gate capacitance ($Cgg=\u2202Qg/\u2202Vg$) at zero drain-source bias for a representative BP FET device. We find that the model agrees well with numerical results as the device transitions from weak to strong inversion regimes. To predict the gate-channel capacitance in strong inversion accurately, the model uses $CET=16$ nm, while the physical oxide thickness $tox=15$ nm in Sentaurus. For this dataset, $\alpha e=3$ and $\alpha h=8$. The values of $\alpha e$ and $\alpha h$ are found to be in the range of 3–8 for all datasets examined in this work. These values are also similar in magnitude to the values reported in our earlier work for graphene-based FETs and bulk silicon FETs.^{24}

The threshold voltage in BP FETs shows a strong dependence on $Vds$, which for long-channel BP FETs is not due to the drain-induced barrier lowering (DIBL) but due to electron injection or leakage at the drain end.^{18,39,40} The drain leakage exacerbates at high electric fields and is particularly significant for thick channel devices that have a narrower bandgap (less than 0.7 eV for thickness greater than 6 nm).^{39} As such, the on–off current ratio degrades by over two orders of magnitude as $|Vds|$ increases (see Sec. IV on Results). The $Vds$-dependence of threshold voltage due to electron injection at the drain end and the degradation in on–off current ratio with increase in $Vds$ can be described via

where $Ve0(h0)$, $\Delta 1e(h)$, and $\Delta 2e$ are extracted from experimental calibration. The threshold voltage model uses seven parameters, out of which the parameters $Ve0(h0)$, $\Delta 1e(h)$, and $\Delta 2e$ can be extracted from the transfer curves by measuring the change in the threshold voltage with $Vdsi$ (see Appendix A for details), and the parameters $\alpha e$ and $\alpha h$ are calibrated from C-V data, if available. Otherwise, we tweak them in the range of 3–8 to match the transition between the on and off characteristics of the transistor.

### B. Velocity model

where $mi\u2217$ is the conductivity effective mass of carriers along the transport (armchair) direction, and $mDOS,i$ is the DOS effective mass of carriers given as $mDOS=mzigzag\xd7marmchair$.^{42} Here, $mzigzag$ and $marmchair$ are the effective mass of carriers along the zigzag and armchair directions, respectively. The $Fsat,i$ function in Eq. (2), describing the output-current transition from the linear regime to the saturation regime, can be physically derived by following the emission-diffusion theory of metal-semiconductor diodes outlined in Ref. 37 and applying it to FETs,^{35}

Here, $Ti=\lambda i/(\lambda i+li)$ is the channel transmission coefficient, wherein $\lambda i$ is the carrier mean free path, and $li$ is the critical back-scattering length of carriers at the VS point. In equilibrium ($Vds\u21920$), $li$ approaches the effective channel length, while it decreases with an increase in the drain-source bias. For nondegenerate carrier concentration, $\lambda i=2\mu i\varphi t/vT,i$, where $\mu i$ is the mobility of carriers.^{34} Details of mobility and its temperature dependence are presented in Sec. III A. The critical length $li$ is a function of the channel potential along the transport direction $x$ and is given as^{34}

where $L$ is the effective channel length (gated region) and $V(x)$ is the potential profile along the transport direction $x$, and $Vi(0)$ is the potential at the beginning of the channel. Once $Vi(x)$ is known, $li$ is known and, therefore, the channel transmission coefficient and the maximum injection velocity under high bias. For drift-diffusive transport, as is the case in BP FETs analyzed in our paper, $Vi(x)$ is obtained via gradual-channel approximation along $x$ and applying velocity saturation as

where $Vgt,i=Qx0i/Cg$, $Ecrit,i$ is the critical electric field at which the velocity of carriers saturates, and the parameter $\eta i$ is specified as

where $\beta i$ is an empirical parameter, typically in the range of 1.5–2.5, and is obtained upon calibration with experimental data. Equations (10)–(13) lead to the following equation for the critical length in the drift-diffusion regime:^{34}

The critical length normalized to the effective channel length is plotted in Fig. 2 as a function of $Vds$ for both electrons and holes. At all values of $Vgs$, $li\u2192L$ as $Vds\u21920$ as expected in equilibrium transport conditions. However, the behavior of critical length vs $Vgs$ is different for electrons and holes. For large negative $Vgs$ values, hole current dominates the total device current, while electron current is negligible. Hence, the ratio $lh/L$ decreases as $Vgs$ becomes negative indicating that the transport in the device is nearly unipolar p-type. On the other hand, as $Vgs$ increases and the electron current becomes significant, $le/L$ decreases, mimicking the nonequilibrium electron transport conditions in the device.

### C. Contact resistance

The parasitic resistances associated with source and drain contacts in BP FETs are bias-dependent, nonlinear resistances. This is because at the interface between the contact metal and the BP channel, a Schottky barrier is formed, which shows a strong bias dependence,^{43–45} as gate voltage changes the tunneling efficiency due to band bending at the metal and BP interface.^{46}

The current, $ISB$, through a Schottky barrier (SB) for a given voltage drop $VSB$ across it is given as^{43,47}

where $Ajun$ is the effective junction area of the SB contact, $AR$ is the Richardson constant, and $\varphi B$ is the SB height. For purely thermionic emission at the SB contact, the nonideality factor $\eta SB=1$. For other types of current conduction, such as thermally assisted tunneling and Fowler-Nordheim tunneling, $\eta SB>1$. Other factors that may lead to $\eta SB>1$ include bias dependence and image force lowering of the SB height, generation and recombination of carriers at the SB contact, and in-homogeneity of the junction.^{47} Equation (15) shows that the SB current is exponentially dependent on the SB height and the voltage drop across the barrier. The analytic model of contact resistance developed here captures these key aspects of current conduction through an SB contact. Moreover, in back-gated BP FETs we study, the region underneath the contacts is intrinsically p-doped. The application of a large negative gate voltage increases the hole doping under the contacts, which leads to a narrower barrier for hole injection and reduces the contact resistance corresponding to the hole branch.

To fully understand the effect of SB contacts on current conduction, we focus on the various transport paths in Fig. 3 that the carriers, once injected from the contacts, can take through the BP channel. In Fig. 3, path 1 for electrons and path 2 for holes are transparent for all values of $Vgs$ when $Vds>0$. For $Vds>0$, $Vgs>0$, and $Vds<Vgs$, paths labeled as $3$ and $4$ are unavailable for electron conduction. Likewise, for hole, paths labeled as $5$ and $6$ are cutoff. Hence, we can conclude that for $Vds<Vgs$, the only way carriers can be transported between the contacts is through the paths labeled as 1 and 2 in Fig. 3. Results of this discussion are summarized in Table I.

Carrier . | From . | To . | Path . | Possibility . |
---|---|---|---|---|

Electron | Source | Drain | 1 | ✓ |

Electron | Source | Drain | 3 | ✗ |

Electron | Drain | Source | 4 | ✗ |

Hole | Drain | Source | 2 | ✓ |

Hole | Source | Drain | 5 | ✗ |

Hole | Drain | Source | 6 | ✗ |

Carrier . | From . | To . | Path . | Possibility . |
---|---|---|---|---|

Electron | Source | Drain | 1 | ✓ |

Electron | Source | Drain | 3 | ✗ |

Electron | Drain | Source | 4 | ✗ |

Hole | Drain | Source | 2 | ✓ |

Hole | Source | Drain | 5 | ✗ |

Hole | Drain | Source | 6 | ✗ |

Next, we consider the case when $Vgs<Vds$, $Vds>0$, and $Vgs>0$. In this case, in addition to paths labeled as 1 and 2 in Fig. 3, there exist path $3$ for electron conduction and path $6$ for hole conduction. Results are summarized in Table II.

Carrier . | From . | To . | Path . | Possibility . |
---|---|---|---|---|

Electron | Source | Drain | 1 | ✓ |

Electron | Source | Drain | 3 | ✓ |

Electron | Drain | Source | 4 | ✗ |

Hole | Drain | Source | 2 | ✓ |

Hole | Source | Drain | 5 | ✗ |

Hole | Drain | Source | 6 | ✓ |

Carrier . | From . | To . | Path . | Possibility . |
---|---|---|---|---|

Electron | Source | Drain | 1 | ✓ |

Electron | Source | Drain | 3 | ✓ |

Electron | Drain | Source | 4 | ✗ |

Hole | Drain | Source | 2 | ✓ |

Hole | Source | Drain | 5 | ✗ |

Hole | Drain | Source | 6 | ✓ |

When $Vds<0$, $Vgs<0$, and $Vds>Vgs$, only path 2 for electrons and path 1 for holes are transparent, while for $Vds<0$, $Vgs<0$, and $Vds<Vgs$, besides paths labeled as 1 and 2 in Fig. 3, there exist path 4 for electron conduction and path 5 for hole conduction. Furthermore, when $Vds>0$ and $Vgs<0$, the possible carrier transport paths are the same as the results reported in Table I.

To further understand the nonlinearity of contact resistances, we consider the numerical simulations for back- and top-gated BP FETs using Synopsys Sentaurus. In Fig. 4, electrostatic potential distribution for $Vgs=\u22121.1$ and $Vds=\u22120.1$ V and $\u22122.5$ V is shown. Results in Figs. 4(a) and 4(b) show that the carriers paths between the source/drain contacts depend on the applied bias voltages for the back-gated device. For $Vds>Vgs$, only paths 1 and 2 labeled in Fig. 3 are transparent for current conduction. However, additional conduction paths for both carrier types contribute to current when $Vds<Vgs$. On the other hand, for top-gated structure, as shown in Figs. 4(c) and 4(d), carrier transport path becomes independent of the drain voltage, confirming the discussion in Tables I and II.

An appropriate bias-dependent model of the SB contact resistance in BP FETs must comprehend the distinct behavior for $Vgs>Vds$ and $Vgs<Vds$ regimes of transport. Note that the bias dependence of SB contact resistance is exacerbated for hole conduction. This is because of the use of Ti as the metal contact in the experimental BP FET devices examined here. Compared to metals, such as Pd or Ni, Ti has a much lower work function, which gives rise to a larger SB height and, therefore, large contact resistance values.^{10,48} Also, the effect of contact resistance corresponding to hole conduction becomes especially significant in short-channel devices in which the contact resistance could easily dominate the total drain-source resistance, limiting the maximum available current through the device.^{49} Here, we assume that $Relec$ is a bias-independent, linear resistance. This is justified because of the p-type background doping in the devices under study.

The hole branch contact resistance assuming both ohmic resistance and the formation of the SB barrier at the metal/BP interface is given as

where $a=a1FFR+a2(1\u2212FFR).$

In the above set of equations, we use the logistic function $FFR$ to model the drain-bias dependence of various parameters. This function is similar to the $FFe$/$FFh$ logistic function used in Eq. (6) and is given as

where $\gamma $ is used to adjust the sharpness of transition between two different bias regions, and $V0$ is the drain-source voltage at which additional paths for current conduction between the contacts are introduced (see Fig. 3 and the discussion). $V0$ is given as

where $V01$ and $V02$ are fitting parameters and $FFh$ is given in Eq. (6b). The model for $Rhole$ described here captures the exponential dependence on $Vgsi$ as expected for SB contacts. Moreover, the model can also explain the drain-bias dependence of the contact resistance, following the discussion pertinent to Fig. 3. The contact resistance model introduces 10 parameters: $Relec$, $Rohmic,1$, $Rohmic,2$, $R01$, $R02$, $V01$, $V02$, $\gamma $, $a1$, and $a2$. The methodology to extract the contact resistance parameters from experimental data is presented in Appendix A.

## III. EFFECT OF TEMPERATURE

The main model parameters that are affected by temperature include the mobility and threshold voltage. Below, we discuss the models to capture the temperature dependence of these parameters.

### A. Mobility

Experimentally measured mobility of carriers in 2D materials is generally much lower than that predicted theoretically based on phonon-dominated collision. Mobility degradation in 2D materials results from defects and charged impurities at room temperature.^{52–56} Previous published works indicate that the mobility of carriers in monolayer^{53,57} and multilayer^{58,59} BP in both armchair and zigzag directions^{53} follows a power law relationship with respect to temperature ($T$). That is, $\mu \u221dT\u2212\xi \mu $ for $T>100$ K^{22,53,57,58} with $\xi \mu $ typically in the range of $\u22120.4$ to 1.2 based on the type and concentration of carriers, BP crystal orientation, and the dielectric environment of the sample.^{53} In this paper, we use a constant carrier mobility model with temperature dependence given as

Here, $\mu 298,e/h$ is the value of the mobility of carriers at 298 K, $\xi \mu ,e/h$ is in the range of 0.1–0.5. The weakened temperature dependence of carrier mobility found for devices analyzed in this work could be attributed to ionized impurity-dominated scatterings. The effect of carrier concentration on mobility, as demonstrated in recent experimental work,^{59} is studied in Appendix B. However, for validation of the model against experimental data, we consider a constant carrier mobility model as in Eq. (19). This allows us to restrict the number of model parameters without compromising the quality of the model fits while also providing reasonable estimates of extracted parameters (see Sec. IV on results).

### B. Threshold voltage

The temperature dependence of the threshold voltage is mainly due to the shift in minimum conductivity point because of the frozen charge traps from the oxide and the oxide/BP interface.^{60,61} This behavior is introduced in the parameter $Ve0(h0)$ used in Eq. (7) according to

where $Ve0(h0),298$ is the value of the parameter at $T=$ 298 K, $\xi Ve0(h0)$ captures the temperature dependence of $Ve0(h0)$, and it is on the order of $10\u22123$ V/K. In thick channel devices, current at high $|Vds|$ is virtually independent of temperature implying that the tunneling current component dominates at high $|Vds|$. The temperature dependence of current at high $|Vds|$ can be reproduced by modeling $\Delta 1e(h)$ in Eq. (7) as a linear function of temperature according to

where $\Delta 1e(h),298$ is the value of $\Delta 1e(h)$ at room temperature, $\xi \Delta 1e(h)$ captures the temperature dependence of $\Delta 1e(h)$, and it is on the order of $10\u22124$ 1/K. The linear variation of $Ve(h)$ with temperature reproduces experimental and numerical results accurately as discussed in Sec. IV.

## IV. RESULTS

The analytic I-V model developed here has 29 parameters, out of which 15 (8) parameters correspond to hole (electron) conduction, and 6 parameters (3 parameters for each carrier type) model the temperature sensitivity of the current conduction. These models can be obtained through a systematic experimental validation methodology as explained in Appendix A.

There are four datasets available from experimental measurements and numerical simulation using the TCAD simulation tool Sentaurus from Synopsys.^{50} The first two datasets correspond to back-gated BP FETs with Schottky source/drain contacts, fabricated at the University of Minnesota.^{18} In these two datasets, BP flakes are exfoliated from the bulk crystal and transferred onto the local back gates. The thickness of the BP flakes forming the device are 7.3 nm and 8.1 nm. The gate dielectric is HfO$2$ with a thickness of 15 nm, and the gate metal is Ti(10 nm)/Pd(40 nm). The source and drain contacts are Ti(10 nm)/Au(90 nm). The device is passivated using 20–30 nm of Al$2$O$3$ to protect BP from atmospheric degradation. Additional fabrication details are given in Ref. 18. The third and fourth dataset are obtained numerically by simulating top- and back-gated BP FETs with Schottky source/drain contacts to show model validity for both FETs structures. All physically measured parameters or in the case of simulations, parameters selected in the simulator, are listed in Table III. Model parameters extracted for all datasets are listed in Tables IV and V.

Parameters . | Dataset 1 . | Dataset 2 . | Dataset 3 . | Dataset 4 . |
---|---|---|---|---|

Width ($\mu $m) | 7.32 | 3.16 | 3.16 | 3.16 |

Length ($\mu $m) | 1 | 0.3 | 0.3 | 0.3 |

BP thickness (nm) | 7.3 | 8.1 | 8.1 | 8.1 |

Back-oxide thickness (nm) | 15 | 15 | 15 | 15 |

Top-oxide thickness (nm) | … | … | 10 | … |

Bandgap (eV) | … | 0.69 | 0.69 | 0.69 |

Electron effective mass (unitless) | 0.15 (armchair) | 0.15 (armchair) | 0.15 (armchair) | 0.15 (armchair) |

1.18 (zigzag) | 1.18 (zigzag) | 1.18 (zigzag) | 1.18 (zigzag) | |

Hole effective mass (unitless) | 0.14 (armchair) | 0.14 (armchair) | 0.14 (armchair) | 0.14 (armchair) |

0.89 (zigzag direction) | 0.89 (zigzag) | 0.89 (zigzag) | 0.89 (zigzag) | |

BP dielectric constant (unitless) | 8.3 | 8.3 | 8.3 | 8.3 |

BP electron affinity (eV) | … | 3.655 | 3.655 | 3.655 |

Oxide dielectric constant (unitless) | 16.6 | 16.6 | 22 | 22 |

Electron critical electric field (V/m) | $13\xd7105$ | $13\xd7105$ | $13\xd7105$ | 13$\xd7105$ |

Hole critical electric field (V/m) | $14\xd7105$ | $14\xd7105$ | $14\xd7105$ | 14$\xd7105$ |

Parameters . | Dataset 1 . | Dataset 2 . | Dataset 3 . | Dataset 4 . |
---|---|---|---|---|

Width ($\mu $m) | 7.32 | 3.16 | 3.16 | 3.16 |

Length ($\mu $m) | 1 | 0.3 | 0.3 | 0.3 |

BP thickness (nm) | 7.3 | 8.1 | 8.1 | 8.1 |

Back-oxide thickness (nm) | 15 | 15 | 15 | 15 |

Top-oxide thickness (nm) | … | … | 10 | … |

Bandgap (eV) | … | 0.69 | 0.69 | 0.69 |

Electron effective mass (unitless) | 0.15 (armchair) | 0.15 (armchair) | 0.15 (armchair) | 0.15 (armchair) |

1.18 (zigzag) | 1.18 (zigzag) | 1.18 (zigzag) | 1.18 (zigzag) | |

Hole effective mass (unitless) | 0.14 (armchair) | 0.14 (armchair) | 0.14 (armchair) | 0.14 (armchair) |

0.89 (zigzag direction) | 0.89 (zigzag) | 0.89 (zigzag) | 0.89 (zigzag) | |

BP dielectric constant (unitless) | 8.3 | 8.3 | 8.3 | 8.3 |

BP electron affinity (eV) | … | 3.655 | 3.655 | 3.655 |

Oxide dielectric constant (unitless) | 16.6 | 16.6 | 22 | 22 |

Electron critical electric field (V/m) | $13\xd7105$ | $13\xd7105$ | $13\xd7105$ | 13$\xd7105$ |

Hole critical electric field (V/m) | $14\xd7105$ | $14\xd7105$ | $14\xd7105$ | 14$\xd7105$ |

. | Dataset 1 . | Dataset 2 . | ||
---|---|---|---|---|

Model parameter . | Electron . | Hole . | Electron . | Hole . |

Carrier mobility at 298 K, $\mu 298$ (cm$2$/V s) | 17.5 | 60 | 14.5 | 58 |

Temperature dependence for mobility, $\xi \mu $ (unitless) | 0.2 | 0.5 | 0.2 | 0.5 |

Nonideality factor, $n$ (unitless) | 8 | 7 | 7 | 4.8 |

$Ve(h)$ at $Vdsi=0$ V at 298 K, $Ve0(h0),298$ (V) | 0.76 | $\u22120.71$ | 1.11 | $\u22120.68$ |

Temperature dependence for $Ve0(h0)$, $\xi Ve0(h0)$ (V/K) | $10\u22123$ | $\u221210\u22123$ | $1.2\xd710\u22123$ | 0 |

Shift in ambipolar point at high $|Vds|$, $\Delta 1$ (unitless) | $\u22120.35$ | 0.55 | 0.037 | 0.18 |

Temperature dependence for $\Delta 1$, $\xi \Delta 1$ (1/K) | 0 | 0 | $5.3\xd710\u22124$ | $2.3\xd710\u22125$ |

Threshold voltage adjustment parameter at high $|Vds|$, $\Delta 2$ (1/V) | 0.1 | … | 0.14 | … |

Shift in the $Vte(h)$ in subthreshold and strong inversion, $\alpha $ (unitless) | 3 | 6 | 3 | 6 |

Empirical parameter for $Fsat$, $\beta $ (unitless) | 1.5 | 1.5 | 1.5 | 1.5 |

Contact resistance at 298 K, $Relec(hole)$ ($\Omega $ m) | 0 | … | 0 | … |

Contact ohmic resistance, $Rohmic1(2)$ ($\Omega $ m) | … | $9\xd710\u22124$(0) | … | $7\xd710\u22124$(0) |

Schottky barrier resistance coefficient, $R01(2)$ ($\Omega $ m) | … | 0.6(0.07) | … | 0.1($0.03$) |

Gate voltage dependence for Schottky barrier resistance, $a1(2)$ (1/V) | … | 2.8(1.15) | … | 2.1(1.46) |

Transition point between two different transport paths, $V01(2)$ (V) | … | 0($\u2212$0.8) | … | 0($\u2212$1) |

Adjustment the smoothness of $Rhole$, $\gamma $ (V) | … | 0.3 | … | 0.46 |

. | Dataset 1 . | Dataset 2 . | ||
---|---|---|---|---|

Model parameter . | Electron . | Hole . | Electron . | Hole . |

Carrier mobility at 298 K, $\mu 298$ (cm$2$/V s) | 17.5 | 60 | 14.5 | 58 |

Temperature dependence for mobility, $\xi \mu $ (unitless) | 0.2 | 0.5 | 0.2 | 0.5 |

Nonideality factor, $n$ (unitless) | 8 | 7 | 7 | 4.8 |

$Ve(h)$ at $Vdsi=0$ V at 298 K, $Ve0(h0),298$ (V) | 0.76 | $\u22120.71$ | 1.11 | $\u22120.68$ |

Temperature dependence for $Ve0(h0)$, $\xi Ve0(h0)$ (V/K) | $10\u22123$ | $\u221210\u22123$ | $1.2\xd710\u22123$ | 0 |

Shift in ambipolar point at high $|Vds|$, $\Delta 1$ (unitless) | $\u22120.35$ | 0.55 | 0.037 | 0.18 |

Temperature dependence for $\Delta 1$, $\xi \Delta 1$ (1/K) | 0 | 0 | $5.3\xd710\u22124$ | $2.3\xd710\u22125$ |

Threshold voltage adjustment parameter at high $|Vds|$, $\Delta 2$ (1/V) | 0.1 | … | 0.14 | … |

Shift in the $Vte(h)$ in subthreshold and strong inversion, $\alpha $ (unitless) | 3 | 6 | 3 | 6 |

Empirical parameter for $Fsat$, $\beta $ (unitless) | 1.5 | 1.5 | 1.5 | 1.5 |

Contact resistance at 298 K, $Relec(hole)$ ($\Omega $ m) | 0 | … | 0 | … |

Contact ohmic resistance, $Rohmic1(2)$ ($\Omega $ m) | … | $9\xd710\u22124$(0) | … | $7\xd710\u22124$(0) |

Schottky barrier resistance coefficient, $R01(2)$ ($\Omega $ m) | … | 0.6(0.07) | … | 0.1($0.03$) |

Gate voltage dependence for Schottky barrier resistance, $a1(2)$ (1/V) | … | 2.8(1.15) | … | 2.1(1.46) |

Transition point between two different transport paths, $V01(2)$ (V) | … | 0($\u2212$0.8) | … | 0($\u2212$1) |

Adjustment the smoothness of $Rhole$, $\gamma $ (V) | … | 0.3 | … | 0.46 |

. | Dataset 3 . | Dataset 4 . | ||
---|---|---|---|---|

Model parameter . | Electron . | Hole . | Electron . | Hole . |

Carrier mobility at 298 K, $\mu 298$ (cm$2$/V s) | 140 | 160 | … | … |

Temperature dependence for mobility, $\xi \mu $ (unitless) | 0.1 | 0.5 | … | … |

Nonideality factor, $n$ (unitless) | 3 | 2 | 2.2 | 2.3 |

$Ve(h)$ at $Vdsi=0$ V at 298 K, $Ve0(h0),298$ (V) | 0.77 | $\u22120.58$ | 0.79 | $\u22120.83$ |

Temperature dependence for $Ve0(h0)$, $\xi Ve0(h0)$ (V/K) | 0 | $\u22123.6\xd710\u22124$ | … | … |

Shift in ambipolar point at high $|Vds|$, $\Delta 1$ (unitless) | 0.1 | 0.1 | … | … |

Temperature dependence for $\Delta 1$, $\xi \Delta 1$ (1/K) | 0 | $1.5\xd710\u22123$ | … | … |

Threshold voltage adjustment parameter at high $|Vds|$, $\Delta 2$ (1/V) | 0.08 | … | … | … |

Shift in the $Vte(h)$ in subthreshold and strong inversion, $\alpha $ (unitless) | 3 | 3 | 3 | 8 |

Empirical parameter for $Fsat$, $\beta $ (unitless) | 1.5 | 1.5 | 1.5 | 1.5 |

Contact resistance at 298 K, $Relec(hole)$ ($\Omega $ m) | 0 | … | … | … |

Contact ohmic resistance, $Rohmic1(2)$ ($\Omega $ m) | … | $4\xd710\u22124$(…) | … | … |

Schottky barrier resistance coefficient, $R01(2)$ ($\Omega $ m) | … | 0.1(…) | … | … |

Gate voltage dependence for Schottky barrier resistance, $a1(2)$ (1/V) | … | 4(…) | … | … |

Transition point between two different transport paths, $V01(2)$ (V) | … | … | … | … |

Adjustment the smoothness of $Rhole$, $\gamma $ (V) | … | … | … | … |

. | Dataset 3 . | Dataset 4 . | ||
---|---|---|---|---|

Model parameter . | Electron . | Hole . | Electron . | Hole . |

Carrier mobility at 298 K, $\mu 298$ (cm$2$/V s) | 140 | 160 | … | … |

Temperature dependence for mobility, $\xi \mu $ (unitless) | 0.1 | 0.5 | … | … |

Nonideality factor, $n$ (unitless) | 3 | 2 | 2.2 | 2.3 |

$Ve(h)$ at $Vdsi=0$ V at 298 K, $Ve0(h0),298$ (V) | 0.77 | $\u22120.58$ | 0.79 | $\u22120.83$ |

Temperature dependence for $Ve0(h0)$, $\xi Ve0(h0)$ (V/K) | 0 | $\u22123.6\xd710\u22124$ | … | … |

Shift in ambipolar point at high $|Vds|$, $\Delta 1$ (unitless) | 0.1 | 0.1 | … | … |

Temperature dependence for $\Delta 1$, $\xi \Delta 1$ (1/K) | 0 | $1.5\xd710\u22123$ | … | … |

Threshold voltage adjustment parameter at high $|Vds|$, $\Delta 2$ (1/V) | 0.08 | … | … | … |

Shift in the $Vte(h)$ in subthreshold and strong inversion, $\alpha $ (unitless) | 3 | 3 | 3 | 8 |

Empirical parameter for $Fsat$, $\beta $ (unitless) | 1.5 | 1.5 | 1.5 | 1.5 |

Contact resistance at 298 K, $Relec(hole)$ ($\Omega $ m) | 0 | … | … | … |

Contact ohmic resistance, $Rohmic1(2)$ ($\Omega $ m) | … | $4\xd710\u22124$(…) | … | … |

Schottky barrier resistance coefficient, $R01(2)$ ($\Omega $ m) | … | 0.1(…) | … | … |

Gate voltage dependence for Schottky barrier resistance, $a1(2)$ (1/V) | … | 4(…) | … | … |

Transition point between two different transport paths, $V01(2)$ (V) | … | … | … | … |

Adjustment the smoothness of $Rhole$, $\gamma $ (V) | … | … | … | … |

For all datasets that are discussed in this section, we refer to the terminal that is grounded as the source terminal ($Vs=0$ V) and the drain terminal is biased at negative voltages ($Vd<0$ V). Note that this terminology of labeling the source/drain contacts is different from that followed conventionally, but it does not impact the interpretation of the model and the results. The implementation of the model in Verilog-A, SPICE circuit simulations, and Gummel symmetry test is presented in Appendix C.

### A. Dataset 1

Dataset 1 corresponds to back-gated BP FETs with $L=1$ $\mu $m, $W=7.32$ $\mu $m, $tBP=7.3$ nm, $tox,back=15$ nm, with rotational angle $42\xb0$ (rotational angle $0\xb0$ is defined along the zigzag direction and $90\xb0$ is along the armchair direction). Model fits to experimental transfer curves ($Ids$-$Vgs$) and transconductance ($gm=\u2202Ids/\u2202Vgs$) for a broad range of $Vds$ values at 298 K are shown in Figs. 5(a) and 5(b). Figures 5(c) and 5(d) show the model fits to measured output curves ($Ids$-$Vds$) and output conductance ($gd=\u2202Ids/\u2202Vds$) for a broad range of $Vgs$ values at 298 K. The model provides an excellent fit to the measured data and has the required smoothness of current derivatives as expected of compact models. We further validate this dataset using the model at $T=280$ K, 240 K, 200 K, and 180 K. Results are shown in Fig. 6.

### B. Dataset 2

The second dataset corresponds to back-gated BP FETs with $L=0.3$ $\mu $m, $W=3.16$ $\mu $m, $tBP=8.1$ nm, $tox,back=15$ nm, and rotational angle of $40\xb0$. To understand the importance of contact resistances, we plot the model results against experimental $Ids$-$Vgs$ data assuming constant $Rhole$ = 21 $\xd710\u22124$ $\Omega $ m in Fig. 7. The figure clearly shows that a constant $Rhole$ is inadequate to explain the experimental data. A higher value of $Rhole$ for $Vgs>\u22122$ V and a lower value of $Rhole$ for $Vgs<\u22122$ V can significantly improve the fit quality. Therefore, a bias-dependent nonlinear model of $Rhole$ is necessitated in this case as discussed in Sec. II C.

Model fits to the experimental I-V, transconductance, and output conductance data for dataset 2 are shown in Figs. 8(a)–8(d). Model parameters extracted from this fit are listed in Table IV. The validation of the model at various temperatures $T=280,240,200$, and 180 K is shown in Fig. 9. The model provides an excellent match with experimental data over broad bias and temperature range for this device.

Figure 10 shows the negative shift of transfer curves along the $Vgs$-axis at low temperatures, which is mainly due to the frozen charge traps at the oxide and the oxide/BP interface.^{60} This behavior is accurately modeled using Eq. (20). Figure 10 also indicates that at high $|Vds|$, current is independent of temperature. This is because the potential barrier between the drain contact and the channel becomes thinner on the drain side at high $|Vds|$, causing electron tunneling current, which is temperature insensitive, to dominate the total device current in this bias regime.

### C. Datasets 3 and 4

The last two datasets are obtained through numerical simulations using the TCAD tool Sentaurus for top-gated and back-gated BP FETs with channel length $L=0.3$ $\mu $m, $W=3.16$ $\mu $m, $tBP=8.1$ nm, top-oxide thickness, $tox,top=10$ nm, back-oxide thickness $tox,back=15$ nm, and source/drain contact length $Lcont=1.16$ $\mu $m. In Sentaurus, the Poisson and current continuity equations are solved self-consistently for both the contact and channel regions under the drift-diffusion approximation. Since Sentaurus does not provide parameters for 2D materials, we use previously published results obtained from experimental and theoretical calculations summarized in Table III.^{12,18,25} The Schottky contact is defined as a boundary condition between the contacts and the semiconductor. The work function of the metal contact is chosen as $\varphi m=4.046$ eV, and the electron affinity of BP is $\chi BP=3.655$ eV.^{18} Besides the thermionic emission, thermally assisted and direct tunneling through the SB are also considered using a nonlocal tunneling model.^{62} To include the effects of quantum confinement without solving Poisson-Schrödinger equations, we use the density-gradient method, which is an approximation of quantum confinement effects coupled with transport equations.^{63} In this simulation, we ignore the effects of traps and carrier generation and recombination.

Model fits to the numerical I-V for dataset 3 (top-gated) at various temperatures are shown in Fig. 11. Figure 1 in Sec. II A shows the net gate charge and the gate capacitance ($Cgg=\u2202Qg/\u2202Vg$) at zero drain-source bias for dataset 4. These figures show that our model agrees well with simulation results.

For all datasets, we observe that $\mu hole>\mu elec$, which is in agreement with experimental results previously reported in various BP FETs.^{12,18,42,64} The extracted values of carrier mobility, unfortunately, are notably smaller than those reported in prior studies.^{12,22,58,65} Yet, the mobility values of these devices are very close to those extracted from $gm$ measurements of similar devices fabricated at the University of Minnesota.^{40,59} The reason for low mobility values in these samples is attributed to a combination of interface scattering, remote phonon scattering from the gate dielectric and the substrate, defects, and charged impurity scattering.^{39,53} For all datasets, the on–off current ratio at all temperatures degrades by at least two orders of magnitude as $|Vds|$ increases from 0.1 V to 2.5 V. This degradation is due, in part, to significant electron injection from the drain end.

We expect that the theoretical potential of BP dictated by its finite bandgap, high carrier mobility, and anisotropic band structure will be realized through appropriate technology improvements. Specifically, thinning the BP channel to increase its bandgap and the use of low work function metals ^{16} can improve the electrostatic integrity of the device and increase its on–off current ratio-characteristics that are desirable for digital switching applications. A comprehensive assessment of the theoretical limits of BP FET technology requires quantum modeling of transport based on Monte Carlo or nonequilibrium Green’s function and is beyond the scope of this work.

## V. CONCLUSION

In this work, we develop a virtual-source emission-diffusion analytic model to describe ambipolar current conduction in drift-diffusive BP transistors over a broad range of bias and temperature values. The model comprehends the nonlinearity and the bias dependence of Schottky source/drain contacts, which is necessary to explain the I-V behavior of fabricated BP transistors. The model can also capture the low-temperature transport in BP transistors. To accomplish this, key parameters such as the carrier mobility and threshold voltage are modeled as simple temperature-dependent functions. The accuracy of the model is demonstrated by applying to top- and back-gated BP transistors with channel lengths of 300–1000 nm with temperature from 298 K to 180 K. The smoothness of the current and its derivatives is guaranteed in the model, thus satisfying a key criterion for compact models. Since the model is based on threshold voltage-based current calculation, it does not require a self-consistent solution based on surface potential. As such, the model is computationally less expensive and suitable to simulate and optimize BP-based circuits.

## ACKNOWLEDGMENTS

Elahe Yarmoghaddam and Shaloo Rakheja acknowledge the funding support of the National Science Foundation (NSF) through Grant No. CCF1565656. Shaloo Rakheja also acknowledges the funding support of NYU Wireless through the Industrial Affiliates Program. Nazila Haratipour and Steven J. Koester acknowledge primary support from the NSF through the University of Minnesota MRSEC under Award No. DMR-1420013. Portions of this work were conducted in the Minnesota Nano Center, which is supported by the NSF through the National Nanotechnology Coordinated Infrastructure under Award No. ECCS-1542202.

### APPENDIX A: EXTRACTION METHODOLOGY OF MODEL PARAMETERS

The parameter extraction begins by identifying the threshold voltage parameters. The values of $Ve0(h0)$, $\Delta 1e(h)$, and $\Delta 2e$ are determined by matching the ambipolar points from experimentally measured transfer curves and the model at all $Vds$ values. As shown in Fig. 5(a) in Sec. IV A, the ambipolar point is a strong function of $Vds$, which varies from $0.06$ V for $Vds=\u22120.1$ V to $\u22121.14$ V for $Vds=\u22122.5$ V.

For devices analyzed here, the on–off current ratio degrades by at least 100$\xd7$ as $|Vds|$ increases from 0.1 V to 2.5 V. This implies that the device is nearly always on at high $Vds$. This behavior can be explained by observing that even when the current due to hole conduction drops, the electron branch current increases, preventing the device from turning off at high $Vds$. By using an appropriate model of the electron and hole threshold voltages as in Eq. (7), we can capture the on–off device behavior in both equilibrium and off-equilibrium transport conditions. With $\Delta 2e\u22600$, we can model the decrease in $Vte$ at large negative $Vds$ values. Figure 12 shows the transfer characteristics of dataset 2 using $\Delta 2e=0$, while other parameters are equal to those listed in Table IV. The figure shows a poor match between the model and experimental data at highly negative $Vds$, implying a nonzero value is required for the parameter $\Delta 2e$ in all devices under study. The transfer curve experimental data are used to obtain the values of the nonideality factor ($n$) at low $Vds$.

We use experimental data on $gm$ measurements to extract the value of the effective mobility of carriers in the fabricated devices $(\mu =LgmWCgVds)$, where $gm$ is the peak transconductance at low $Vds$.^{39,40} For critical electric field $Ecrit,i$, we have used values reported in Ref. 51, which systematically studies saturation effects in 7–11 nm thick BP films as a function of carrier density, temperature, electric field, and the crystalline direction. Here, $Ecrit,i$ is considered to be independent of temperature since the temperature dependence of both carrier mobility and saturation velocity^{51} is rather weak. This allows us to reduce the number of model parameters without compromising the quality of model fits.

The process to find optimal parameters in $Rhole$, which is nonlinear and bias-dependent, is as follows. First, we choose a large negative value of $Vgs$ such that the term $exp\u2061(aVgsi)\u21920$ [see Eq. (16)]. This allows us to extract appropriate values of $Rohmic,1(2)$ using the measured output characteristics. On the other hand, at very low $|Vgs|$, $exp(aVgsi)$ in Eq. (16) approaches unity. Using experimental I-V data in this regime, we can identify the value of $R01(2)$. The parameters $a1(2)$, $V01(2)$, and $\gamma $ are empirical in nature and extracted by minimizing the least-square error between the model and experimental I-V data. Finally, we use the I-V measurements at different temperatures to identify the temperature dependence of key model parameters, namely, mobility and threshold voltage. The extracted temperature coefficients lie in the expected range based on previous experimental and theoretical predictions.

The empirical parameters in this model are $\beta e/h$ and $\alpha e/h$. According to previous published works for the VS model,^{24,66,67} $\beta e/h$ is tweaked in the range of 1.5–2.5, while $\alpha e/h$ is tweaked within 3–8.

Note that there is some interaction between model parameters, especially between mobility and contact resistance, and therefore, it is important to address the issue of uniqueness of model fits. The interaction between mobility and contact resistance can be decoupled (to a large extent) by noting that mobility affects the I-V curves in the linear regime of operation. We use $gm$ measurements at low bias to extract the mobility at various gate-source voltages. Once the dependence of $\mu $ on $Vgs$ (or carrier concentration) is known, we can then proceed to extracting contact resistances under high-bias conditions. Note that a change in low-field mobility does affect the mean free path and channel transmission, and therefore, the saturation current. Thus, it is important to combine the information from measurement data and theoretical results in the development of the mobility model, which we have done in Appendix B. The values of mobility reported in our work show good agreement with experimental observations.^{59} The extracted values of contact resistances in our model are also similar to those reported in the literature, with similar FET geometries and contact materials.^{39,48}

### APPENDIX B: CARRIER CONCENTRATION DEPENDENT MOBILITY MODEL

In recent experimental work, the effect of carrier concentration on both hole and electron mobility in back-gated BP FETs was explored over a broad range of temperatures.^{59} As a result of electrostatic screening in the sample, the carrier mobility was found to increase with an increase in carrier concentration for all temperatures ranging from 77 K to 295 K. To handle the variation in carrier mobility with carrier concentration, Eq. (19) is modified as

where $Be/h$ and $Q0,e/h$ are fitting parameters depending on the oxide dielectric constant, BP crystal orientation, and impurity density,^{68} $Qx0,e/h$ is also given in Eq. (4).

The updated model introduces an additional two fitting parameters, which can be determined from measurement data such as in Ref. 59. Model parameters corresponding to dataset 2 at 298 K are extracted using Eq. (B1) for the hole mobility. The fitting results are shown in Fig. 13. Here, only the parameters corresponding to the contact resistances ($Rohmic1(2)$, $R01(2)$, and $a1(2)$) are tweaked, while the remainder model parameters are the same as those reported in Table IV. With a positive value of $Be/h$, the carrier mobility increases at higher $Vgs$ (higher carrier concentration); consequently, the mean free path and channel transmission, and therefore, injection velocity and the saturation current change, necessitating a slightly larger value of the hole contact resistance to match the experimental $Ids$ data in on-state.

### APPENDIX C: VERILOG-A AND CIRCUIT SIMULATION

The model is implemented in Verilog-A to perform circuit simulations in SPICE. In Fig. 14, the model is used to simulate the dc behavior of a BP-based inverter and the transient behavior of a BP-based 3-stage ring oscillator circuit. Transistor parameters corresponding to dataset 2 are used for these simulations. The results show the mathematical robustness of the model, ease of computation, and the capability to handle large-scale circuit simulations.

The model developed in this paper also satisfies the Gummel symmetry test (GST) required for physically accurate and well-behaved compact models. According to the GST, the model must have a symmetric formulation around $Vds=0$ V. Additionally, the higher order derivatives of current must be continuous.^{69} To demonstrate that the model passes the GST, the source and drain are biased differentially ($Vx$ and $\u2212Vx$), while the gate terminal has a fixed voltage. Results of the GST reported in Fig. 14(c) verify that the model and its higher order derivatives are symmetric around $Vx=0$ V.

## REFERENCES

*et al.*, in

*2016 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2016), pp. 5.5.1–5.5.4.

*2017 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2017), pp. 31.4.1–31.4.4.

*2014 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2014), pp. 35.1.1–35.1.4.

*2017 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2017), pp. 5–7.

*et al.*,

*et al.*,

*2013 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)*(IEEE, 2013), pp. 184–187.