Abstract

We propose a simple model, derived from Pao-Sah theory, valid in all modes from weak to strong inversion, to calculate the drain current in Metal Oxide Semiconductor Field Effect Transistor (MOSFET). The Pao-Sah double integral is decomposed into single integrals with limits of integration calculated from Taylor polynomials of inverse functions. The solution is presented analytically wherever possible, and the integration is made from simple numerical methods (Simpson, Romberg) or adaptative algorithms and can be implemented in simple C-program or in usual mathematical software. The transconductance and the diffusion current are also calculated with the same model.

1. Introduction

The Metal Oxide Semiconductor Field Effect Transistor (MOSFET), first proposed in 1926 by Lilienfeld [1], is considered one of the most widely used electronic devices, particularly in digital integrated circuits due to the attractivity of Complementary Metal Oxide Semiconductor (CMOS) logic, a device using energy only during transition states. Since the early days and its simplified model of a conductive channel when the gate voltage 𝑉𝑔 is above a threshold voltage 𝑉𝑇 [2], the MOSFET has certainly been one of the electronic devices which received the most extensive attention from the microelectronic design community. The requirement in the knowledge of the MOSFET working with a gate voltage just below threshold, when the device is still conducting, has induced a complete modeling of drain current in all conductive modes from strong to weak inversion. The first attempt was in 1966, from Pao and Sah [3], who gave an expression of the drain current in a double integration format derived from the equation of the inversion charges  versus surface potential. Pierret and Shields gave in 1983 a single integral expression based on the derivative of the electric field [4]. In 1995, Persi and Gildenblat proposed a computational calculation of the double integral by numerical treatment of carrier concentration and surface potential [5]. Recently, a complete history of MOSFET modeling [6] has given a reference for users of MOS transistors.

The MOSFET modeling is now so much well covered and addressed in BSIM, EKV, and PSP compact models [8] that yet-another-paper on the topic of MOSFET models among the thousand of previous works could appear unnecessary. But, if the state-of-the-art MOSFET models [8] are always the reference in Computer Aided Design (CAD), it could also be interesting to present a semianalytic resolution of Pao-Sah integral which can be implemented in usual commercial software, giving most results quite instantaneously.

Although most of the resolutions of Poisson-Boltzmann equation use finite element software [9–11], we prefer to focus our attention on the physics of the surface potential in order to estimate the effect of gate/drain bias. In the same time, an analytic resolution has the advantage to highlight the influence of the different physical parameters on the derived characteristics. We previously used a similar method in the analytic description of scanning capacitance microscopy based on silicon surface potential [12].

Under thermal equilibrium conditions, mobile charge densities are exponential functions of the potential distribution and this leads to a nonlinear differential equation for the potential πœ‘(π‘₯) in the form of the well known Poisson-Boltzmann equation. Under the Gradual Channel Approximation (GCA) [13], this equation is solved analytically and gives the exact 1-D electric field 𝐹(π‘₯). The gate voltage 𝑉𝑔  versus surface potential πœ‘π‘† explicit equation is inversed using Taylor expansion with an iterative step compatible with the necessary accuracy of integral simulation. The same algorithm is used to estimate the surface potential  versus the channel potential 𝑉(𝑦) induced by drain bias.

We qualify our approach being β€œsemianalytic” in the sense that we solve all equations analytically up to the point where analytic calculations can be done, then decompose the Pao Sah double integral into single integrals and finally we converge to an approximated solution by using simple iterative algorithms which can easily be implemented in usual commercial mathematical worksheet tools or even encoded in simple C programs. This approach is equally well suited for calculation (however not restricted to) of the drift current, the transconductance and the diffusion current, as is further discussed in next sections, since those quantities share the same type of integration methodologies.

2. Surface Potential in MOSFET

In this paper, like in most other papers dealing with analytic description of surface potential, we base our derivations on the gradual channel approximation [13]. In this method, the electric field magnitude in the 𝑦 direction parallel to the conducting channel is assumed to be much smaller than in the vertical π‘₯ direction toward silicon bulk allowing the formulation div𝐹=𝑑𝐹/𝑑π‘₯ and 𝐹[π‘₯,𝑦]=βˆ’π‘‘[πœ‘(π‘₯)]/𝑑π‘₯. The Poisson-Boltzmann equation can then be solved analytically according to the 1-D model of Nicollian and Brews [14] from the complete charge density: 𝜌(π‘₯,𝑦)=π‘ž[𝑝(π‘₯,𝑦)βˆ’π‘›(π‘₯,𝑦)+π‘π·βˆ’π‘π΄].

In the following, electrons distribution in the channel will be considered out of thermal equilibrium with the introduction of the quasi-Fermi energy [15]. The drift and diffusion electron currents are shown in Figure 1. Source and bulk will always be considered to be grounded (except in Section 9) and gate voltage 𝑉𝑔 is defined relative to flat band.

The presentation relates to a 𝑝-type semiconductor (𝑛-MOSFET), but this is not restrictive and could easily be extended to 𝑛-type semiconductor. Upon usual notation conventions (see Nomenclature section), the bulk Fermi level is πœ‘π‘. πœ‘(π‘₯,𝑦) is the potential in every points of the semiconductor material and 𝑉(𝑦) is the channel voltage which defines the quasi-Fermi level offset from equilibrium for electrons. Upon application of a drain bias 𝑉𝐷, 𝑉(𝑦) ranges from 0 at the source side (𝑦=0) to 𝑉𝐷 at the drain side (𝑦=𝐿).π‘ˆπ‘‡=𝐾𝑇/π‘ž is the thermal voltage.

