Semiconductor Device Modeling and Simulation for Electronic Circuit Design

Abstract

This chapter covers different methods of semiconductor device modeling for electronic circuit simulation. It presents a discussion on physics-based analytical modeling approach to predict device operation at specific conditions such as applied bias (e.g., voltages and currents); environment (e.g., temperature, noise); and physical characteristics (e.g., geometry, doping levels). However, formulation of device model involves trade-off between accuracy and computational speed and for most practical operation such as for SPICE-based circuit simulator, empirical modeling approach is often preferred. Thus, this chapter also covers empirical modeling approaches to predict device operation by implementing mathematically fitted equations. In addition, it includes numerical device modeling approaches, which involve numerical device simulation using different types of commercial computer-based tools. Numerical models are used as virtual environment for device optimization under different conditions and the results can be used to validate the simulation models for other operating conditions.

1. Introduction

Researchers are devoting their time and efforts on the development of efficient and high-speed device models as the requirement for faster, smaller circuits and systems are becoming more and more stringent. These models take into account semiconductor devices such as MOSFETs and analyze the device features through rigorous testing and optimization of the device characteristics. The device models are developed employing various mathematical calculations and simulation processes before being adopted into practical circuit and system-level simulation steps. This book chapter takes a deeper dive into different types of semiconductor device modeling and simulation techniques by leveraging their full potential.

To understand the characteristics of a device, analytical modeling is performed under various environmental and operating conditions by varying the device parameters. Physics-based analytical models are developed by analyzing the behaviors of a device based on the fundamental physics by solving rigorous mathematical equations governing the underlying device physics. This enables selection of optimal device features and predictability of the device operation under various conditions to simulate its performance in electronic circuits and systems. Although this type of modeling provides accurate results, the calculation scheme is not fast enough for higher level analysis tools including circuit simulators such as SPICE. For fast calculation of device characteristics, empirical models are commonly used. Empirical modeling is a computer-based modeling scheme involving experimental data. Unlike analytical modeling approach, this specific technique depends on simulation of device behavior by implementing mathematically fitted equations. This results in a trade-off between the accuracy of device models and the computational speed. Empirical modeling of a device plays a leading role in circuit simulation using tools such as SPICE, which is a general-purpose analog electronic circuit simulator. It has the capability to analyze and predict the device behavior in a circuit design by solving a combination of theoretical and empirical models. However, to explore full potential of a device in implementing efficient circuits, researchers are dedicating much time in implementing numerical models using commercial device simulators such as technology computer-aided design (TCAD) tools. Numerical modeling optimizes the characteristics of a device by simulating the device in virtual environment. This type of modeling of a device involves solutions to a set of coupled partial differential equations. The core concept of this type of modeling involves analyzing various mathematical techniques to provide a realistic solution for predictability and operation of the device. Creating virtual environment under different conditions to simulate various device features and characteristics enables researchers to achieve a grasp on the functionality of a certain device and its characteristics in circuit environment.

2. Physics-based and empirical compact modeling for circuit simulation

2.1 Physics-based models

Physics-based models are developed by simplifying coupled nonlinear partial differential equations that describe the physics of the semiconductor devices. These models need to be robust, accurate, and computationally efficient and are often preferred in analog and RF circuit design. Although modern industry standard models often have some empirical fitting parameters, their core is usually physics-based and most parameters have some physical significance. Since these models try to capture the essence of underlying physics, they usually have better scalability and predictive power compared to their empirical counterparts.

2.1.1 Types of physics-based models

Physics-based compact models can be divided into several subcategories. The most widely used models in industry are BSIM3 [1] and BSIM4 [2]. These models are threshold-based, source-referenced, characterized by inversion charge proportional to overdrive voltage, and use interpolation between strong and weak inversion. In addition, some empirical fitting parameters are used in these models to cover wide bias, geometry, and temperature ranges.

Surface potential-based models such as PSP [3] and HiSIM [4] have recently emerged as other important physics-based compact modeling paradigms. PSP was developed as a collaboration between Gildenblat et al. and Philips corporation and is characterized by both inversion charge and drain current being computed from the surface potential φ_s . This is a symmetric and body-referenced model whose core idea rests upon the symmetric linearization of surface potential, which provides an efficient solution for surface potential with 10 pV accuracy.

Another approach is known as inversion charge-based model such as EKV named after its founders (Enz-Krmmencher-Vittoz) [5] and ACM (Advanced Compact model) [6]. EKV was originally developed at CSEM and EPFL at Switzerland and has found use in low-power and RF integrated circuit (IC) design due to some of the shortcomings of BSIM-based models in these applications. It is a body-referenced and symmetric model based on channel charge linearization. There is also a companion model for hand analysis.

In addition to these popular choices, there are other models such as, analytical model for ballistic field-effect transistors [7], MIT Virtual Source model [8] etc. Although these models are not used in industry, they involve a relatively small set of parameters with strong physical significance and can provide useful insight into the underlying physical phenomena of nanoscale transistors.

2.1.2 Brief overview of previous works

A number of works have been performed over the years on physics-based modeling of transistors. Gummel in [9] developed a model based on finite difference method for solving model equations to provide information about internal parameters such as potential and electric field distribution along with the terminal characteristics. This approach was modified in [10] using a new discretization technique for ensuring convergence. Building upon Scharfetter-Gummel algorithm, Slotboom [11] proposed a new model using two new artificial variables for linearization of differential equations facilitating implementation in computer-aided design (CAD) programs.

