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.

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 0.3 eV, while in its monolayer form, the bandgap increases to 2 eV.9,18,19 Additionally, the hole mobility in BP is reported to be as high as 1000 cm2/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 10 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., EEc<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, ϕ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 transport27–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 1.5V for transfer characteristics and drain-source bias of 1.2V for output characteristics). Models based on the Landauer transport theory32,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 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.

In transistors with ambipolar current conduction, the net drain-source current, Ids, for a given device width, W, is given as24 

Ids=Ielec+Ihole.
(1)

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,

Ielec=WQx0,evT,eFsat,evinj,e,
(2a)
Ihole=WQx0,hvT,hFsat,hvinj,h.
(2b)

Note all quantities calculated at the top-of-the-barrier or the virtual source (VS) are denoted with the subscript x0. 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 theory37 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 as24 

Vgsi=Vgs(IelecRelec+IholeRhole),
(3a)
Vdsi=Vds2(IelecRelec+IholeRhole),
(3b)

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.

Using the threshold voltage-based model of VS charges

Qx0e=Cgneϕtlog(1+exp(VgsiVteneϕt)),
(4a)
Qx0h=Cgnhϕtlog(1+exp(Vdgi+Vthnhϕt)),
(4b)

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=VdsiVgsi. The gate-channel capacitance Cg=ϵox/CET, where ϵox is the oxide permittivity and CET is the capacitance equivalent thickness of the oxide given as CET=tox+zϵox/ϵch. Here, tox is the physical oxide thickness, z is the finite depth of the wave-function of carriers in the channel (measured from insulator-channel interface), and ϵch 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 ϕt between weak and strong inversion regimes.24 To accommodate this shift, the threshold voltage of electron and hole branches is modeled as

Vte=VeαeϕtFFe,
(5a)
Vth=Vh+αhϕtFFh.
(5b)

The factors αe and α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,

FFe=11+exp(Vgsi(Veαeϕt2)αeϕt),
(6a)
FFh=11+exp(Vdgi+(Vh+αhϕt2)αhϕt).
(6b)

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=Qg/Vg) 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, αe=3 and αh=8. The values of αe and α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 

FIG. 1.

Model fit to the numerical data obtained from TCAD simulation at room temperature for Vds=0 V for BP FET with channel length L=0.3μm, width W=3.16μm, BP thickness tBP=8.1 nm, and oxide thickness tox=15 nm: (a) net gate charge (Qg=Qx0,eQx0,h) as a function of Vgs and (b) gate capacitance (Cgg=Qg/Vgs) as a function of Vgs. Solid lines are model fits, while symbols correspond to numerical data. For model, CET=16 nm, ne(h)=2.2(2.3), Ve0(h0)=0.79(0.83) V, and αe(h)=3(8).

FIG. 1.

Model fit to the numerical data obtained from TCAD simulation at room temperature for Vds=0 V for BP FET with channel length L=0.3μm, width W=3.16μm, BP thickness tBP=8.1 nm, and oxide thickness tox=15 nm: (a) net gate charge (Qg=Qx0,eQx0,h) as a function of Vgs and (b) gate capacitance (Cgg=Qg/Vgs) as a function of Vgs. Solid lines are model fits, while symbols correspond to numerical data. For model, CET=16 nm, ne(h)=2.2(2.3), Ve0(h0)=0.79(0.83) V, and αe(h)=3(8).

Close modal

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

Ve=Ve0Δ1eVdsiΔ2eVdsi2,
(7a)
Vh=Vh0Δ1hVdsi,
(7b)

where Ve0(h0), Δ1e(h), and Δ2e are extracted from experimental calibration. The threshold voltage model uses seven parameters, out of which the parameters Ve0(h0), Δ1e(h), and Δ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 αe and α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.

The thermal velocity vT,i (i=e/h for electrons/holes) in Eq. (2) is given as34,41

vT,i=2kBTmi2πmDOS,i2,
(8)

where mi 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×marmchair.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 

Fsat,i=(Ti2Ti)(1exp(Vdsiϕt)1+Ti2Tiexp(Vdsiϕt)).
(9)

Here, Ti=λi/(λi+li) is the channel transmission coefficient, wherein λi is the carrier mean free path, and li is the critical back-scattering length of carriers at the VS point. In equilibrium (Vds0), li approaches the effective channel length, while it decreases with an increase in the drain-source bias. For nondegenerate carrier concentration, λi=2μiϕt/vT,i, where μ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 as34 

li=0Lexp(Vi(x)Vi(0)ϕt)dx,
(10)

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

Vi(x)=Vgt1,i(11xL(1ηi2)),
(11)
Vgt1,i=2Vgt,i1+2Vgt,iLEcrit,i,
(12)

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

ηi=1VdsiVgt1,i(1+(VdsiVgt1,i)βi)1/βi,
(13)

where β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 