In the following, we use the reduced potentials 𝑒𝑏=πœ‘π‘/π‘ˆπ‘‡,𝑒(π‘₯,𝑦)=πœ‘(π‘₯,𝑦)/π‘ˆπ‘‡ and πœ‰(𝑦)=𝑉(𝑦)/π‘ˆπ‘‡ instead of πœ‘π‘,πœ‘(π‘₯,𝑦) and 𝑉(𝑦), respectively. 𝑒𝑆(𝑦)=𝑒(π‘₯,𝑦)|π‘₯=0,𝐹𝑆(𝑦)=𝐹(π‘₯,𝑦)|π‘₯=0,𝑛𝑆(𝑦)=𝑛(π‘₯,𝑦)|π‘₯=0 are shorthand notations for surface quantities.

In the first MOSFET modeling, charge densities have been written by introducing the quasi-Fermi level only in 𝑛(π‘₯,𝑦) with the expressions:𝑛(π‘₯,𝑦)=𝑛𝑖𝑒𝑒(π‘₯,𝑦)βˆ’πœ‰(𝑦),(1)𝑝(π‘₯)=π‘›π‘–π‘’βˆ’π‘’(π‘₯,𝑦)𝑁,(2)𝐴=π‘›π‘–π‘’βˆ’π‘’π‘π‘,(3)𝐷=𝑛𝑖𝑒𝑒𝑏.(4)

After solving the 1-D Poisson equation 𝑑𝐹(π‘₯,𝑦)/𝑑π‘₯=𝜌(π‘₯,𝑦), these expressions give the 1-D electric field by=ξƒŽπΉ(π‘₯,𝑦)2πΎπ‘‡π‘›π‘–πœ€π‘†Γ—ξ”π‘’βˆ’πœ‰(𝑦)𝑒𝑒(π‘₯,𝑦)βˆ’π‘’π‘’π‘ξ€Ύ+π‘’βˆ’π‘’(π‘₯,𝑦)βˆ’[π‘’βˆ’π‘’π‘βˆ’π‘’π‘’π‘]ξ€Ίπ‘’π‘ξ€»βˆ’π‘’(π‘₯,𝑦)βˆ’π‘’βˆ’π‘’π‘,(5) and its derivative𝑑𝐹(π‘₯,𝑦)π‘‘πœ‰(𝑦)=βˆ’2πΎπ‘‡π‘›π‘–πœ€π‘†π‘’βˆ’πœ‰(𝑦)𝑒𝑒(π‘₯,𝑦)βˆ’π‘’π‘’π‘ξ€»2𝐹(π‘₯,𝑦).(6)

As 𝑒𝑒𝑏β‰ͺ𝑒𝑒(π‘₯,𝑦) in the channel region, this expression is equivalent to𝑑𝐹(π‘₯,𝑦)π‘‘πœ‰(𝑦)=βˆ’2πΎπ‘‡πœ€π‘†π‘›(π‘₯,𝑦)2𝐹(π‘₯,𝑦),(7)

This later expression gives a simple path to the Pao-Sah integral [4].

But, as this has recently been emphasized and extensively discussed in [16], since quasi-neutrality imposes that 𝑛(π‘₯β†’βˆž)=𝑁𝐷, (4) must be corrected by𝑁𝐷=π‘›π‘–π‘’π‘’π‘βˆ’πœ‰(𝑦).(8)

With this exact formulation of the donor density, the electric field gradient from the 1-D Poisson equation (5) now reads𝑑𝐹(π‘₯)=𝑑π‘₯π‘žπ‘›π‘–πœ€π‘†ξ€·π‘’βˆ’π‘’(π‘₯,𝑦)βˆ’π‘’π‘’(π‘₯,𝑦)βˆ’πœ‰(𝑦)+π‘’π‘’π‘βˆ’πœ‰(𝑦)βˆ’π‘’βˆ’π‘’π‘ξ€Έ,(9) which gives the electric fieldπΉξƒŽ(π‘₯,𝑦)=2πΎπ‘‡π‘›π‘–πœ€π‘†βˆšπ‘’βˆ’πœ‰(𝑦)𝐻[𝑒][𝑒](π‘₯,𝑦)+𝐺(π‘₯,𝑦),(10) where we have defined𝐻(𝑒)=𝑒𝑒+𝑒𝑒𝑏𝑒𝑏,βˆ’π‘’βˆ’1𝐺(𝑒)=π‘’βˆ’π‘’βˆ’π‘’βˆ’π‘’π‘ξ€Ίπ‘’π‘ξ€».βˆ’π‘’+1(11)

In this form, the derivative𝑑𝐹(π‘₯,𝑦)π‘‘πœ‰(𝑦)=βˆ’2πΎπ‘‡πœ€π‘†π‘’βˆ’πœ‰(𝑦)𝐻[]𝑒(π‘₯,𝑦)2𝐹(π‘₯,𝑦),(12) is no longer with a numerator proportional to 𝑛(π‘₯,𝑦) as in (7) and the Pao-Sah double integral must now be studied by other methods as outlined in Section 5. In the following, the surface electric field 𝐹𝑆(𝑦)=𝐹(π‘₯,𝑦)|π‘₯=0 is given by substituting 𝑒 and 𝑒(π‘₯,𝑦) by 𝑒𝑆(𝑦) in (10)-(11).

3. The Surface Potential Dependence to Gate Bias

From electrostatic considerations [14], the gate voltage is expressed by𝑉𝑔=𝛾0ξ”π‘’βˆ’πœ‰(𝑦)𝐻𝑒𝑆𝑒(𝑦)+𝐺𝑆(𝑦)+π‘ˆπ‘‡ξ€Ίπ‘’π‘†ξ€»(𝑦)βˆ’π‘’(𝑏).(13)