Early pioneers of device modeling include Pao and Sah who came up with the classic double integral drain current expression and explored different characteristics of transistor action [12, 13]. However, these formulas are computation intensive and CAD implementation requires a simplified model. Brews in [14] proposed a charge sheet model that assumes the inversion layer to be a conducting plane of zero thickness.

As the transistors kept getting smaller, the significance of subthreshold leakage conduction became apparent and the authors in [15] explored basic charge equation to derive a model that covers conduction mechanism from subthreshold to strong inversion. A one-dimensional model reported in [16] incorporates the dependence of subthreshold conduction on drain voltage, substrate bias, and temperature. The effect of terminal voltages on subthreshold conduction was captured in another model proposed by Taylor [17], which also explored the potential adverse two-dimensional effects on drain conductance. Later, Taylor unified the existing short- and long-channel models into an analytical model in [18]. Moreover, with the increasing use of MOSFETs in the domain of analog circuit design, there has been a growing need for a better small signal CAD models. Liu et al. presented a first-order and a second-order large-signal MOSFET models and derived corresponding small-signal models where the parameters are associated with bias condition and process technology [19].

Several researchers have explored the effect of electric field, doping density, and temperature on carrier mobility. Arora et al. developed an analytical expression of carrier mobilities in silicon as a function of doping concentration and temperature based on experimental data and modified Brooks-Herring theory of mobility [20]. Masetti et al. explored the dependence of carrier mobility on dopant concentration in silicon [21]. Later, Reggiani et al. proposed a unified model that combines mobility dependence on multiple factors such as doping, temperature, and electric field [22]. Another major contribution was the development of a universal expression for carrier generation and recombination, which was independently reported by Hall [23] and Shockley-Read [24] in 1952. Carrier lifetime is the most important parameter affecting the rate of generation-recombination and its dependence on temperature and electric field was analyzed by Schenk in [25].

The first-generation CAD models refer to the three built-in models (Level 1, 2, and 3) in the SPICE2 program [26]. Sheu et al. developed CSIM (Compact Short-channel IGFET Model) [27] and later improved and converted it into its now-famous successor called “Berkeley Short-channel IGFET Model” (BSIM) [28]. This was the beginning of the second generation of SPICE models. Technical considerations underlying the evolution of these models will be explored in greater detail in Section 2.5.

2.2 Empirical model

Empirical modeling is an alternative to physics-based modeling, which is used for rapid and accurate circuit simulation. The basic methodology is to take experimental or simulated device data as input and then create a model capable of accurately reproducing the complex nonlinear behavior of the semiconductor devices under different operating conditions using some sort of special functions such as Black-Box model or table look-up models. Usually, these methods are flexible and can be applied to different types of transistors [29].

2.2.1 Types of empirical models

The most commonly used empirical models utilize look-up table and quadratic or higher order polynomials for interpolation between data points. The data can be obtained from TCAD simulations or experimental results. Some simulators internally use table look-up methods for faster simulation. Currently, there are some mixed-device circuit simulators, which can generate table look-up models from TCAD simulations for carrying out circuit simulation. These methods can have some serious drawbacks such as wiggles in I-V characteristics, nonphysical negative resistances, limited applicability in exponentially varying regions of operation resulting from the use of low-order polynomial interpolation, inability to accommodate changes in temperature and geometry, limited capability for modeling statistical variation and noise, and lack of predictive capabilities for future technologies.

Another method involves using special functions such as hyperbolic tangent in an ingenious way to match the shape of experimental data. Parameters tend to be fitting coefficients but they may have some sort of first-order physical significance. A different approach, known as Black-Box modeling, leverages sophisticated numerical packages that take arbitrary sets of data and automatically generate some form of mathematical model that best fits the supplied data. This approach has mostly been confined to complex interconnect structures and system-level modeling. The boundary between physics-based and empirical models has gradually become fuzzier in recent years since most modern industry standard physics-based models now have some empirical content.

2.2.2 Brief overview of previous works

The table look-up method was implemented for digital circuit simulation in the timing simulator MOTIS [30], which was used to obtain timing information on the propagation of signals through circuits. Rofougaran et al. reported a FET model consisting of a main look-up table, along with a coarse 3-D sub-table for including substrate effects and another table for interpolating between different channel lengths [31]. This was used to resolve some of the limitations faced by analytical models in capturing short-channel effects. Another table-based method was used in [32] to address the requirements of analog circuit design, which demands greater accuracy in reproducing small-signal parameters, large-signal nonlinearities, subthreshold conduction, body effects, and bias-dependent capacitances.

Meijer proposed an n-dimensional model with first-order continuity in [33] where each table model could replicate the DC I-V characteristics of two basic physical device models, namely, the Ebers-Moll model for bipolar transistor and the GLASMOST model for MOSFET with higher accuracy and faster simulation time compared to physics-based compact models. Bourenkov et al. implemented ingeniously combined exponential and polynomial interpolation into a blending function to accurately evaluate drain current in the moderate inversion region [34]. The user can choose among several available interpolation schemes based on speed, memory, and accuracy requirement and the scheme worked well for DC, transient, and AC analyses. Root et al. used table data to develop non-quasi-static FET models for RF simulation [35]. In these models, the characteristics of GaAs FET devices were determined by state functions, which define nonlinear relationships for the three-terminal lumped elements. An array of s-parameters, measured over a wide range of terminal biases, is used to determine state functions, which dictate the characteristics of GaAs devices.