li=2L(Vgt1,iϕt)2(1ηi2)(exp(Vgt1,iϕt(ηi1))×(1Vgt1,iϕtηi)(1Vgt1,iϕt)).
(14)

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, liL as Vds0 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.

FIG. 2.

(a) Hole and (b) electron critical lengths normalized to channel length as a function of Vds for Vgs=2.5 V (solid line), Vgs=1 V (dashed line), and Vgs=0.5 V (dashed-dotted line).

FIG. 2.

(a) Hole and (b) electron critical lengths normalized to channel length as a function of Vds for Vgs=2.5 V (solid line), Vgs=1 V (dashed line), and Vgs=0.5 V (dashed-dotted line).

Close modal

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 as43,47

ISB=AjunART2exp(ϕBϕt)[exp(VSBηSBϕt)1],
(15)

where Ajun is the effective junction area of the SB contact, AR is the Richardson constant, and ϕB is the SB height. For purely thermionic emission at the SB contact, the nonideality factor ηSB=1. For other types of current conduction, such as thermally assisted tunneling and Fowler-Nordheim tunneling, ηSB>1. Other factors that may lead to η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.

FIG. 3.

Possible carrier transport paths in the channel. Note that the figure is not drawn to scale.

FIG. 3.

Possible carrier transport paths in the channel. Note that the figure is not drawn to scale.

Close modal
TABLE I.

Possibility of different carrier transport paths in a Schottky barrier back-gated MOSFET for Vds>0, Vgs>0, and Vds<Vgs.

CarrierFromToPathPossibility
Electron Source Drain ✓ 
Electron Source Drain ✗ 
Electron Drain Source ✗ 
Hole Drain Source ✓ 
Hole Source Drain ✗ 
Hole Drain Source ✗ 
CarrierFromToPathPossibility
Electron Source Drain ✓ 
Electron Source Drain ✗ 
Electron Drain Source ✗ 
Hole Drain Source ✓ 
Hole Source Drain ✗ 
Hole Drain Source ✗ 

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.

TABLE II.

Possibility of different carrier transport paths in a Schottky barrier back-gated MOSFET for Vds>0, Vgs>0, and Vds>Vgs.

CarrierFromToPathPossibility
Electron Source Drain ✓ 
Electron Source Drain ✓ 
Electron Drain Source ✗ 
Hole Drain Source ✓ 
Hole Source Drain ✗ 
Hole Drain Source ✓ 
CarrierFromToPathPossibility
Electron Source Drain ✓ 
Electron Source Drain ✓ 
Electron Drain Source ✗ 
Hole Drain Source ✓ 
Hole Source Drain ✗ 
Hole Drain Source ✓ 

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=1.1 and Vds=0.1 V and 2.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.

FIG. 4.

Surface plot of the electrostatic potential for the back-gated simulated device with L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, and Lcont=1.16μm for (a) Vgs=1.1 V and Vds=0.1 V and (b) Vgs=1.1 V and Vds=2.5 V, and top-gated structure with L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, tox,top=10 nm, and Lcont=1.16μm for (c) Vgs=1.1 V and Vds=0.1 V and (d) Vgs=1.1 V and Vds=2.5 V. Solid and dashed lines represent the flow of electrons and holes, respectively, through the channel.

FIG. 4.

Surface plot of the electrostatic potential for the back-gated simulated device with L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, and Lcont=1.16μm for (a) Vgs=1.1 V and Vds=0.1 V and (b) Vgs=1.1 V and Vds=2.5 V, and top-gated structure with L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, tox,top=10 nm, and Lcont=1.16μm for (c) Vgs=1.1 V and Vds=0.1 V and (d) Vgs=1.1 V and Vds=2.5 V. Solid and dashed lines represent the flow of electrons and holes, respectively, through the channel.

Close modal

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

Rhole=Rohmic,1FFR+Rohmic,2(1FFR)Ohmiccontactresistance+(R01FFR+R02(1FFR))exp(aVgsi)Schottkybarrierresistance,
(16)

where a=a1FFR+a2(1FFR).

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

FFR=11+exp(Vdsi+V0γ),
(17)

where γ 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

V0=V01FFh+V02(1FFh),
(18)

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, γ, a1, and a2. The methodology to extract the contact resistance parameters from experimental data is presented in  Appendix A.

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.

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 monolayer53,57 and multilayer58,59 BP in both armchair and zigzag directions53 follows a power law relationship with respect to temperature (T). That is, μTξμ for T>100 K22,53,57,58 with ξμ typically in the range of 0.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

μe/h=μ298,e/h(298T)ξμ,e/h.
(19)

Here, μ298,e/h is the value of the mobility of carriers at 298 K, ξμ,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).

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

Ve0(h0)=Ve0(h0),298ξVe0(h0)(T298),
(20)