The exponential expression in (13) allows to separate 𝑒(𝑠) and πœ‰(𝑦)π»ξ€Ίπ‘’πœ‰(𝑦)=ln𝑆(𝑦)𝐸𝑒𝑆(𝑒𝑦)βˆ’πΊπ‘†(𝑦),(14) with the derivativeπ‘‘πœ‰(𝑦)=𝑒𝑑𝑒(𝑠)𝑑𝐻𝑆(𝑦)/𝑑𝑒𝑆(𝑦)𝐻𝑒𝑆(ξ€»βˆ’ξ€Ίπ‘’π‘¦)𝑑𝐸𝑆(𝑦)/𝑑𝑒𝑆𝑒(𝑦)βˆ’π‘‘πΊπ‘†ξ€»(𝑦)/𝑑𝑒𝑆(𝑦)𝐸𝑒𝑆(𝑒𝑦)βˆ’πΊπ‘†(ξ€»,𝑦)(15) in which we introduce the dimensionless quantity𝐸(𝑒)=π‘‰π‘”βˆ’π‘ˆπ‘‡ξ€Ίπ‘’βˆ’π‘’π‘ξ€»π›Ύ0ξƒ°2.(16)

Several approximations to (13) have been reported. For instance, the equation𝑉𝑔2ξƒŽ=Ξ¨(𝑠)+𝛾𝑛𝑖𝑁𝐴2ξ‚ΈexpΞ¨(𝑠)βˆ’π‘‰(𝑦)π‘ˆπ‘‡ξ‚Ή+Ξ¨(𝑠)π‘ˆπ‘‡,(17) is often used in practice in a large variety of surface potential-based models with the band-bending: Ξ¨(𝑠)=πœ‘(𝑠)βˆ’πœ‘π‘. It provides a sufficiently accurate solution for the surface potential in inversion mode but generally lacks accuracy near the flat-band conditions: 𝑒(𝑠)=𝑒𝑏.

4. The Surface Potential Dependence to Drain Bias

The explicit relation πœ‰(𝑦)=𝑓[𝑒𝑆(𝑦)] is given above from (14). According to the importance of this expression in the Pao-Sah double integral, we must pay a particular attention to πœ‰(𝑦)=𝑓[𝑒𝑆(𝑦)] and to the inverse function 𝑒𝑆(𝑦)=π‘“βˆ’1[πœ‰(𝑦)] which will be used in the second integration of the Pao-Sah double integral.

Unfortunately, (14) cannot be inversed by mathematical function and needs a special treatment which is presented next, after having fixed the limits of 𝑒𝑆(𝑦) compatible with the denominator 𝐸[𝑒𝑆(𝑦)]βˆ’πΊ[𝑒𝑆(𝑦)].

4.1. Boundary Limits of Surface Potential at Constant Gate Voltage

Equation (14) has only physical meaning when 𝐸[𝑒𝑆(𝑦)]βˆ’πΊ[𝑒𝑆(𝑦)]>0. At a constant 𝑉𝑔 value, the reduced surface potential 𝑒𝑆(𝑦) ranges between 𝑒Low (lower value), solution of (13) with πœ‰(𝑦)=0𝑉𝑔=𝛾0𝐻𝑒Low𝑒+𝐺Lowξ€Έ+π‘ˆπ‘‡ξ€·π‘’Lowβˆ’π‘’π‘ξ€Έ,(18) and between 𝑒Up (upper value), solution of (13) with πœ‰(𝑦)β†’βˆž leading to π‘’βˆ’πœ‰(𝑦)β†’0𝑉𝑔=𝛾0𝐺𝑒Upξ€Έ+π‘ˆπ‘‡ξ€·π‘’Upβˆ’π‘’π‘ξ€Έ,(19)𝑒Low and π‘’π‘ˆ are defined by implicit relations to the gate voltage 𝑉𝑔. Following a methodology previously described in [17], 𝑒Low and 𝑒Up can easily be obtained iteratively from first order Taylor expansion by setting 𝑉𝑔=π‘˜β‹…π›Ώ, in which 𝛿 is the sample step size and π‘˜ an integer. For a given step size and for π‘˜ varying from 0 to 𝑉𝑔/𝛿, 𝑒Low and 𝑒Up at iteration π‘˜ are calculated values at previous iteration according to 𝑒Low,π‘˜+1=𝑒Low,π‘˜+𝛿𝑑𝑒||||𝑑𝑉𝑔𝑒=𝑒Low,π‘˜,(20)𝑑𝑒||||𝑑𝑉𝑔𝑒=𝑒Low,π‘˜=1π‘ˆπ‘‡ξƒ¬π›Ύ1+0π‘ˆπ‘‡π‘’βˆ’π‘’π‘βˆ’π‘’π‘’π‘βˆ’π‘’βˆ’π‘’Low,π‘˜+𝑒𝑒Low,π‘˜2√𝐻(𝑒Low,π‘˜)+𝐺(𝑒Low,π‘˜)ξƒ­βˆ’1,(21)𝑒Up,π‘˜+1=𝑒Up,π‘˜+𝛿𝑑𝑒||||𝑑𝑉𝑔𝑒=𝑒Up,π‘˜,(22)𝑑𝑒||||𝑑𝑉𝑔𝑒=𝑒Up,π‘˜=1π‘ˆπ‘‡ξƒ¬π›Ύ1+0π‘ˆπ‘‡π‘’βˆ’π‘’π‘βˆ’π‘’βˆ’π‘’Up,π‘˜2√𝐺(𝑒Up,π‘˜)ξƒ­βˆ’1,(23)

In equations (20) and (22), the initial value π‘˜=0 corresponds to flat band (𝑉𝑔=0), and the starting conditions are 𝑒Low,0=𝑒Up,0=𝑒𝑏.

In (20), 𝑒Low,π‘˜ reaches the threshod voltage 𝑉𝑇 when: π‘˜=π‘˜π‘‡=int(𝑉𝑇/𝛿); with the function int = (number down to the nearest integer). As a result, 𝑉𝑇=1.2819 if 𝑁𝐴=1017cmβˆ’3 and 𝑑ox=10 nm and π‘˜π‘‡=12819.  𝑒Low,π‘˜π‘‡=15.712=βˆ’π‘’π‘ when the step is 𝛿=10βˆ’4 and the corresponding band bending is Ψ𝑆=π‘ˆπ‘‡(𝑒Low,π‘˜π‘‡βˆ’π‘’π‘)=βˆ’2πœ‘π‘.