Yanilmaz et al. extended Bernstein approximation technique to multidimensional variation diminishing interpolation and modeled DC I-V and Q-V characteristics of a MOSFET [36]. An important attribute of this model is the preservation of continuity and monotonicity, which are important for convergence in Newton-Raphson algorithm commonly used in most modern circuit simulators. Similarly, the authors in [37] have used Lagrange interpolation and Bernstein approximation techniques for modeling multi-gate SOI transistors. Shima et al. employed monotonic piecewise cubic interpolation on stored table data obtained from a 2-D TCAD simulator to evaluate operating points of MOS transistors [38]. Authors in [39] used a tableau-style quadratic spline formulation that ensured the continuity of the model function and its derivative and proposed a new data-compression method for efficient storage of coefficients. The authors in [40] have used multidimensional linear and cubic spline interpolation methods for modeling SOI four-gate transistor (G4FET).

2.3 CAD model for SPICE

2.3.1 Criteria for a good SPICE model

A comprehensive list of attributes of a good CAD model for circuit simulation can be found in [41]. In this section, a brief discussion on the most salient features of an excellent compact model is presented. To design analog and mixed analog/digital circuits effectively, CAD models have to meet some basic requirements based on I-V characteristics, charges, leakage currents etc.; should provide continuous results that must have physical sense; should meet the requirement for intrinsic and extrinsic effects; and should provide accurate prediction for temperature range of interests etc. In addition, the model must have criteria for any combination of channel length and width values from the minimum specified upward bound and should provide a flag when it is used outside the limit of validity. Finally, and perhaps most importantly, a compact model must be computationally efficient, sufficiently accurate, and numerically robust.

2.3.2 Model formulation

The model developers need to satisfy certain criteria beyond model simulation. The general considerations include choice of parameters [42], choice of references [43], choice of modeling expressions [44], out of range behavior [42], charge versus capacitance formulation [45], and non-quasi-static analysis [46]. Furthermore, there are two important points, namely the choice of smoothing function and the judicious use of conditionals, which are the topics of the following discussion.

Smoothing function is an important tool to build a single expression by combining several expressions. Sometimes it is hard to model a device with general expressions for different regions of operation. The resultant single expression from smoothing function can be used for all the regions. However, the smoothing functions reveal particular expression in particular regions while maintaining continuity at boundaries between regions.

For example, the traditional MOSFET drain current, I_DS with an abrupt transition from non-saturation to saturation operation at a drain-source voltage, V′ _DS can be written as shown in the following equation [42],

$$I_{\mathit{DS}}=\frac{W}{L}\mu C_{\mathit{ox}}^{\prime }\left[\left(V_{\mathit{GS}}-V_{\mathit{TH}}\right)-\frac{\alpha }{2}{V}_{\mathit{DS},\mathit{eff}}\right]$$ E1