where Ve0(h0),298 is the value of the parameter at T= 298 K, ξVe0(h0) captures the temperature dependence of Ve0(h0), and it is on the order of 103 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 Δ1e(h) in Eq. (7) as a linear function of temperature according to

Δ1e(h)=Δ1e(h),298ξΔ1e(h)(T298),
(21)

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

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 HfO2 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 Al2O3 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.

TABLE III.

Physical parameters of various devices analyzed in this work.12,18,50,51 The value of bandgap and BP electron affinity are not reported for dataset 1.

ParametersDataset 1Dataset 2Dataset 3Dataset 4
Width (μm) 7.32 3.16 3.16 3.16 
Length (μm) 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×105 13×105 13×105 13×105 
Hole critical electric field (V/m) 14×105 14×105 14×105 14×105 
ParametersDataset 1Dataset 2Dataset 3Dataset 4
Width (μm) 7.32 3.16 3.16 3.16 
Length (μm) 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×105 13×105 13×105 13×105 
Hole critical electric field (V/m) 14×105 14×105 14×105 14×105 
TABLE IV.

Model parameters for experimental dataset 1 (L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm) and dataset 2 (L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm) obtained from Ref. 18.

Dataset 1Dataset 2
Model parameterElectronHoleElectronHole
Carrier mobility at 298 K, μ298 (cm2/V s) 17.5 60 14.5 58 
Temperature dependence for mobility, ξμ (unitless) 0.2 0.5 0.2 0.5 
Nonideality factor, n (unitless) 4.8 
Ve(h) at Vdsi=0 V at 298 K, Ve0(h0),298 (V) 0.76 0.71 1.11 0.68 
Temperature dependence for Ve0(h0), ξVe0(h0) (V/K) 103 103 1.2×103 
Shift in ambipolar point at high |Vds|, Δ1 (unitless) 0.35 0.55 0.037 0.18 
Temperature dependence for Δ1, ξΔ1 (1/K) 5.3×104 2.3×105 
Threshold voltage adjustment parameter at high |Vds|, Δ2 (1/V) 0.1 … 0.14 … 
Shift in the Vte(h) in subthreshold and strong inversion, α (unitless) 
Empirical parameter for Fsat, β (unitless) 1.5 1.5 1.5 1.5 
Contact resistance at 298 K, Relec(hole) (Ω m) … … 
Contact ohmic resistance, Rohmic1(2) (Ω m) … 9×104(0) … 7×104(0) 
Schottky barrier resistance coefficient, R01(2) (Ω 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(0.8) … 0(1) 
Adjustment the smoothness of Rhole, γ (V) … 0.3 … 0.46 
Dataset 1Dataset 2
Model parameterElectronHoleElectronHole
Carrier mobility at 298 K, μ298 (cm2/V s) 17.5 60 14.5 58 
Temperature dependence for mobility, ξμ (unitless) 0.2 0.5 0.2 0.5 
Nonideality factor, n (unitless) 4.8 
Ve(h) at Vdsi=0 V at 298 K, Ve0(h0),298 (V) 0.76 0.71 1.11 0.68 
Temperature dependence for Ve0(h0), ξVe0(h0) (V/K) 103 103 1.2×103 
Shift in ambipolar point at high |Vds|, Δ1 (unitless) 0.35 0.55 0.037 0.18 
Temperature dependence for Δ1, ξΔ1 (1/K) 5.3×104 2.3×105 
Threshold voltage adjustment parameter at high |Vds|, Δ2 (1/V) 0.1 … 0.14 … 
Shift in the Vte(h) in subthreshold and strong inversion, α (unitless) 
Empirical parameter for Fsat, β (unitless) 1.5 1.5 1.5 1.5 
Contact resistance at 298 K, Relec(hole) (Ω m) … … 
Contact ohmic resistance, Rohmic1(2) (Ω m) … 9×104(0) … 7×104(0) 
Schottky barrier resistance coefficient, R01(2) (Ω 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(0.8) … 0(1) 
Adjustment the smoothness of Rhole, γ (V) … 0.3 … 0.46 
TABLE V.

Model parameters for numerical datasets 3 and 4. L=0.3μm, W=3.16μm, tBP=8.1 nm, back-oxide thickness tox,back=15 nm, and source/drain contact length, Lcont=1.16μm. Top-oxide thickness for dataset 3 is tox,top=10 nm.

Dataset 3Dataset 4
Model parameterElectronHoleElectronHole
Carrier mobility at 298 K, μ298 (cm2/V s) 140 160 … … 
Temperature dependence for mobility, ξμ (unitless) 0.1 0.5 … … 
Nonideality factor, n (unitless) 2.2 2.3 
Ve(h) at Vdsi=0 V at 298 K, Ve0(h0),298 (V) 0.77 0.58 0.79 0.83 
Temperature dependence for Ve0(h0), ξVe0(h0) (V/K) 3.6×104 … … 
Shift in ambipolar point at high |Vds|, Δ1 (unitless) 0.1 0.1 … … 
Temperature dependence for Δ1, ξΔ1 (1/K) 1.5×103 … … 
Threshold voltage adjustment parameter at high |Vds|, Δ2 (1/V) 0.08 … … … 
Shift in the Vte(h) in subthreshold and strong inversion, α (unitless) 
Empirical parameter for Fsat, β (unitless) 1.5 1.5 1.5 1.5 
Contact resistance at 298 K, Relec(hole) (Ω m) … … … 
Contact ohmic resistance, Rohmic1(2) (Ω m) … 4×104(…) … … 
Schottky barrier resistance coefficient, R01(2) (Ω 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, γ (V) … … … … 
Dataset 3Dataset 4
Model parameterElectronHoleElectronHole
Carrier mobility at 298 K, μ298 (cm2/V s) 140 160 … … 
Temperature dependence for mobility, ξμ (unitless) 0.1 0.5 … … 
Nonideality factor, n (unitless) 2.2 2.3 
Ve(h) at Vdsi=0 V at 298 K, Ve0(h0),298 (V) 0.77 0.58 0.79 0.83 
Temperature dependence for Ve0(h0), ξVe0(h0) (V/K) 3.6×104 … … 
Shift in ambipolar point at high |Vds|, Δ1 (unitless) 0.1 0.1 … … 
Temperature dependence for Δ1, ξΔ1 (1/K) 1.5×103 … … 
Threshold voltage adjustment parameter at high |Vds|, Δ2 (1/V) 0.08 … … … 
Shift in the Vte(h) in subthreshold and strong inversion, α (unitless) 
Empirical parameter for Fsat, β (unitless) 1.5 1.5 1.5 1.5 
Contact resistance at 298 K, Relec(hole) (Ω m) … … … 
Contact ohmic resistance, Rohmic1(2) (Ω m) … 4×104(…) … … 
Schottky barrier resistance coefficient, R01(2) (Ω 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, γ (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.

Dataset 1 corresponds to back-gated BP FETs with L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm, with rotational angle 42° (rotational angle 0° is defined along the zigzag direction and 90° is along the armchair direction). Model fits to experimental transfer curves (Ids-Vgs) and transconductance (gm=Ids/Vgs) 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=Ids/Vds) 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.

FIG. 5.

Model fit to the experimental dataset 1 at room temperature: (a) transfer characteristics (Ids-Vgs), (b) transconductance (gm=Ids/Vgs), (c) output characteristics (Ids-Vds), and (d) output conductance (gd=Ids/Vds). Solid lines are model fits, while symbols correspond to experimental data. Dataset 1: L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm, and T=298 K.

FIG. 5.

Model fit to the experimental dataset 1 at room temperature: (a) transfer characteristics (Ids-Vgs), (b) transconductance (gm=Ids/Vgs), (c) output characteristics (Ids-Vds), and (d) output conductance (gd=Ids/Vds). Solid lines are model fits, while symbols correspond to experimental data. Dataset 1: L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm, and T=298 K.

Close modal
FIG. 6.

Model fit to the experimental dataset 1 at T= (a) 280 K, (b) 240 K, (c) 200 K, and (d) 180 K for Vds:0.1,0.2,1.5,2, and 2.5 V. Solid lines are model fits, while symbols correspond to experimental data. Dataset 1: L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm at different temperatures.

FIG. 6.

Model fit to the experimental dataset 1 at T= (a) 280 K, (b) 240 K, (c) 200 K, and (d) 180 K for Vds:0.1,0.2,1.5,2, and 2.5 V. Solid lines are model fits, while symbols correspond to experimental data. Dataset 1: L=1μm, W=7.32μm, tBP=7.3 nm, tox,back=15 nm at different temperatures.

Close modal

The second dataset corresponds to back-gated BP FETs with L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, and rotational angle of 40°. To understand the importance of contact resistances, we plot the model results against experimental Ids-Vgs data assuming constant Rhole = 21 ×104Ω 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>2 V and a lower value of Rhole for Vgs<2 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.

FIG. 7.

Transfer characteristics (Ids-Vgs) of dataset 2 with constant Rhole and Relec at room temperature. Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

FIG. 7.

Transfer characteristics (Ids-Vgs) of dataset 2 with constant Rhole and Relec at room temperature. Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

Close modal

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.

FIG. 8.

Model fit to the experimental dataset 2 at room temperature: (a) transfer characteristics (Ids-Vgs), (b) transconductance (gm=Ids/Vgs), (c) output characteristics (Ids-Vds), and (d) output conductance (gd=Ids/Vds). Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, and T=298 K.

FIG. 8.

Model fit to the experimental dataset 2 at room temperature: (a) transfer characteristics (Ids-Vgs), (b) transconductance (gm=Ids/Vgs), (c) output characteristics (Ids-Vds), and (d) output conductance (gd=Ids/Vds). Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, and T=298 K.

Close modal
FIG. 9.

Model fit to the experimental dataset 2 at T= (a) 280 K, (b) 240 K, (c) 200 K, and (d) 180 K for Vds:0.1,1,1.5,2, and 2.5 V. Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm at different temperatures.

FIG. 9.

Model fit to the experimental dataset 2 at T= (a) 280 K, (b) 240 K, (c) 200 K, and (d) 180 K for Vds:0.1,1,1.5,2, and 2.5 V. Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm at different temperatures.

Close modal

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.

FIG. 10.

Experimental data for transfer characteristics (Ids-Vgs) of dataset 2 for Vds=0.1,2.5 V at different temperatures. Solid line, dotted line, dashed-dotted line, and dashed line correspond to the T=298 K, 240 K, 200 K, and 180 K, respectively. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm at different temperatures.

FIG. 10.

Experimental data for transfer characteristics (Ids-Vgs) of dataset 2 for Vds=0.1,2.5 V at different temperatures. Solid line, dotted line, dashed-dotted line, and dashed line correspond to the T=298 K, 240 K, 200 K, and 180 K, respectively. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm at different temperatures.

Close modal

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μm, W=3.16μ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μ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 ϕm=4.046 eV, and the electron affinity of BP is χ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=Qg/Vg) at zero drain-source bias for dataset 4. These figures show that our model agrees well with simulation results.

FIG. 11.

Model fit to the numerical dataset 3: (a) output characteristics (Ids-Vds) at room temperature, and transfer characteristics (Ids-Vgs) at T= (b) 298 K, (c) 240 K, and (d) 200 K. Solid lines are model fits, while symbols correspond to numerical data. Dataset 3: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, tox,top=10 nm, Lcont=1.16μm.

FIG. 11.

Model fit to the numerical dataset 3: (a) output characteristics (Ids-Vds) at room temperature, and transfer characteristics (Ids-Vgs) at T= (b) 298 K, (c) 240 K, and (d) 200 K. Solid lines are model fits, while symbols correspond to numerical data. Dataset 3: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm, tox,top=10 nm, Lcont=1.16μm.

Close modal

For all datasets, we observe that μhole>μ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.

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.

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.

The parameter extraction begins by identifying the threshold voltage parameters. The values of Ve0(h0), Δ1e(h), and Δ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=0.1 V to 1.14 V for Vds=2.5 V.

For devices analyzed here, the on–off current ratio degrades by at least 100× 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 Δ2e0, we can model the decrease in Vte at large negative Vds values. Figure 12 shows the transfer characteristics of dataset 2 using Δ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 Δ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.

FIG. 12.

Transfer characteristics (Ids-Vgs) for dataset 2, at room temperature, by setting Δ2e=0 in Eq. (7). Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

FIG. 12.

Transfer characteristics (Ids-Vgs) for dataset 2, at room temperature, by setting Δ2e=0 in Eq. (7). Solid lines are model fits, while symbols correspond to experimental data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

Close modal

We use experimental data on gm measurements to extract the value of the effective mobility of carriers in the fabricated devices (μ=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 velocity51 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(aVgsi)0 [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 γ 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 βe/h and αe/h. According to previous published works for the VS model,24,66,67βe/h is tweaked in the range of 1.5–2.5, while α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 μ 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

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

μe/h=μ298,e/h(1+Qx0,e/hQ0,e/h)Be/h(298T)ξμ,e/h,
(B1)

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.

FIG. 13.

Transfer characteristics (Ids-Vgs) for dataset 2, at room temperature, by using Eq. (B1) for hole mobility. μ0,h=58 cm2/V s, Bh=1.5, Q0,h=0.03 C/m2, Rohmic1(2)=0.0017(0.0017)Ω m, R01(2)=0.11(0.017)Ω m, a1(2)=2.2(1.2) 1/V, all other fitting parameters are equal to the values mentioned in Table IV. Solid lines are model fits, while symbols correspond to numerical data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

FIG. 13.

Transfer characteristics (Ids-Vgs) for dataset 2, at room temperature, by using Eq. (B1) for hole mobility. μ0,h=58 cm2/V s, Bh=1.5, Q0,h=0.03 C/m2, Rohmic1(2)=0.0017(0.0017)Ω m, R01(2)=0.11(0.017)Ω m, a1(2)=2.2(1.2) 1/V, all other fitting parameters are equal to the values mentioned in Table IV. Solid lines are model fits, while symbols correspond to numerical data. Dataset 2: L=0.3μm, W=3.16μm, tBP=8.1 nm, tox,back=15 nm.

Close modal

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.

FIG. 14.

(a) DC simulation of inverter: output voltage as a function of input voltage for VDD=1 V. The inset shows schematic of an inverter, (b) transient simulation for 3-stage ring oscillator where CL=3 fF. The inset indicates the schematic of the 3-stage ring oscillator (transistor parameters are the same as dataset 2), (c) results of the Gummel symmetry test (GST) for the (i) current, (ii) first, (iii) second, and (iv) third order derivatives of current with respect to Vx for Vgs=1 V with arbitrary units.

FIG. 14.

(a) DC simulation of inverter: output voltage as a function of input voltage for VDD=1 V. The inset shows schematic of an inverter, (b) transient simulation for 3-stage ring oscillator where CL=3 fF. The inset indicates the schematic of the 3-stage ring oscillator (transistor parameters are the same as dataset 2), (c) results of the Gummel symmetry test (GST) for the (i) current, (ii) first, (iii) second, and (iv) third order derivatives of current with respect to Vx for Vgs=1 V with arbitrary units.

Close modal

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 Vx), 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.

1.
E.
Yarmoghaddam
and
S.
Rakheja
,
J. Appl. Phys.
122
,
083101
(
2017
).
2.
S.
Rakheja
,
V.
Kumar
, and
A.
Naeemi
,
Proc. IEEE
101
,
1740
(
2013
).
3.
K. S.
Novoselov
,
A. K.
Geim
,
S.
Morozov
,
D.
Jiang
,
M.
Katsnelson
,
I.
Grigorieva
,
S.
Dubonos
, and
A. A.
Firsov
,
Nature
438
,
197
(
2005
).
4.
Q. H.
Wang
,
K.
Kalantar-Zadeh
,
A.
Kis
,
J. N.
Coleman
, and
M. S.
Strano
,
Nat. Nanotechnol.
7
,
699
(
2012
).
5.
R.
Fivaz
and
E.
Mooser
,
Phys. Rev.
163
,
743
(
1967
).
6.
M.
Houssa
,
E.
Scalise
,
K.
Sankaran
,
G.
Pourtois
,
V.
Afanas’ev
, and
A.
Stesmans
,
Appl. Phys. Lett.
98
,
223107
(
2011
).
7.
P.
Vogt
,
P.
De Padova
,
C.
Quaresima
,
J.
Avila
,
E.
Frantzeskakis
,
M. C.
Asensio
,
A.
Resta
,
B.
Ealet
, and
G.
Le Lay
,
Phys. Rev. Lett.
108
,
155501
(
2012
).
8.
E.
Bianco
,
S.
Butler
,
S.
Jiang
,
O. D.
Restrepo
,
W.
Windl
, and
J. E.
Goldberger
,
ACS Nano
7
,
4414
(
2013
).
9.
F.
Xia
,
H.
Wang
,
D.
Xiao
,
M.
Dubey
, and
A.
Ramasubramaniam
,
Nat. Photonics
8
,
899
(
2014
).
10.
S.
Das
,
M.
Demarteau
, and
A.
Roelofs
,
ACS Nano
8
,
11730
(
2014
).
11.
Y.
Yi
,
X.-F.
Yu
,
W.
Zhou
,
J.
Wang
, and
P. K.
Chu
,
Mater. Sci. Eng. R Rep.
120
,
1
(
2017
).
12.
J.
Qiao
,
X.
Kong
,
Z.-X.
Hu
,
F.
Yang
, and
W.
Ji
,
Nat. Commun.
5
,
4475
(
2014
).
13.
A. K.
Geim
and
K. S.
Novoselov
,
Nat. Mater.
6
,
183
(
2007
).
14.
E.
Hwang
and
S. D.
Sarma
,
Phys. Rev. B
77
,
115449
(
2008
).
15.
X.
Han
,
H. M.
Stewart
,
S. A.
Shevlin
,
C. R. A.
Catlow
, and
Z. X.
Guo
,
Nano Lett.
14
,
4607
(
2014
).
16.
Y.
Du
,
H.
Liu
,
Y.
Deng
, and
P. D.
Ye
,
ACS Nano
8
,
10035
(
2014
).
17.
H.
Huang
,
Q.
Xiao
,
J.
Wang
,
X.-F.
Yu
,
H.
Wang
,
H.
Zhang
, and
P. K.
Chu
,
NPJ 2D Mater. Appl.
1
,
20
(
2017
).
18.
N.
Haratipour
,
S.
Namgung
,
S.-H.
Oh
, and
S. J.
Koester
,
ACS Nano
10
,
3791
(
2016
).
19.
L.
Yang
,
G.
Qiu
,
M.
Si
,
A.
Charnas
,
C.
Milligan
,
D.
Zemlyanov
,
H.
Zhou
,
Y.
Du
,
Y.
Lin
, and
W.
Tsai
et al., in 2016 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2016), pp. 5.5.1–5.5.4.
20.
S.
Bohloul
,
L.
Zhang
,
K.
Gong
, and
H.
Guo
,
Appl. Phys. Lett.
108
,
033508
(
2016
).
21.
Z.
Luo
,
J.
Maassen
,
Y.
Deng
,
Y.
Du
,
R. P.
Garrelts
,
M. S.
Lundstrom
,
D. Y.
Peide
, and
X.
Xu
,
Nat. Commun.
6
,
8572
(
2015
).
22.
F.
Xia
,
H.
Wang
, and
Y.
Jia
,
Nat. Commun.
5
,
4458
(
2014
).
23.
D.
Yin
and
Y.
Yoon
,
J. Appl. Phys.
119
,
214312
(
2016
).
24.
S.
Rakheja
,
Y.
Wu
,
H.
Wang
,
T.
Palacios
,
P.
Avouris
, and
D. A.
Antoniadis
,
IEEE Trans. Nanotechnol.
13
,
1005
(
2014
).
25.
A. V.
Penumatcha
,
R. B.
Salazar
, and
J.
Appenzeller
,
Nat. Commun.
6
,
8948
(
2015
).
26.
I. S.
Esqueda
,
H.
Tian
,
X.
Yan
, and
H.
Wang
,
IEEE Trans. Electron Devices
64
,
5163
(
2017
).
27.
L.
Wang
,
S.
Peng
,
W.
Wang
,
G.
Xu
,
Z.
Ji
,
N.
Lu
,
L.
Li
,
Z.
Jin
, and
M.
Liu
,
J. Appl. Phys.
120
,
084509
(
2016
).
28.
L.
Wang
,
Y.
Li
,
X.
Feng
,
K.-W.
Ang
,
X.
Gong
,
A.
Thean
, and
G.
Liang
, in 2017 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2017), pp. 31.4.1–31.4.4.
29.
J.
Cao
,
S.
Peng
,
W.
Liu
,
Q.
Wu
,
L.
Li
,
D.
Geng
,
G.
Yang
,
Z.
Ji
,
N.
Lu
, and
M.
Liu
,
J. Appl. Phys.
123
,
064501
(
2018
).
30.
Y.
Taur
,
J.
Wu
, and
J.
Min
,
IEEE Trans. Electron Devices
63
,
2550
(
2016
).
31.
E. G.
Marin
,
S. J.
Bader
, and
D.
Jena
,
IEEE Trans. Electron Devices
65
,
1239
(
2018
).
32.
R.
Landauer
,
IBM J. Res. Dev.
1
,
223
(
1957
).
33.
S.
Datta
,
Electronic Transport in Mesoscopic Systems
(
Cambridge University Press
,
1997
).
34.
S.
Rakheja
,
M.
Lundstrom
, and
D.
Antoniadis
, in 2014 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2014), pp. 35.1.1–35.1.4.
35.
M.
Lundstrom
,
S.
Datta
, and
X.
Sun
,
IEEE Trans. Electron Devices
62
,
4174
(
2015
).
36.
T.
Agarwal
,
A.
Szabo
,
M.
Bardon
,
B.
Sorée
,
I.
Radu
,
P.
Raghavan
,
M.
Luisier
,
W.
Dehaene
, and
M.
Heyns
, in 2017 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2017), pp. 5–7.
37.
C.
Crowell
and
S.
Sze
,
Solid State Electron.
9
,
1035
(
1966
).
38.
M. S.
Lundstrom
and
D. A.
Antoniadis
,
IEEE Trans. Electron Devices
61
,
225
(
2014
).
39.
N.
Haratipour
,
M. C.
Robbins
, and
S. J.
Koester
,
IEEE Electron Device Lett.
36
,
411
(
2015
).
40.
N.
Haratipour
and
S. J.
Koester
,
IEEE Electron Device Lett.
37
,
103
(
2016
).
41.
S.
Rakheja
,
M. S.
Lundstrom
, and
D. A.
Antoniadis
,
IEEE Trans. Electron Devices
62
,
2786
(
2015
).
42.
H.
Liu
,
Y.
Du
,
Y.
Deng
, and
D. Y.
Peide
,
Chem. Soc. Rev.
44
,
2732
(
2015
).
43.
H.-Y.
Chang
,
W.
Zhu
, and
D.
Akinwande
,
Appl. Phys. Lett.
104
,
113504
(
2014
).
44.
H.-M.
Chang
,
A.
Charnas
,
Y.-M.
Lin
,
D. Y.
Peide
,
C.-I.
Wu
, and
C.-H.
Wu
,
Sci. Rep.
7
,
16857
(
2017
).
45.
Y.
Du
,
A. T.
Neal
,
H.
Zhou
, and
D. Y.
Peide
,
J. Phys. Condens. Matter
28
,
263002
(
2016
).
46.
H.
Liu
,
A. T.
Neal
, and
P. D.
Ye
,
ACS Nano
6
,
8563
(
2012
).
47.
Y.
Xu
,
C.
Cheng
,
S.
Du
,
J.
Yang
,
B.
Yu
,
J.
Luo
,
W.
Yin
,
E.
Li
,
S.
Dong
,
P.
Ye
et al.,
ACS Nano
10
,
4895
(
2016
).
48.
N.
Haratipour
,
S.
Namgung
,
R.
Grassi
,
T.
Low
,
S.-H.
Oh
, and
S. J.
Koester
,
IEEE Electron Device Lett.
38
,
685
(
2017
).
49.
A.
Valletta
,
A.
Daami
,
M.
Benwadih
,
R.
Coppard
,
G.
Fortunato
,
M.
Rapisarda
,
F.
Torricelli
, and
L.
Mariucci
,
Appl. Phys. Lett.
99
,
271
(
2011
).
50.
S. D. U. Guide, Inc., “Version M-2016.12, Synopsys,” Mountain View, CA (2016).
51.
X.
Chen
,
C.
Chen
,
A.
Levi
,
L.
Houben
,
B.
Deng
,
S.
Yuan
,
C.
Ma
,
K.
Watanabe
,
T.
Taniguchi
,
D.
Naveh
et al.,
ACS Nano
12
,
5003
(
2018
).
52.
G.
Gaddemane
,
W. G.
Vandenberghe
,
M. L.
Van de Put
,
S.
Chen
,
S.
Tiwari
,
E.
Chen
, and
M. V.
Fischetti
,
Phys. Rev. B
98
,
115416
(
2018
).
53.
Z.-Y.
Ong
,
G.
Zhang
, and
Y. W.
Zhang
,
J. Appl. Phys.
116
,
214505
(
2014
).
54.
D.
Jariwala
,
V. K.
Sangwan
,
L. J.
Lauhon
,
T. J.
Marks
, and
M. C.
Hersam
,
ACS Nano
8
,
1102
(
2014
).
55.
B.
Radisavljevic
, and
A.
Kis
,
Nat. Mater.
12
,
815
(
2013
).
56.
D.
Ovchinnikov
,
A.
Allain
,
Y.-S.
Huang
,
D.
Dumcenco
, and
A.
Kis
,
ACS Nano
8
,
8174
(
2014
).
57.
Y.
Trushkov
and
V.
Perebeinos
,
Phys. Rev. B
95
,
075436
(
2017
).
58.
L.
Li
,
Y.
Yu
,
G. J.
Ye
,
Q.
Ge
,
X.
Ou
,
H.
Wu
,
D.
Feng
,
X. H.
Chen
, and
Y.
Zhang
,
Nat. Nanotechnol.
9
,
372
(
2014
).
59.
N.
Haratipour
,
Y.
Liu
,
R. J.
Wu
,
S.
Namgung
,
P. P.
Ruden
,
K. A.
Mkhoyan
,
S.-H.
Oh
, and
S. J.
Koester
,
IEEE Trans. Electron Devices
65
,
1
9
(
2018
).
60.
X.
Li
,
R.
Grassi
,
S.
Li
,
T.
Li
,
X.
Xiong
,
T.
Low
, and
Y.
Wu
,
Nano Lett.
18
,
26
(
2017
).
61.
T.
Li
,
Z.
Zhang
,
X.
Li
,
M.
Huang
,
S.
Li
,
S.
Li
, and
Y.
Wu
,
Appl. Phys. Lett.
110
,
163507
(
2017
).
62.
G.
Arutchelvan
,
C. J. L.
de la Rosa
,
P.
Matagne
,
S.
Sutar
,
I.
Radu
,
C.
Huyghebaert
,
S.
De Gendt
, and
M.
Heyns
,
Nanoscale
9
,
10869
(
2017
).
63.
N.
Pons
,
F.
Triozon
,
M.-A.
Jaud
,
R.
Coquand
,
S.
Martinie
,
O.
Rozeau
,
Y.-M.
Niquet
,
V.-H.
Nguyen
,
A.
Idrissi-El Oudrhiri
, and
S.
Barraud
, in 2013 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD) (IEEE, 2013), pp. 184–187.
64.
P.
Chen
,
N.
Li
,
X.
Chen
,
W.-J.
Ong
, and
X.
Zhao
,
2D Mater.
5
,
014002
(
2017
).
65.
H.
Liu
,
A. T.
Neal
,
Z.
Zhu
,
Z.
Luo
,
X.
Xu
,
D.
Tománek
, and
P. D.
Ye
,
ACS Nano
8
,
4033
(
2014
).
66.
K.
Li
, and
S.
Rakheja
,
J. Appl. Phys.
123
,
184501
(
2018
).
67.
A.
Khakifirooz
,
O. M.
Nayfeh
, and
D.
Antoniadis
,
IEEE Trans. Electron Devices
56
,
1674
(
2009
).
68.
N.
Ma
and
D.
Jena
,
Phys. Rev. X
4
,
011043
(
2014
).
69.
U.
Radhakrishna
, Ph.D. thesis, Massachusetts Institute of Technology, 2016.