As this has been previously shown [12], the error induced by first order Taylor expansion depends on the step size 𝛿. Equations (20) and (22) show that the relative error is in the same order of magnitude as the step 𝛿. For instance, a step 𝛿=10βˆ’10 gives a relative error less than 10βˆ’10.

Figure 2 shows 𝑒Low and 𝑒Up plots  versus 𝑉𝑔. The surface potential has well defined upper and lower limits in strong inversion, however, those limits are almost confounded in weak inversion (small channel voltage drop), which might give issues when 𝑒Low and 𝑒Up are used as integration limits in Pao-Sah integral (see Section 7).

4.2. The Inversion of πœ‰(𝑦)=𝑓[𝑒𝑆(𝑦)]

Equation (14) is an explicit relationship between 𝑒𝑆(𝑦) and πœ‰(𝑦). As this has been done in previous section with the calculation of the lower (𝑒Low) and upper (𝑒Up) limits of the surface potential, the inversion of (14) can be obtained by a first order Taylor expansion of the inverse function 𝑒𝑆(𝑦)=π‘“βˆ’1[πœ‰(𝑦)] by setting πœ‰(𝑦,π‘š)=π‘šβ‹…β„Ž in which β„Ž is the sample step size and π‘š an integer:𝑒𝑆,π‘š+1=𝑒𝑆,π‘š+β„Žπ‘‘π‘’π‘†||||π‘‘πœ‰π‘’π‘†=𝑒𝑆,π‘š,(24)𝑑𝑒𝑆||||π‘‘πœ‰π‘’π‘†=𝑒𝑆,π‘š=𝑒𝑒𝑆,π‘šβˆ’π‘’π‘’π‘π»ξ€·π‘’π‘†,π‘šξ€Έ+2π‘ˆπ‘‡ξ€½π‘‰π‘”βˆ’π‘ˆπ‘‡ξ€Ίπ‘’π‘†,π‘šβˆ’π‘’π‘ξ€»ξ€Ύ/𝛾20+π‘’βˆ’π‘’π‘βˆ’π‘’βˆ’π‘’π‘†,π‘šπΈξ€Ίπ‘’π‘†,π‘šξ€»ξ€Ίπ‘’βˆ’πΊπ‘†,π‘šξ€»ξƒ°βˆ’1.(25)

In (24), initial value for 𝑒𝑆 is 𝑒𝑆,0=𝑒Low,π‘˜; 𝑒Low,π‘˜ from (20) with π‘˜=𝑉𝑔/π‘ˆπ‘‡, and iteration stops at π‘š=𝑀=𝑉(𝑦)/β„Žπ‘ˆπ‘‡.

Figure 3 shows [𝑒𝑆(𝑦),𝑉(𝑦)] plots along the channel in strong inversion (𝑉𝑔=5 V). (a) is derived from (14) by sampling 𝑒𝑆(𝑦) and (b) is derived from iterative (24) (step β„Ž=0.01). Curves (a) and (b) merge with an error less than 10βˆ’4.

Note that this representation is 𝑒𝑆(𝑦)  versus 𝑉(𝑦) and not 𝑒𝑆(𝑦)=𝑔(𝑦) which would need the knowledge of the variation of 𝑉(𝑦)  versus 𝑦. Such variations suppose a complete 2-D resolution of the Poisson equation, but as is discussed later in Section 4, under the gradual channel approximation, a precise knowledge of 𝑉(𝑦)  versus 𝑦 is not needed for generating the main device current voltage and other related characteristics.

Figure 4 shows a complete set of curves 𝑒𝑆(𝑦)  versus 𝑉(𝑦) from 𝑉𝑔=1.2 up to 5 V. 𝑒𝑆(𝑦)=𝑒𝑆,π‘š is calculated from (24) and πœ‰(𝑦) from 𝑒𝑆,π‘š in (14). This figure illustrates the pinch-off voltage in strong inversion which appears when 𝑒𝑆(𝑦)βˆ’πœ‰(𝑦) becomes to decrease.

In terms of surface electron density 𝑛𝑆(𝑦)=𝑛𝑖𝑒𝑒𝑆(𝑦)βˆ’πœ‰(𝑦) (Figure 5), the difference between strong and weak inversion is evident. The drift current dominates when surface electron density is almost constant and diffusion takes place when there is a 𝑑𝑛𝑆/𝑑𝑦 gradient.

In strong inversion, the transition between 𝑛𝑆(𝑦)=𝑐𝑑𝑒 and the 𝑛(𝑦)π›Όπ‘’βˆ’πœ‰(𝑦) regimes happens for a voltage 𝑉𝑑(𝑉𝑔). The inset plots in Figure 5 shows that 𝑉𝑑(𝑉𝑔)=0 when 𝑉𝑔 is equal to the threshold voltage 𝑉𝑇=βˆ’2πœ‘π‘βˆš+𝛾|2πœ‘π‘|/π‘ˆπ‘‡β‰ˆ1.28 V defined in the charge sheet model [18].

5. The Pao-Sah Double Integral

Under the gradual channel approximation [13] the drift drain-source current density varies from bulk silicon toward gate oxide-silicon interface and along the channel:𝐽(π‘₯,𝑦)=π‘žπœ‡π‘›π‘›(π‘₯,𝑦)𝐹(π‘₯,𝑦).(26)

Or, in terms of potential𝐽[]𝑒(π‘₯,𝑦),πœ‰(𝑦)=π‘žπœ‡π‘›π‘›π‘–π‘’[𝑒(π‘₯,𝑦)βˆ’πœ‰(𝑦)]𝐹[]𝑒(π‘₯,𝑦),πœ‰(𝑦).(27)