where, $V_{\mathit{DS},\mathit{eff}}=\left\{\begin{array}{c}{V}_{\mathit{DS}},\\ V_{\mathit{DS}}^{\prime },\end{array}\right.\begin{array}{c}{V}_{\mathit{DS}}\le V_{\mathit{DS}}^{\prime }\\ V_{\mathit{DS}}>V_{\mathit{DS}}^{\prime }\end{array}$

where μ is the mobility, W and L are the channel width and length respectively, C’_ox is the oxide capacitance, V_TH is the threshold voltage, V_GS is the gate-source voltage, and α is a fitting parameter. Now, rather than using two conditional terms for $V_{\mathit{DS},\mathit{eff}}$ , the equation can be simplified and converted into a single expression using smoothing function. The following equations are based on [47, 48],

$$V_{\mathit{DS},\mathit{eff}}=V_{\mathit{DS}}^{\prime }-0.5\ast \left[V_{\mathit{DS}}^{\prime }-V_{\mathit{DS}}-\delta +\sqrt{{\left(V_{\mathit{DS}}^{\prime }-V_{\mathit{DS}}-\delta \right)}^2+4\delta V_{\mathit{DS}}^{\prime }}\right]$$ E2

and

$$V_{\mathit{DS},\mathit{eff}}=\frac{V_{\mathit{DS}}}{{\left[1+{\left[\frac{\left|V_{\mathit{DS}}\right|}{V_{\mathit{DS}}^{\prime }}\right]}^A\right]}^{1/A}}$$ E3

Here $\delta$ and A are fitting parameters. Figure 1 shows how these functions make the transition region smooth and conserve higher order continuity. Similarly, other functions are used for smooth transition between weak- and strong-inversion [49].

Figure 1.

Illustration of two smoothing functions for ensuring higher order continuity of the model as shown in Eq. (1-3). (Here, $V_{\mathit{DS}}^{\prime }$ =1V, $\delta$ =0.1 V, A=2).

Sometimes, it is hard to maintain perfectly smooth model without using any conditionals. For efficiency, conditional models are used to include or exclude block models that are turned on or off by using different model parameters [47]. For example, let us consider the following common limiting function:

$$y=\frac{1}{A}ln\left[1+exp\left(\mathit{Ax}\right)\right]$$ E4

This function asymptotically approaches zero for large negative x and approaches x for large positive x. But it may cause a numeric overflow for large positive x, which can be solved using the following equivalent conditional equation:

$$y=\left\{\begin{array}{l}\frac{1}{A}ln\left[1+exp\left(\mathit{Ax}\right)\right],x\le 0\\ x+\frac{1}{A}ln\left[1+exp\left(-\mathit{Ax}\right)\right],x>0\end{array}\right.$$ E5

2.3.3 Evolution of different generations of CAD models for SPICE

The “Simulation Program with Integrated Circuit Emphasis” (SPICE) was developed in the Electronics Research Laboratory at the University of California, Berkeley in 1973 [50]. However, the real popularity of SPICE comes after the release of SPICE 2 in 1975, which included three first-generation MOSFET models for the first time [51]. The first-generation compact models for MOSFET include Level 1, Level 2, and Level 3 models. The second-generation models are BSIM, BSIM2, and HSPICE Level 28. The third-generation models include Level 7, Level 48, BSIM3, and other advanced models. The following section presents the chronological progression from the Level 1 model to Level 3 and finally to BSIM models.

Level 1 model: This first SPICE model for MOS transistor, published in 1968, is often referred to as the Shichman-Hodges model [52]. It is the simplest compact model that is only accurate for long-channel (more than ∼10 μm) MOSFETs with uniform doping. This threshold voltage-based model assumes that, when the gate-to-source voltage, V_GS , is greater than or equal to the threshold voltage, V_TH , then the carrier concentration at the surface under the gate oxide gets inverted in polarity with respect to the substrate. The threshold voltage, V_TH as expressed in Eq. (6), is a function of the body effect parameter, γ, the potential difference between the source and the substrate, V_BS , and the bulk potential, 2ϕ_p [53]. In Eq. (6), V_T0 represents the threshold voltage when V_BS = 0.

$$V_{\mathit{TH}}=V_{\mathit{\text{TO}}}+\gamma \left(\sqrt{2{\varphi }_p-V_{\mathit{BS}}}-\sqrt{2{\varphi }_p}\right)$$ E6

In this model, the MOSFET operation is divided into three regions. When V_GS < V_TH , the drain-to-source current, I_DS , is zero. This is called the cutoff region. When V_GS > V_TH , the MOSFET is turned on. In this state, if the drain-to-source voltage, V_DS < V_GS − V_TH , the device is said to be in the linear region and if V_DS > V_GS − V_TH , then the device is operating in the saturation region. The drain current, I_DS in linear and saturation regions is expressed in Eqs. (7) and (8), respectively [53],

$$I_{\mathit{DS}}=K_P\frac{W}{L-2X_{\mathit{jl}}}\left(V_{\mathit{GS}}-V_{\mathit{TH}}-\frac{V_{\mathit{DS}}}{2}\right)V_{\mathit{DS}}\left(1+\lambda V_{\mathit{DS}}\right)$$ E7

$$I_{\mathit{DS}}=\frac{K_P}{2}\frac{W}{L-2X_{\mathit{jl}}}{\left(V_{\mathit{GS}}-V_{\mathit{TH}}\right)}^2\left(1+\lambda V_{\mathit{DS}}\right)$$ E8

In Eqs. (7) and (8), the term, X_jl , refers to the lateral diffusion in source and drain regions, W is the channel width, and K_P is the conduction parameter. The effective channel length is then given by (L−X_jl), where, L is the actual channel length. The term (1+λV_DS) introduces an empirical correction of the conductance in the saturation region, where λ represents the channel length modulation parameter. These five parameters, K_P, V_T0, γ, λ, and 2ϕ_p , characterize the model and can be directly specified in SPICE simulation or can be calculated from physical device parameters including substrate doping, substrate permittivity, gate oxide capacitance, and thickness.

Level 2 model: The Level 1 model does not include the effect of channel dimensions on the threshold voltage. However, the experimental data show that when the channel length is small enough to be comparable with the source and the drain depletion regions, then the relationship between the channel dimension and the threshold voltage is no longer negligible. John E. Meyer addressed this effect along with some other second-order effects to provide a more accurate model for smaller sized devices, which is considered as the Level 2 model [54]. The threshold voltage V_TH is modified in this model, as shown in Eq. (9), by introducing a change in the body effect parameter, γ, as shown in Eq. (10) [53],

$$V_{\mathit{TH}}=V_{\mathit{FB}}+2{\varphi }_p+\gamma \prime \sqrt{2{\varphi }_p-V_{\mathit{BS}}}+\frac{{ϵ}_s\mathit{\delta \pi }}{4C_{\mathit{ox}}^{\prime }W}\left(2{\varphi }_p-V_{\mathit{BS}}\right)$$ E9

$$\gamma \prime =\gamma \left[1-\frac{X_j}{2L_{\mathit{eff}}}\left(\sqrt{1+\frac{2W_S}{X_j}}+\sqrt{1+\frac{2W_D}{X_j}}-2\right)\right]$$ E10

where V_FB represents the flat-band voltage (the amount of applied gate voltage when the channel region and the substrate have the same amount of carrier concentration of the same polarity), ε_s is the substrate permittivity, δ is the width-effect parameter, C’_ox is the oxide capacitance, X_j represents the doping depth in the source and the drain regions, W_S and W_D are the depletion widths at the source and the drain regions, respectively. W_S and W_D both are functions of the substrate bias and the bulk potential. Level 1 model also assumes that the fixed charge in the depleted channel region is independent of the channel-to-substrate voltage. This assumption becomes erroneous for large V _DS , as a result of the significant difference of the depletion widths between the drain and the source regions. Taking this effect into account, the drain current, I_DS expression in the linear region is modified, as presented in Eq. (11), which gives current values close to those of the Level 1 model for smaller V_DS .

$$\begin{array}{ll}{I}_{\mathit{DS}}=& \frac{K_P}{\left(1-\lambda V_{\mathit{DS}}\right)}\frac{W}{\left(L-2X_{\mathit{jl}}\right)}\left\{\left(V_{\mathit{GS}}-V_{\mathit{FB}}-2{\varphi }_p-\frac{V_{\mathit{DS}}}{2}\right)V_{\mathit{DS}}\right.\\ & \left.-\frac{2}{3}\gamma \left[{\left(V_{\mathit{DS}}-V_{\mathit{BS}}+2{\varphi }_p\right)}^{1.5}-{\left(-V_{\mathit{BS}}+2{\varphi }_p\right)}^{1.5}\right]\right\}\end{array}$$ E11

The drain voltage in the saturation region is expressed in Eq. (12). The drain current in the saturation region is obtained from Eq. (11) at V_DS = V_D,sat . A significant modification from Level 1 is the variation of the drain current with γ even if V_BS = 0. The drain-to-source voltage at saturation can be expressed as,

$$V_{D,\mathit{sat}}=V_{\mathit{GS}}-V_{\mathit{FB}}-2{\varphi }_p+{\gamma }^2\left[1-\sqrt{1+\frac{2}{{\gamma }^2}\left(V_{\mathit{GS}}-V_{\mathit{FB}}\right)}\right]$$ E12

The Level 2 model offers other modifications including nonzero drain current in the weak-inversion region when V_GS < V_TH , and voltage-dependence of K_P , reflecting the change of the carrier mobility with the gate and the drain voltages.

Level 3 model: The basic equations of the Level 3 model were proposed by Dang in 1979 [55]. This model can successfully predict the characteristic of a device with a channel length down to 2 μm. In the linear region, the drain current expression of Level 2 is simplified in this model by using the Taylor series expansion, as expressed in Eq. (13),

$$I_{\mathit{DS}}=\beta \left(V_{\mathit{GS}}-V_{\mathit{TH}}-\frac{1+F_B}{2}{V}_{\mathit{DS}}\right)V_{\mathit{DS}};F_B=\frac{\gamma F_S}{2\sqrt{2{\varphi }_p-V_{\mathit{BS}}}}+F_n$$ E13

The short-channel effect influences the empirical parameters, F_S and β, while the narrow channel influences another fitting parameter F_n . Level 3 uses a hypothesis, similar to Level 2, in formulating the threshold voltage. The threshold voltage proposed in Level 3 is presented in Eq. (14). However, the influence of the difference in the depletion width between the source and the drain regions, at higher V_DS , is now empirically expressed with the parameter, σ.

$$V_{\mathit{TH}}=V_{\mathit{FB}}+2{\varphi }_p-\sigma V_{\mathit{DS}}+\gamma F_S\sqrt{2{\varphi }_p-V_{\mathit{BS}}}+F_n\left(2{\varphi }_p-V_{\mathit{BS}}\right)$$ E14

This model handles the gate voltage dependence of the surface mobility, μ_s , and the drain-to-source electric field dependence of the effective mobility, μ_eff , in a simpler way, as expressed in Eq. (15), by using the mobility modulation parameter, θ.

$${\mu }_s=\frac{\mu }{1+\theta \left(V_{\mathit{GS}}-V_{\mathit{TH}}\right)};{\mu }_{\mathit{eff}}=\frac{{\mu }_s}{1+{\mu }_s\left(\raisebox{1ex}{$V_{\mathit{DS}}$}\!\left/ \!\raisebox{-1ex}{$v_{\mathit{max}}{L}_{\mathit{eff}}$}\right.\right)}$$ E15

In Eq. (15), the term v_max is the velocity limit reached by the carrier when V_DS = V_D,sat . The Level 3 version of the saturation voltage V_D,sat is presented in Eq. (16).

$$V_{D,\mathit{sat}}=\frac{V_{\mathit{GS}}-V_{\mathit{TH}}}{1+F_B}+\frac{v_{\mathit{max}}{L}_{\mathit{eff}}}{{\mu }_s}-\sqrt{{\left(\frac{V_{\mathit{GS}}-V_{\mathit{TH}}}{1+F_B}\right)}^2+{\left(\frac{v_{\mathit{max}}{L}_{\mathit{eff}}}{{\mu }_s}\right)}^2}$$ E16

BSIM1 model: To bring higher accuracy in modeling shorter channel devices (down to 1 μm), the Berkeley Short Channel IGFET Model 1 (BSIM1) was introduced in 1987 [28]. Similar to Level 2 and Level 3 models, this model incorporates both the strong- and weak-inversion components of the drain current. However, unlike the previous models, BSIM1 used numerical approximation for modeling the dependence of the drain current on the substrate bias to speed up the simulation process. An automated parameter extraction program was designed to extract the model parameters. BSIM1 incorporates an improved formulation of short-channel effects to deal with short-channel devices. To enhance the scalability of the model, several fitting parameters are introduced for each of the model parameters. However, this model offers no significant improvement in scalability and due to a large number of fitting parameters, it has failed to gain popularity among circuit designers.

BSIM2 model: BSIM2 was developed in 1990 for submicron devices by improving some aspects of BSIM1 including model continuity, output conductance, and subthreshold current [56]. It employs cubic spline to achieve smoother transition between weak- and strong-inversion and between linear and saturation regions. BSIM2 uses more parameters for improved accuracy but the model has not solved the issue regarding the difficult parameter extraction process. The user still cannot get a set of parameters that are valid for a range of device sizes and has to deal with many sets of model parameters, each covering a limited range of device geometries.

BSIM3 model: BSIM3 was developed in 1994 from a coherent quasi-two-dimensional analysis of the MOSFET by addressing these issues of previous models [49]. To ensure good scalability, this model explicitly considers the effect of device size and process variation. The second version, BSIM3v2, was released with better accuracy and scalability than the previous version. However, this second version still suffers from the discontinuity, such as negative conductance and glitches in transconductance over drain current versus gate voltage plot at the boundary of weak- and strong-inversion regions. The third version, BSIM3v3.0, was introduced to eliminate all the shortcomings of the previous versions [57]. The significant distinction of BSIM3v3.0 is the introduction of the single-equation approach and prediction-ability that enabled statistical analysis [58]. The improvement of this third version of BSIM3 continued to three consecutive subversions, from BSIM3v3.0 to BSIM3v3.2 to incorporate significant developments including the introduction of a new charge-capacitance model considering the quantization effect, improved threshold model, substrate current model, and non-quasi-static model [59].

BSIM4 model: BSIM 4 was released in 2000 to support sub 130-nmm CMOS technologies and the growth of high-speed analog, mixed-signal, and RF integrated circuits [60]. A significant addition is a holistic noise model for the channel thermal noise and induced gate noise [61]. An intrinsic input resistance model is used to obtain accuracy at high-frequency operation. It also introduces a charge layer thickness model incorporating novel quantum effects [62]. Moreover, it includes the first model that takes into account the drain leakage current resulting from direct-tunneling [63]. Other improvements include pocket implant effect, layout-dependent factors including mechanical stress and well-proximity and enabled modeling of high-k metal-gate stacks and non-silicon materials.

BSIM-bulk model: There have been recent efforts to combine the best features of BSIM4 and EKV models. This has resulted in the newest addition to the Bulk MOSFET model from BSIM group and is named BSIM-bulk (formerly known as BSIM6). In contrast to its threshold voltage-based predecessors, BSIM-bulk is a charge-based and body-referenced model, which passes the symmetry test for DC and AC, correctly predicts harmonic slope, and exhibits accurate results for RF and analog simulations [64]. In addition, bulk charge effect has been modeled analytically to improve the model accuracy for transconductance and output conductance [65]. There is also a new NQS (non-quasi-static) model effective up to the millimeter wave regime.

2.4 Modeling of novel transistors and emerging devices

To continue the advancement of semiconductor industry beyond the imminent demise of Moore’s Law, researchers all around the world have been exploring new technologies and consequently, a multitude of novel devices have emerged in recent years. This book chapter has been focused mainly on the development of compact models for bulk MOSFET. These models have been modified and extended for several novel transistors. For example, for SOI (silicon-on-insulator) devices, PSP-based model [66] and BSIM-based model [67, 68] have been reported. From the BSIM group, BSIM-IMG model for independent multi-gate MOSFET operation [69] and BSIM-CMG for common multi-gate transistor, that is, FinFET have been developed [70].

In addition to different nanoscale transistors, the lessons learned during development of different generations of compact models have also been successfully used to model many emerging devices such as metal-oxide-based resistive random-access memory (RRAM) devices [71, 72], insulator-metal-transition devices [73], biomolecular memristors [74] and memcapacitors [75] etc.

3. Device simulation

For any newly developed semiconductor device, different parameters need to be optimized before the device model can be adopted for practical circuit applications. However, with rapid development of new device and process technologies, optimization of semiconductor manufacturing processes guided by experimental approach becomes very time-consuming and expensive. As an alternative, device simulation can allow optimization of the device parameters in a virtual environment in a fast and cost-effective way to verify the device models for a newly developed semiconductor technology. In general, the tools used for numerical device simulation include three major components: (i) simulation of the fabrication process, (ii) simulation of the device characteristics, and (iii) simulation of the device for circuit applications. Figure 2 shows the basic hierarchy of process, device and circuit simulation. Process simulation is closely coupled to the original device simulation because the behavior of the fabricated device is significantly related to the overall processing steps. Process simulation is implemented by mimicking the original device fabrication steps that include several lithography steps as well as ion implantation, diffusion, annealing, and oxidation steps. This simulation is based on the measurement of doping profile using secondary ion mass spectroscopy (SIMS), topography provided by transmission electron microscopy (TEM), the process recipe, and the lithography masks. Simulation of the device characteristics is based on the overall geometry of the device coupled with the simulated profile from the process simulation. This part of the simulation mainly focuses on the device output and transfer current-voltage (I-V) characteristics, capacitance-voltage (C-V) characteristics, or frequency response. The compact model intended to be used for the circuit simulation can be developed based on the output of the process and device simulation data optimized using different device parameters in TCAD tools.

Figure 2.

Hierarchy of process, device, and circuit simulation.

3.1 Formulation of device simulation models

Semiconductor device simulation usually follows two different approaches such as semiclassical approach and quantum mechanical formulations. In the semiclassical approach, Boltzmann transport equation is used to model the carrier transport in the semiconductor devices by developing the drift-diffusion model along with the energy transport model. The solution of Boltzmann transport equation is usually obtained from doping profile of the device structure. The electrostatic potential, which depends on the dopant profile, is obtained from the Poisson’s equation. In order to incorporate quantum mechanical transport phenomena in the modern device models, Boltzmann transport equation is coupled with Schrödinger equation and Wigner function [76, 77]. Device simulation using the classical approach is discussed in next section, which is then followed by a discussion of the quantum mechanical approach.

3.1.1 Device simulation using semiclassical approach

In the semiclassical approach, the Boltzmann transport equations to describe the transport of electron and holes in a device can be written as,

$$\frac{\partial f}{\partial t}+\mathit{v}.{\nabla }_{r}f-\frac{q\mathit{E}}{\mathrm{\hslash }}\bullet {\nabla }_{\mathit{k}}f=Q\left(f\right)$$ E17

where f(r,k,t) stands for carrier distribution, which is a function of space, r, momentum, k, and time, t; v is the velocity; and E is the electric field. Q(f) denotes the collision operator as a function of the carrier distribution that takes into account different scattering phenomena of the particle in the presence of impurities, phonon, interfaces, and scattering from other sources. However, solving the equation in this form can be computationally intense and thus usual practice is to solve it by applying approximate methods. One such method is called method of moments, which yields a set of differential equations as a function of time and space after multiplying a weight function with each term [78, 79]. The definition of the moments of the distribution function is written as,

$$⟨\Phi ⟩=\frac{1}{4{\pi }^3}\int \mathit{\Phi f}\left(r,k,t\right)d^{3}k$$ E18

The drift-diffusion model can be obtained, by applying the method of moment in the Boltzmann transport equation. For the derivation of the drift-diffusion model, the first two moments such as Φ₀ = 1 and Φ₁ = ℏk of the distribution function are multiplied and integrated over k space. The integral of the collision operator Q can be approximated by applying macroscopic relaxation time (RTA) along with further simplification using the relationship of parabolic dispersion. Using this method, the following set of differential equations for the drift-diffusion model can be obtained,

$$\nabla .J_n=\mathit{qR}+q\frac{\partial n}{\partial t}$$ E19

$$\nabla .J_p=-\mathit{qR}-q\frac{\partial p}{\partial t}$$ E20

$$J_n=\mathit{qn}{\mu }_{n}E+q{D}_n\nabla n$$ E21

$$J_p=\mathit{qp}{\mu }_{p}E-q{D}_p\nabla p$$ E22

where J is the current density, R is the net recombination rate, n and p are the concentration of electron and hole, respectively, μ is the mobility, E is the electric field, and D is the diffusion coefficient. These set of equations for the drift-diffusion model combined with the Poisson’s equation build the basic platform for the semiconductor device simulation. The charge transport in the semiconductor device with respect to the electrostatic potential is obtained by using the Poisson’s equation as given by,

$$\nabla \bullet \left({\epsilon }_s\nabla \phi \right)=q\left(n-p-N_D-N_A\right)$$ E23

where φ is the electrostatic potential, ε_s is the dielectric permittivity, and N_D and N_A are the doping concentration of donor and acceptor atoms, respectively.

3.1.2 Device simulation using quantum mechanical approach

The basic simulation model using the semiclassical approach does not consider the quantum mechanical properties present in the semiconductor device. However, with the continuous reduction of the device size and due to the dual wave-particle nature of electron, quantum mechanical effects in the semiconductor devices have become significant and thus need to be incorporated in the model. There have been several procedures developed to formulate the quantum mechanical transport in modern electronic devices. Figure 3 shows a flow chart representing the relationships among the relevant formulation of quantum models including Schrödinger equation, transfer-matrix, density matrix, Green’s functions, Wigner function, and path integral approaches [80].

Figure 3.

Flow chart illustrating the formulation of the quantum mechanical transport.

One of the most common ways to incorporate the quantum mechanical effect with the already developed semiclassical model involves coupling of Boltzmann transport equation with Schrödinger equation inside a self-consistent Schrödinger-Poisson loop. In this approach, the carrier concentration and the electrostatic potential are obtained by using Schrödinger equation and Poisson equation, respectively, in several iterations until a self-consistent solution is obtained. Boltzmann transport equations are then solved using the derived carrier concentration [81].

Another approach involves using the Wigner function along with Boltzmann transport equation to derive Boltzmann-Wigner equation [82],

$$\frac{\partial f}{\partial t}+v.{\nabla }_{r}f-\frac{q}{\mathrm{\hslash }}\sum _{\alpha =0}^{\infty }\frac{{\left(-1\right)}^{2\alpha }}{4^{\alpha }\left(2n+1\right)!}{\nabla }_r^{2n+1}V\left(r\right)\bullet {\nabla }_k^{2n+1}f={\left(\frac{\partial f}{\partial t}\right)}_C$$ E24

where V denotes an external potential. This equation can be reduced to the classical Boltzmann equation by considering the term α = 0. On the other hand, considering α = 1 results in the density-gradient model as given by [83, 84]

$$\frac{\partial f}{\partial t}+\frac{\mathrm{\hslash }\bullet k}{m^{\ast }}{\nabla }_{r}f-\frac{1}{\mathrm{\hslash }}{\nabla }_r\left(V\left(r\right)-\frac{{\mathrm{\hslash }}^2}{12m^{\ast }}{\nabla }_r^{2}ln\left(n\right)\right){\nabla }_{k}f={\left(\frac{\partial f}{\partial t}\right)}_C$$ E25

where m^* is the effective mass of the carrier. By using this equation, the quantum drift-diffusion model can be obtained in a similar way by applying the method of moments described in the previous section [85].

$$n=N_{c}exp\left(\frac{E_f-E_c-\Lambda }{k_{B}T}\right)$$ E26

$$J_n=-{\mu }_n{k}_{B}T\nabla n-{\mu }_{n}n\nabla \left(E_c-k_{B}Tln{N}_c+\Lambda \right)$$ E27

$$\Lambda =-\frac{\gamma {\mathrm{\hslash }}^2}{12m^{\ast }}\left({\nabla }^{2}lnn+\frac{1}{2}{\left(\nabla lnn\right)}^2\right)=-\frac{\gamma {\mathrm{\hslash }}^2}{6m^{\ast }}\frac{{\nabla }^2\sqrt{n}}{\sqrt{n}}$$ E28

where γ and Λ are used as correction factors, E_f and E_c are the Fermi level and the conduction energy band respectively, k_B is the Boltzmann constant, and T is the temperature. The implementation of the density-gradient model represents local quantum effect, which is in many cases more convenient to implement in a numerical device simulator than Schrödinger-Poisson equation that depends on nonlocal quantities [86, 87, 88, 89, 90]. This is due to the fact that for most numerical simulations, the overall device structure is divided into infinitesimal meshes and the electrical properties are then numerically calculated considering propagation along the meshes. However, there have been some studies that show that although this approach is suitable to model the carrier concentration of the inversion layer of the MOSFET, it fails to model the tunneling currents, which represent other important quantum mechanical effects present in modern devices [85]. Application of nonequilibrium Green’s function is another approach to formulate quantum system that has nonvanishing boundary conditions for Schrödinger equations. However, the solution of the Green’s function quantum transport equation is very complex and several assumptions and approximations are applied to simplify the derivation [91, 92] and the carrier concentration can be calculated as,

$$D=\frac{m^{\ast }{k}_{B}T}{2{\pi }^2{\mathrm{\hslash }}^2}\int A\left(E\right)ln\left(1+exp\left(\frac{E_f-E_c}{k_{B}T}\right)\right)\mathit{dE}$$ E29

3.2 TCAD device simulation

The numerical models described in the previous section are very computationally intensive and iterative in nature. For this reason, device simulation is usually performed using technology computer-aided design (TCAD) tools, which can provide excellent predicting capability of the device properties allowing virtual prototyping and optimization [93, 94, 95]. Most of the modern TCAD packages consist of several tools to implement device processing, device design, parameter extraction, device and circuit simulation, and data visualization as shown in Figure 4. Several of these tools can be combined together based on designer preference in order to study the impact of any single step on the overall system performance. Among these, the process simulation tool creates a virtual environment to emulate the original fabrication and processing recipes that allow the process engineer to study each processing step on the device characteristics and thus facilitate fine tuning of their recipe to optimize the device performance. Electrical, thermal, and optical properties of semiconductor devices are analyzed using another dedicated tool for device simulation. Most TCAD device simulation tools implement finite element methods and some of them have the capability to extract SPICE model parameters for implementation in circuit applications. Physical equations and material properties are incorporated for better prediction while considering the convergence speed. The simulated data are stored in a standard format, which can be represented by the visualization tools for further analysis.

Figure 4.

Principle tools inside TCAD device simulation suites with input and output files.

4. Conclusion

In the era of rapid development of device technology and electronic circuit design, semiconductor device modeling plays a vital role. To build sophisticated electronic circuits, one needs to take a deeper dive into understanding the features and the characteristics of the device through various modeling techniques. In this chapter, different approaches of device modeling for electronic circuit design have been discussed. To understand the behavior of any type of device, the first and foremost approach is to understand the underlying physics for the device operation. Physics-based analytical models provide accurate description of a device formulation using mathematical equations with different operating parameters. However, with the development of more complicated device structure, this type of technique is becoming more complex and computationally inefficient. To achieve a computationally efficient process, empirically fitted parameters are introduced in the original physics-based equations. In order to predict the device behavior in a circuit application, SPICE simulation with efficient, accurate, and robust device model becomes essential. With continuous research breakthroughs and introduction of new physical phenomena, compact device models have evolved from the first-generation Level 1 model to more sophisticated BSIM models. The newly developed device models include a number of device parameters not only to describe the fundamental device characteristics but to include other secondary phenomena that have become significant with continuous device scaling and modification. In addition to the development of efficient device models for wide range of circuit applications, it is also important for device simulation to optimize the fabrication parameters, prior to actual processing. This is to keep pace with the rapid development of newly emerging devices. This book chapter attempts to provide a brief overview of different aspects and methods for device modeling and simulation for electronic circuit design.