With the introduction of the inversion charge density𝑄invξ€œ=π‘žπ‘₯𝑑0𝑛(π‘₯,𝑦)𝑑π‘₯=π‘žπ‘›π‘–π‘ˆπ‘‡ξ€œπ‘’π‘†0(𝑦)π‘’π‘’βˆ’πœ‰(𝑦)𝐹(𝑒,πœ‰(𝑦))𝑑𝑒,(28) where we recall 𝐹=βˆ’π‘ˆπ‘‡(𝑑𝑒/𝑑π‘₯) and π‘₯𝑑 the inversion length in π‘₯ direction; the expression for the drain-source drift current follows𝐼drift=πœ‡π‘›π‘ŠπΏξ€œπ‘‰π·0𝑄inv𝑑𝑉(𝑦)=πœ‡π‘›π‘ŠπΏπ‘ˆπ‘‡ξ€œ0πœ‰(𝐿)𝑄invπ‘‘πœ‰(𝑦),(29) in which a total channel width π‘Š is assumed. The Pao-Sah double integral then reads𝐼drift=π‘žπ‘›π‘–πœ‡π‘›π‘ŠπΏπ‘ˆ2π‘‡ξ€œ0πœ‰(𝐿)ξ‚»ξ€œπ‘’π‘†0π‘’π‘’βˆ’πœ‰(𝑦)𝐹(𝑒,πœ‰(𝑦))π‘‘π‘’π‘‘πœ‰(𝑦).(30)

Since πœ‰(𝑦) can be expressed in terms of 𝑒𝑆(𝑦) or alternatively 𝑒𝑆(𝑦) defined as function of πœ‰(𝑦), the Pao-Sah double integral can be reduced into separate single integrals depending on the choice we make for the integration variable.

5.1. Solution from Surface Potential 𝑒(𝑠)

πœ‰(𝑦) is a function of 𝑒𝑆(𝑦) from (14). In (30), π‘‘πœ‰(𝑦) is replaced byπ‘‘πœ‰(𝑦)=π‘‘πœ‰(𝑦)𝑑𝑒𝑆(𝑦)𝑑𝑒𝑆(𝑦),(31) in which π‘‘πœ‰(𝑦)/𝑑𝑒𝑆(𝑦) from (15), is function of 𝑒𝑆(𝑦).

The Pao-Sah double integral can be expressed in terms of two single iterated integrals [19], and the drift current is given by𝐼drift=π‘žπ‘›π‘–π‘ŠπΏπ‘ˆ2π‘‡Γ—ξ€œπ‘’π‘†π‘’(𝐿)𝑆(0)ξ‚»ξ€œπ‘’π‘†0(𝑦)πœ‡neο¬€π‘’π‘’βˆ’πœ‰(𝑦)𝐹(𝑒,πœ‰(𝑦))π‘‘π‘’π‘‘πœ‰(𝑦)𝑑𝑒𝑆(𝑦)𝑑𝑒𝑆(𝑦).(32)

The integral in braces is integrated first and gives a function of 𝑒𝑆(𝑦) and the final integral is integrated with respect to 𝑒𝑆(𝑦) from 𝑒𝑆(0) to 𝑒𝑆(𝐿).

In (32):(i)𝑒𝑆(𝑦) is the surface potential along the channel.(ii)πœ‰(𝑦) is a function of 𝑒𝑆(𝑦) from  (14).(iii)π‘‘πœ‰(𝑦)/𝑑𝑒𝑆(𝑦) is given by  (15).(iv)𝑒𝑆(0) is the surface potential in 𝑦=0. At a constant 𝑉𝑔 bias, it corresponds to 𝑒Low,π‘˜ solution of (20) with π‘˜=𝑉𝑔/𝛿.(v)𝑒𝑆(𝐿) is the surface potential in 𝑦=𝐿. At a constant 𝑉𝑔, it corresponds to 𝑒𝑆,π‘š solution of (24) with π‘š=𝑉𝐷/β„Žπ‘ˆπ‘‡.(vi)πœ‡neff is the effective mobility which has extensively been studied previously [20, 21]. The correction over a constant mobility (πœ‡π‘›=550cm2β‹…Vβˆ’1sβˆ’1) used in this paper, principally depends on 𝐹[𝑒,πœ‰(𝑦)] and is easy to implement in the integral.

5.2. Solution from Channel Potential πœ‰(𝑦)=𝑉(𝑦)/π‘ˆπ‘‡

In the drift current (30), it is possible to define the surface potential 𝑒𝑆  versus πœ‰ from (24). In this case, the channel potential is sampled according to Taylor expansion, and the integral on πœ‰=π‘šβ„Ž must be replaced by a discrete summation𝐼drift=π‘žπ‘›π‘–π‘ˆ2π‘‡π‘ŠπΏπ‘‰π·/β„Žπ‘ˆπ‘‡ξ“π‘š=0ξ‚»ξ€œπ‘’π‘†,π‘š0πœ‡neο¬€π‘’π‘’βˆ’π‘šβ„Žξ‚ΌπΉ(𝑒,πœ‰(𝑦))π‘‘π‘’β„Ž.(33)

This expression contains only one integral and is more suitable for numerical treatment with a simple worksheet; but this leads to a longer calculation time due to multiple loops in the summation.

5.3. Drift Current-Voltage Characteristics

Equation (32) has been calculated with a C-program in a large range of drain and gate voltages. The integration is made from Simpson algorithm.

Figure 6 shows the simulation results in log scale with 𝑁𝐴=1017cmβˆ’3 and 𝑑ox=10 nm. The solid lines correspond to (32) and markers (+) to (33). The comparison between (32) and (33), in weak and strong inversion gives numerical values with accuracy better than 1%.

Figures 7(a) and 7(b) show the simulation results in linear scale. Curve (a) is compared to data reported in [4] (𝑁𝐴=1015cmβˆ’3,𝑑ox=50nm) and curve (b) with more recent data reported in [7] (𝑁𝐴=5.1017cmβˆ’3,𝑑ox=5nm). All these calculations have been performed using a constant mobility: πœ‡π‘›=550cm2Vβˆ’1sβˆ’1.

6. Transconductance

The transconductance of the drift current is defined from the derivative [3]π‘”π‘š=𝑑𝐼drift||||𝑑𝑉𝑔𝑉𝐷=𝑑𝐼drift𝑑𝑒𝑆||||𝑉𝐷𝑑𝑒𝑆||||𝑑𝑉𝑔𝑉𝐷.(34)

The mathematical definition of the integral operator [22] gives the conditionξ€œIf𝑔(π‘₯)=π‘₯0𝑓(𝑒)𝑑𝑒,thenβŸΉπ‘‘π‘”(π‘₯)𝑑π‘₯=𝑓(π‘₯),(35) and, consequentlyπ‘‘π‘‘π‘’π‘†ξ‚»ξ€œπ‘’π‘†0π‘’π‘’βˆ’πœ‰ξ‚Ό=𝑒𝐹(𝑒,πœ‰)π‘‘π‘’π‘’π‘†βˆ’πœ‰πΉξ€·π‘’π‘†ξ€Έ,πœ‰.(36)

The transconductance is calculated at a constant 𝑉𝐷 bias. The derivative of the drain current is only on 𝑒𝑆 and (36) in (32) results in𝑑𝐼drift𝑑𝑒𝑆=π‘žπ‘›π‘–π‘ˆ2π‘‡π‘ŠπΏξ€œπ‘‰π·/π‘ˆπ‘‡0ξƒ―πœ‡neο¬€π‘’π‘’π‘†βˆ’πœ‰πΉξ€·π‘’π‘†ξ€Έξƒ°,πœ‰π‘‘πœ‰.(37)

The transconductance isπ‘”π‘š=π‘žπ‘›π‘–π‘ˆ2π‘‡π‘ŠπΏξ€œπ‘‰π·/π‘ˆπ‘‡0πœ‡neο¬€π‘’π‘’π‘†βˆ’πœ‰πΉξ€·π‘’π‘†ξ€Έ,πœ‰π‘‘π‘’π‘†π‘‘π‘‰π‘”π‘‘πœ‰,(38) which from (31) gives a single integral of the surface potential 𝑒𝑆(𝑦)π‘”π‘š=π‘žπ‘›π‘–π‘ˆ2π‘‡π‘ŠπΏΓ—ξ€œπ‘’π‘†π‘’(𝐿)𝑆(0)πœ‡neο¬€π‘’π‘’π‘†βˆ’πœ‰πΉξ€·π‘’π‘†ξ€Έξ‚΅,πœ‰π‘‘π‘’π‘†(𝑦)π‘‘π‘‰π‘”ξ‚Άξ‚΅π‘‘πœ‰(𝑦)𝑑𝑒𝑆(𝑦)𝑑𝑒𝑆(𝑦),(39) with𝑑𝑒𝑆(𝑦)=1π‘‘π‘‰π‘”π‘ˆπ‘‡βŽ‘βŽ’βŽ’βŽ£π›Ύ1+0π‘ˆπ‘‡π‘’βˆ’πœ‰π‘‘π»/𝑑𝑒+𝑑𝐺/𝑑𝑒2βˆšπ‘’βˆ’πœ‰|||||𝐻+𝐺𝑒=𝑒𝑆(𝑦)⎀βŽ₯βŽ₯βŽ¦βˆ’1.(40)

Figure 8(a) shows [π‘”π‘š,𝑉𝑔] plots in log scale when the source is grounded (𝑉𝑆=0,𝑉𝐷=5 V) and Figure 8(b) shows the simulation results in linear scale using (39) compared to the charge sheet model in the quadratic region when π‘”π‘š=πœ‡π‘›πΆ0𝑉𝐷(π‘Š/𝐿) is a linear function  versus 𝑉𝐷.

In the E.K.V. model [23], the authors introduce the normalized transconductance: 𝐺(𝑖𝐹)=π‘”π‘šπ‘›π‘ˆπ‘‡/𝑖𝐹 with 𝑖𝐹=𝐼𝐹/2π‘›π›½π‘ˆ2𝑇; 𝑛 is the slope factor and 𝛽=πœ‡π‘›πΆoxπ‘Š/𝐿. Figure 8(c) shows [𝐺(𝑖𝐹),𝑖𝐹] plots using (39) and (32) (𝐼𝐹=𝐼drift). The function 𝐺(𝑖𝐹), the asymtotes 𝐺(𝑖𝐹)=1 and 𝐺(𝑖𝐹)=(𝑖𝐹)βˆ’1/2, respectively, in weak and strong inversion are in good agreement with this model.

7. The Accuracy of Pao-Sah Integral

We have integrated the Pao-Sah equation between the limits 𝑒𝑆(0) and 𝑒𝑆(𝐿) given by iterative equations (20) and (24). Figure 9 shows the variations of these integration limits  versus 𝑉𝐷, respectively, in strong (a) and weak inversion (b). The difference between 𝑒𝑆(0) and 𝑒𝑆(𝐿) is sufficient large in strong inversion. However, the situation is not the same in weak inversion. In this regime, the Pao-Sah integral will be integrated on a very small range and the accuracy of the simulation heavily depends on the number of digits used in the floating point arithmetic.

The accuracy of Pao-Sah equation can be increased if the step 𝛿 in (20) is decreased. As an example, in weak inversion, 𝛿=10βˆ’10 gives an accurate value of 𝑒Low,π‘˜ which leads to 𝑉𝑔[𝑒Low,π‘˜] equal to 𝑉𝑔 with an absolute error less than 10βˆ’10.

8. The Diffusion Current

8.1. Integral of Current Diffusion

The diffusion current density of an n-MOSFET is given by𝐽diff(π‘₯,𝑦)=βˆ’πœ‡π‘›πΎπ‘‡π‘‘π‘›(π‘₯,𝑦)𝑑𝑦.(41)

As shown in Figure 1, the gradient of concentration gives a diffusion of electrons moving from source to drain. The diffusion current is in the same direction as the drift current. For the same reason, it is evident that the diffusion current will have a relative contribution higher in weak inversion than in strong inversion. For instance, in strong inversion, a large region of the channel has a constant electrons concentration which does not contribute to the diffusion current.

The total diffusion current is the integral of the diffusion current density𝐼diο¬€ξ€œ=π‘Š0π‘₯𝑑𝐽diff(π‘₯,𝑦)𝑑π‘₯=βˆ’πœ‡π‘›π‘‘πΎπ‘‡π‘Šξ‚»ξ€œπ‘‘π‘¦0π‘₯𝑑𝑛(π‘₯,𝑦)𝑑π‘₯.(42)

Equation (28) gives the equivalence to potentialsξ€œ0π‘₯𝑑𝑛(π‘₯,𝑦)𝑑π‘₯=π‘›π‘–π‘ˆπ‘‡ξ€œ0𝑒(𝑠)π‘’π‘’βˆ’πœ‰πΉ(𝑒,πœ‰)𝑑𝑒,(43) and after integration on 𝑦, (42) becomes𝐼diff=βˆ’π‘žπœ‡π‘›π‘ŠπΏπ‘›π‘–π‘ˆ2π‘‡ξ€œπΏ0π‘‘ξ‚»ξ€œπ‘’π‘†0(𝑦)π‘’π‘’βˆ’πœ‰ξ‚ΌπΉ(𝑒,πœ‰)𝑑𝑒,(44) and reduces to single integrals𝐼diff=π‘žπœ‡π‘›π‘ŠπΏπ‘›π‘–π‘ˆ2π‘‡Γ—ξ‚»ξ€œπ‘’π‘†0(0)π‘’π‘’βˆ’πœ‰(0)𝐹[]ξ€œπ‘’,πœ‰(0)π‘‘π‘’βˆ’π‘’π‘†0(𝐿)π‘’π‘’βˆ’πœ‰(𝐿)𝐹[]ξ‚Ό.𝑒,πœ‰(𝐿)𝑑𝑒(45)

Equation (45) is equivalent to the well known expression𝐼diff=πœ‡π‘›π‘ŠπΏπ‘ˆπ‘‡ξ€Ίπ‘„inv(0)βˆ’π‘„invξ€»(𝐿).(46)

This expression is easily calculated from the exponential variations of 𝑄inv  versus 𝑉𝑔 in weak inversion. Equation (13) gives the exact relation 𝑉𝑔=𝑓[𝑒𝑆(𝑦),πœ‰(𝑦)], and the inversion charge 𝑄inv  versus 𝑒𝑆(𝑦) is given from (28).

Figure 10 shows [𝑄inv,𝑉𝑔] plots according to (13) and (28) for different 𝑉𝐷 bias. 𝑄inv(0) corresponds to πœ‰(𝑦)=0 and 𝑄inv(𝐿) corresponds to the drain bias πœ‰(𝐿)=𝑉𝐷/π‘ˆπ‘‡. 𝑄inv(𝐿) is rapidly negligible in (45). Figure 11 shows 𝐼diff(𝑉𝑔) plots in weak inversion calculated from (45) and compared with (46) from 𝑄inv in (28).

In weak inversion, the diffusion current is well captured by𝐽diff=𝐽0𝑒𝑉𝑔/𝑁𝐾𝑇,(47) with 𝑁=1.38, which is in good agreement with the slope of [𝑄inv,𝑉𝑔] plots (Figure 10).

8.2. Diffusion Plots 𝐼𝑑𝑖𝑓𝑓(𝑉𝐷)

A complete graph 𝐼diff(𝑉𝐷) can be plotted by sampling 𝑒𝑆(𝐿) from (24). Figure 12 shows the contribution of the diffusion current (𝐼diff) compared with the drift current (𝐼𝑑) in strong and weak inversion (𝐼𝑑1>𝐼diff1 and 𝐼𝑑2<𝐼diff2). The diffusion current is negligible in strong inversion and dominant in weak inversion.

9. Effects of Source Bias

In Section 5, the Pao-Sah integral is calculated with a zero bias source. If the source is biased with a voltage 𝑉𝑆𝐡  versus bulk silicon, the reduced drain-source potential πœ‰(𝑦) along the channel varies from πœ‰(0)=𝑉𝑆𝐡/π‘ˆπ‘‡ to πœ‰(𝐿)=𝑉𝐷/π‘ˆπ‘‡. Then, the surface potential 𝑒𝑆,π‘š starts at 𝑒𝑆,π‘šπ‘† given from (24) with π‘šπ‘†=𝑉𝑆𝐡/β„Žπ‘ˆπ‘‡ and stops at π‘š=𝑀=𝑉𝐷𝐡/β„Žπ‘ˆπ‘‡.

The drift and the diffusion currents are always given by (32) and (45) by substituting 𝑒𝑆,0 by 𝑒𝑆,π‘šπ‘† solution of (24) with π‘šπ‘†=𝑉𝑆𝐡/β„Žπ‘ˆπ‘‡.

By setting𝐺𝑒𝑆=π‘žπ‘›π‘–π‘ŠπΏπ‘ˆ2π‘‡ξ€œπ‘’π‘†0πœ‡neο¬€π‘’π‘’βˆ’πœ‰πΉ(𝑒,πœ‰)𝑑𝑒.(48)

Equation (32) can be rewritten with π‘š index in 𝑒𝑆𝐼𝑑=ξ€œπ‘’π‘†,𝑀𝑒𝑆,0𝐺[]𝑒(𝑠)π‘‘πœ‰π‘‘π‘’π‘†π‘‘π‘’π‘†βˆ’ξ€œπ‘’π‘†,π‘šπ‘ π‘’π‘†,0𝐺[]𝑒(𝑠)π‘‘πœ‰π‘‘π‘’π‘†π‘‘π‘’π‘†,𝐼(49)𝑑=𝐼𝐹+𝐼𝑅.(50)

The integrals in (49) correspond, respectively, to the forward 𝐼𝐹 and reverse 𝐼𝑅 currents introduced in the E.K.V model [23]. The same expressions can also be defined in the diffusion current in (45).

10. Transfer Characteristics

The transfer characteristics represented in Figure 13 shows 𝐼=𝐼drift+𝐼diff=𝑓(𝑉𝑔) calculated from (32) and (45) for constant 𝑉𝐷𝑆 and different 𝑉𝑆𝐡. The parameters of 𝐼=𝑓(𝑉𝑔) are 𝑉FB(flatband)=0, 𝑁𝐴=1017cmβˆ’3,𝑑ox=10 nm, 𝑉𝐷𝑆=5 volt and 𝑉𝑆𝐡 from 0 to 1 volt.

As a comparison, Figure 14 shows the curves 𝐼=𝑓(𝑉𝑔) with the same parameters reported in [4]. The conditions are 𝑉FB(flatband)=0.92 for 𝑁𝐴=1014cmβˆ’3and 𝑉FB=0.86 for 𝑁𝐴=1015cmβˆ’3. The oxide thickness is 𝑑ox=13 nm and the drain bias is 𝑉𝐷=1 volt. Our simulations are in good agreement with these results.

11. Conclusion

In previous papers, we presented an analytic resolution of Poisson-Boltzmann equation applicable in semiconductor junctions [12, 24]. We showed that the analytic method, which can be lead as far as possible without approximation (except the hypothesis of Channel Gradual Approximation), is able to illustrate the influence of physical parameters in surface potential and carriers density. By following the same methodology, the drift and diffusion current and the transconductance in MOSFET are given by iterated integrals easily solved quite instantaneously. The iterative treatment by Taylor expansions leads to a reasonable computation efficiency and simulation speed. All simulations using the flowchart as described in Section 5.1 are almost instantaneous with a simple C-program.

The excellent agreement of our results with the standard models [8], can be considered as an accurate tools for users in the complete knowledge of the MOSFET without access to specific CAD software. Moreover, we are well aware that the present work can not be a complete model of the MOSFET in terms of equivalent circuit with resistances and capacitances. But, by a simple calculation, available to a large community, it could give an overview of the complete MOSFET in all inversion modes with a single integral formulation.

Figure 15 shows the user interface for current calculation. The parameters are doping 𝑁𝐴, oxide thickness 𝑑ox, and gate and drain voltages 𝑉𝑔 and 𝑉𝐷.

12. Annex: Taylor Polynomials of Inverse Functions

If a function 𝑦=𝑓(π‘₯) has continuous derivatives up to (𝑛)th order 𝑓(𝑛)(π‘₯), then this function can be expanded in the following Taylor polynomials [25]ξ€·π‘₯𝑓(π‘₯)=𝑓0ξ€Έ+𝑛1ξ€·π‘₯βˆ’π‘₯0𝑛𝑓𝑛!(𝑛)(π‘₯)+𝑅(𝑛),(51) where 𝑅(𝑛) is called the remainder after 𝑛+1 terms.

This expansion converges over a certain range of π‘₯, if limπ‘›β†’βˆžπ‘…(𝑛)=0, and the expansion is called the Taylor Series of 𝑓(π‘₯) expanded about π‘₯0.

In inverse function, if 𝑦=𝑓(π‘₯) is a one-to-one function, then π‘“βˆ’1 is continuous and, if π‘“βˆ’1 has continuous derivatives up to (𝑛)th order, π‘“βˆ’1 can be expanded in the Taylor polynomials. These conditions have been previously verified in the inversion of surface potential [12], with an accuracy compatible with the integral simulations.

The limits of Paoh-Sah integral are based on the inverse function 𝑒𝑆(𝑦)=π‘“βˆ’1[πœ‰(𝑦)], The first derivative variations given by (25) are shown in Figure 16 for different 𝑉𝑔 bias. The first derivative is define in all the [0,𝑉𝐷] region and varies from 1 to β‰ˆ10βˆ’13 without discontinuity. The strong decrease observed for a 𝑉𝐷 value depending on 𝑉𝑔 is due to the asymptotic value in the surface potential 𝑒𝑆(𝑦) when 𝑒𝑆(𝑦) reaches 𝑒Up, solution of 𝐸(𝑒)βˆ’πΊ(𝑒)=0.

Nomenclature

𝑉𝑔: Voltage gate with bulk silicon groundedπœ€π‘† and πœ€ox: Silicon and silicon oxide permittivity𝑑ox: Oxide thickness𝐿: Channel lengthπ‘Š: Channel width𝑛𝑖: Intrinsic carrier concentration in cmβˆ’3𝑁𝐴 and 𝑁𝐷: Dopant concentrations in cmβˆ’3π‘ˆπ‘‡=𝐾𝑇/π‘ž: Thermal voltageπœ‘(𝑏)=πœ‘π‘=βˆ’π‘ˆπ‘‡ln(𝑁𝐴/𝑛𝑖): Bulk potential of 𝑝-doped silicon𝑒(π‘₯)=πœ‘(π‘₯)/π‘ˆπ‘‡: Reduced potential𝑉𝑇=βˆ’2πœ‘π‘βˆš+𝛾|2πœ‘π‘|/π‘ˆπ‘‡: Threshold voltage in strong inversion𝛾0=(1/𝐢ox)√2πΎπ‘‡πœ€π‘†π‘›π‘–: Intrinsic body factor with 𝐢ox=𝑑ox/πœ€ox𝛾=(1/𝐢ox)√2πΎπ‘‡πœ€π‘†π‘π΄: 𝑝-type semiconductor body factor.

Numerical applications use SI units, excepted dopant concentration in cmβˆ’3, πœ€π‘† and πœ€ox in farads.cmβˆ’1.