# Taylor Expansion of Surface Potential in MOSFET: Application to Pao-Sah Integral

Authors: Hugues Murray; Patrick Martin; Serge Bardy
Journal: Active and Passive Electronic Components
 (2010)

---

## 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.


---

## Body

## 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 Vg is above a threshold voltage VT [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φ(x) 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 F(x). The gate voltage Vg  versus surface potential φS 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 V(y) 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 y direction parallel to the conducting channel is assumed to be much smaller than in the vertical x direction toward silicon bulk allowing the formulation divF=dF/dx and F[x,y]=-d[φ(x)]/dx. 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: ρ(x,y)=q[p(x,y)-n(x,y)+ND-NA].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 Vg is defined relative to flat band.Figure 1
The channel scheme in the Gradual Channel Approximation.The presentation relates to ap-type semiconductor (n-MOSFET), but this is not restrictive and could easily be extended to n-type semiconductor. Upon usual notation conventions (see Nomenclature section), the bulk Fermi level is φb. φ(x,y) is the potential in every points of the semiconductor material and V(y) is the channel voltage which defines the quasi-Fermi level offset from equilibrium for electrons. Upon application of a drain bias VD, V(y) ranges from 0 at the source side (y=0) to VD at the drain side (y=L).UT=KT/q is the thermal voltage.In the following, we use the reduced potentialsub=φb/UT,u(x,y)=φ(x,y)/UT and ξ(y)=V(y)/UT instead of φb,φ(x,y) and V(y), respectively. uS(y)=u(x,y)|x=0,FS(y)=F(x,y)|x=0,nS(y)=n(x,y)|x=0 are shorthand notations for surface quantities.In the first MOSFET modeling, charge densities have been written by introducing the quasi-Fermi level only inn(x,y) with the expressions:(1)n(x,y)=nieu(x,y)-ξ(y),(2)p(x)=nie-u(x,y),(3)NA=nie-ub,(4)ND=nieub.After solving the 1-D Poisson equationdF(x,y)/dx=ρ(x,y), these expressions give the 1-D electric field by(5)F(x,y)=2KTniεS×e-ξ(y){eu(x,y)-eub}+e-u(x,y)-[e-ub-eub][ub-u(x,y)]-e-ub,
and its derivative(6)dF(x,y)dξ(y)=-2KTniεSe-ξ(y)[eu(x,y)-eub]2F(x,y).Aseub≪eu(x,y) in the channel region, this expression is equivalent to(7)dF(x,y)dξ(y)=-2KTεSn(x,y)2F(x,y),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 n(x→∞)=ND, (4) must be corrected by(8)ND=nieub-ξ(y).With this exact formulation of the donor density, the electric field gradient from the 1-D Poisson equation (5) now reads(9)dF(x)dx=qniεS(e-u(x,y)-eu(x,y)-ξ(y)+eub-ξ(y)-e-ub),
which gives the electric field(10)F(x,y)=2KTniεSe-ξ(y)H[u(x,y)]+G[u(x,y)],
where we have defined(11)H(u)=eu+eub[ub-u-1],G(u)=e-u-e-ub[ub-u+1].In this form, the derivative(12)dF(x,y)dξ(y)=-2KTεSe-ξ(y)H[u(x,y)]2F(x,y),
is no longer with a numerator proportional to n(x,y) 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 FS(y)=F(x,y)|x=0 is given by substituting u and u(x,y) by uS(y) in (10)-(11).


## 3. The Surface Potential Dependence to Gate Bias

From electrostatic considerations [14], the gate voltage is expressed by(13)Vg=γ0e-ξ(y)H[uS(y)]+G[uS(y)]+UT[uS(y)-u(b)].The exponential expression in (13) allows to separate u(s) and ξ(y)(14)ξ(y)=lnH[uS(y)]E[uS(y)]-G[uS(y)],
with the derivative(15)dξ(y)du(s)=dH[uS(y)]/duS(y)H[uS(y)]-dE[uS(y)]/duS(y)-dG[uS(y)]/duS(y)E[uS(y)]-G[uS(y)],
in which we introduce the dimensionless quantity(16)E(u)={Vg-UT[u-ub]γ0}2.Several approximations to (13) have been reported. For instance, the equation(17)Vg2=Ψ(s)+γ(niNA)2exp[Ψ(s)-V(y)UT]+Ψ(s)UT,
is often used in practice in a large variety of surface potential-based models with the band-bending: Ψ(s)=φ(s)-φb. It provides a sufficiently accurate solution for the surface potential in inversion mode but generally lacks accuracy near the flat-band conditions: u(s)=ub.


## 4. The Surface Potential Dependence to Drain Bias

The explicit relationξ(y)=f[uS(y)] 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 ξ(y)=f[uS(y)] and to the inverse function uS(y)=f-1[ξ(y)] 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 uS(y) compatible with the denominator E[uS(y)]-G[uS(y)].


### 4.1. Boundary Limits of Surface Potential at Constant Gate Voltage

Equation (14) has only physical meaning when E[uS(y)]-G[uS(y)]>0. At a constant Vg value, the reduced surface potential uS(y) ranges between uLow (lower value), solution of (13) with ξ(y)=0(18)Vg=γ0H(uLow)+G(uLow)+UT(uLow-ub),
and between uUp (upper value), solution of (13) with ξ(y)→∞ leading to e-ξ(y)→0(19)Vg=γ0G(uUp)+UT(uUp-ub),uLow and uU are defined by implicit relations to the gate voltage Vg. Following a methodology previously described in [17], uLow and uUp can easily be obtained iteratively from first order Taylor expansion by setting Vg=k·δ, in which δ is the sample step size and k an integer. For a given step size and for k varying from 0 to Vg/δ, uLow and uUp at iteration k are calculated values at previous iteration according to (20)uLow,k+1=uLow,k+δdudVg|u=uLow,k,(21)dudVg|u=uLow,k=1UT[1+γ0UTe-ub-eub-e-uLow,k+euLow,k2H(uLow,k)+G(uLow,k)]-1,(22)uUp,k+1=uUp,k+δdudVg|u=uUp,k,(23)dudVg|u=uUp,k=1UT[1+γ0UTe-ub-e-uUp,k2G(uUp,k)]-1,In equations (20) and (22), the initial value k=0 corresponds to flat band (Vg=0), and the starting conditions are uLow,0=uUp,0=ub.In (20), uLow,k reaches the threshod voltage VT when: k=kT=int(VT/δ); with the function int =(number down to the nearest integer). As a result, VT=1.2819 if NA=1017cm-3 and tox=10 nm and kT=12819.  uLow,kT=15.712=-ub when the step is δ=10-4 and the corresponding band bending is ΨS=UT(uLow,kT-ub)=-2φb.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.Figure2 shows uLow and uUp plots  versus Vg. 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 uLow and uUp are used as integration limits in Pao-Sah integral (see Section 7).Figure 2
UpperuUp and lower uLow limits of the reduced surface potential  versus gate voltage Vg. (a) equation (20) and (b) equation (22). NA=1017cm-3 and tox=10 nm. Vg=VT when (a) crosses the line uLow=-ub.


### 4.2. The Inversion ofξ(y)=f[uS(y)]

Equation (14) is an explicit relationship between uS(y) and ξ(y). As this has been done in previous section with the calculation of the lower (uLow) and upper (uUp) limits of the surface potential, the inversion of (14) can be obtained by a first order Taylor expansion of the inverse function uS(y)=f-1[ξ(y)] by setting ξ(y,m)=m·h in which h is the sample step size and m an integer:(24)uS,m+1=uS,m+hduSdξ|uS=uS,m,(25)duSdξ|uS=uS,m={euS,m-eubH(uS,m)+2UT{Vg-UT[uS,m-ub]}/γ02+e-ub-e-uS,mE[uS,m]-G[uS,m]}-1.In (24), initial value for uS is uS,0=uLow,k; uLow,k from (20) with k=Vg/UT, and iteration stops at m=M=V(y)/hUT.Figure3 shows [uS(y),V(y)] plots along the channel in strong inversion (Vg=5 V). (a) is derived from (14) by sampling uS(y) and (b) is derived from iterative (24) (step h=0.01). Curves (a) and (b) merge with an error less than 10-4.Figure 3
[uS(y),V(y)] plots along the channel. (a) equation (14) and (b) equation (24) (step h=0.01). NA=1017cm-3,tox=10 nm and Vg=5 V.Note that this representation isuS(y)  versus V(y) and not uS(y)=g(y) which would need the knowledge of the variation of V(y)  versus y. 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 V(y)  versus y is not needed for generating the main device current voltage and other related characteristics.Figure4 shows a complete set of curves uS(y)  versus V(y) from Vg=1.2 up to 5 V. uS(y)=uS,m is calculated from (24) and ξ(y) from uS,m in (14). This figure illustrates the pinch-off voltage in strong inversion which appears when uS(y)-ξ(y) becomes to decrease.Figure 4
[uS(y),V(y)] and [uS(y)-ξ(y),V(y)] plots for different Vg bias. NA=1017cm-3 and tox=10 nm.In terms of surface electron densitynS(y)=nieuS(y)-ξ(y) (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 dnS/dy gradient.Figure 5
The electron density[nS(y),V(y)] plots along the channel from (1) and (14). NA=1017cm-3 and tox=10 nm.In strong inversion, the transition betweennS(y)=cte and the n(y)αe-ξ(y) regimes happens for a voltage Vt(Vg). The inset plots in Figure 5 shows that Vt(Vg)=0 when Vg is equal to the threshold voltage VT=-2φb+γ|2φb|/UT≈1.28 V defined in the charge sheet model [18].


## 4.1. Boundary Limits of Surface Potential at Constant Gate Voltage

Equation (14) has only physical meaning when E[uS(y)]-G[uS(y)]>0. At a constant Vg value, the reduced surface potential uS(y) ranges between uLow (lower value), solution of (13) with ξ(y)=0(18)Vg=γ0H(uLow)+G(uLow)+UT(uLow-ub),
and between uUp (upper value), solution of (13) with ξ(y)→∞ leading to e-ξ(y)→0(19)Vg=γ0G(uUp)+UT(uUp-ub),uLow and uU are defined by implicit relations to the gate voltage Vg. Following a methodology previously described in [17], uLow and uUp can easily be obtained iteratively from first order Taylor expansion by setting Vg=k·δ, in which δ is the sample step size and k an integer. For a given step size and for k varying from 0 to Vg/δ, uLow and uUp at iteration k are calculated values at previous iteration according to (20)uLow,k+1=uLow,k+δdudVg|u=uLow,k,(21)dudVg|u=uLow,k=1UT[1+γ0UTe-ub-eub-e-uLow,k+euLow,k2H(uLow,k)+G(uLow,k)]-1,(22)uUp,k+1=uUp,k+δdudVg|u=uUp,k,(23)dudVg|u=uUp,k=1UT[1+γ0UTe-ub-e-uUp,k2G(uUp,k)]-1,In equations (20) and (22), the initial value k=0 corresponds to flat band (Vg=0), and the starting conditions are uLow,0=uUp,0=ub.In (20), uLow,k reaches the threshod voltage VT when: k=kT=int(VT/δ); with the function int =(number down to the nearest integer). As a result, VT=1.2819 if NA=1017cm-3 and tox=10 nm and kT=12819.  uLow,kT=15.712=-ub when the step is δ=10-4 and the corresponding band bending is ΨS=UT(uLow,kT-ub)=-2φb.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.Figure2 shows uLow and uUp plots  versus Vg. 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 uLow and uUp are used as integration limits in Pao-Sah integral (see Section 7).Figure 2
UpperuUp and lower uLow limits of the reduced surface potential  versus gate voltage Vg. (a) equation (20) and (b) equation (22). NA=1017cm-3 and tox=10 nm. Vg=VT when (a) crosses the line uLow=-ub.


## 4.2. The Inversion ofξ(y)=f[uS(y)]

Equation (14) is an explicit relationship between uS(y) and ξ(y). As this has been done in previous section with the calculation of the lower (uLow) and upper (uUp) limits of the surface potential, the inversion of (14) can be obtained by a first order Taylor expansion of the inverse function uS(y)=f-1[ξ(y)] by setting ξ(y,m)=m·h in which h is the sample step size and m an integer:(24)uS,m+1=uS,m+hduSdξ|uS=uS,m,(25)duSdξ|uS=uS,m={euS,m-eubH(uS,m)+2UT{Vg-UT[uS,m-ub]}/γ02+e-ub-e-uS,mE[uS,m]-G[uS,m]}-1.In (24), initial value for uS is uS,0=uLow,k; uLow,k from (20) with k=Vg/UT, and iteration stops at m=M=V(y)/hUT.Figure3 shows [uS(y),V(y)] plots along the channel in strong inversion (Vg=5 V). (a) is derived from (14) by sampling uS(y) and (b) is derived from iterative (24) (step h=0.01). Curves (a) and (b) merge with an error less than 10-4.Figure 3
[uS(y),V(y)] plots along the channel. (a) equation (14) and (b) equation (24) (step h=0.01). NA=1017cm-3,tox=10 nm and Vg=5 V.Note that this representation isuS(y)  versus V(y) and not uS(y)=g(y) which would need the knowledge of the variation of V(y)  versus y. 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 V(y)  versus y is not needed for generating the main device current voltage and other related characteristics.Figure4 shows a complete set of curves uS(y)  versus V(y) from Vg=1.2 up to 5 V. uS(y)=uS,m is calculated from (24) and ξ(y) from uS,m in (14). This figure illustrates the pinch-off voltage in strong inversion which appears when uS(y)-ξ(y) becomes to decrease.Figure 4
[uS(y),V(y)] and [uS(y)-ξ(y),V(y)] plots for different Vg bias. NA=1017cm-3 and tox=10 nm.In terms of surface electron densitynS(y)=nieuS(y)-ξ(y) (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 dnS/dy gradient.Figure 5
The electron density[nS(y),V(y)] plots along the channel from (1) and (14). NA=1017cm-3 and tox=10 nm.In strong inversion, the transition betweennS(y)=cte and the n(y)αe-ξ(y) regimes happens for a voltage Vt(Vg). The inset plots in Figure 5 shows that Vt(Vg)=0 when Vg is equal to the threshold voltage VT=-2φb+γ|2φb|/UT≈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)J(x,y)=qμnn(x,y)F(x,y).Or, in terms of potential(27)J[u(x,y),ξ(y)]=qμnnie[u(x,y)-ξ(y)]F[u(x,y),ξ(y)].With the introduction of the inversion charge density(28)Qinv=q∫0xdn(x,y)dx=qniUT∫0uS(y)eu-ξ(y)F(u,ξ(y))du,
where we recall F=-UT(du/dx) and xd the inversion length in x direction; the expression for the drain-source drift current follows(29)Idrift=μnWL∫0VDQinvdV(y)=μnWLUT∫0ξ(L)Qinvdξ(y),
in which a total channel width W is assumed. The Pao-Sah double integral then reads(30)Idrift=qniμnWLUT2∫0ξ(L){∫0uSeu-ξ(y)F(u,ξ(y))du}dξ(y).Sinceξ(y) can be expressed in terms of uS(y) or alternatively uS(y) defined as function of ξ(y), 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 Potentialu(s)

ξ(y) is a function of uS(y) from (14). In (30), dξ(y) is replaced by(31)dξ(y)=dξ(y)duS(y)duS(y),
in which dξ(y)/duS(y) from (15), is function of uS(y).The Pao-Sah double integral can be expressed in terms of two singleiterated integrals [19], and the drift current is given by(32)Idrift=qniWLUT2×∫uS(0)uS(L){∫0uS(y)μneffeu-ξ(y)F(u,ξ(y))du}dξ(y)duS(y)duS(y).The integral in braces is integrated first and gives a function ofuS(y) and the final integral is integrated with respect to uS(y) from uS(0) to uS(L).In (32):(i)
uS(y) is the surface potential along the channel.(ii)
ξ(y) is a function of uS(y) from  (14).(iii)
dξ(y)/duS(y) is given by  (15).(iv)
uS(0) is the surface potential in y=0. At a constant Vg bias, it corresponds to uLow,k solution of (20) with k=Vg/δ.(v)
uS(L) is the surface potential in y=L. At a constant Vg, it corresponds to uS,m solution of (24) with m=VD/hUT.(vi)
μneff is the effective mobility which has extensively been studied previously [20, 21]. The correction over a constant mobility (μn=550cm2·V-1s-1) used in this paper, principally depends on F[u,ξ(y)] and is easy to implement in the integral.


### 5.2. Solution from Channel Potentialξ(y)=V(y)/UT

In the drift current (30), it is possible to define the surface potential uS  versus ξ from (24). In this case, the channel potential is sampled according to Taylor expansion, and the integral on ξ=mh must be replaced by a discrete summation(33)Idrift=qniUT2WL∑m=0VD/hUT{∫0uS,mμneffeu-mhF(u,ξ(y))du}h.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.Figure6 shows the simulation results in log scale with NA=1017cm-3 and tox=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%.Figure 6
Idrift  versus drain voltage VD from (32) in log scale. NA=1017cm-3,tox=10 nm and W/L=10. Markers (+) corresponds to (33).Figures7(a) and 7(b) show the simulation results in linear scale. Curve (a) is compared to data reported in [4] (NA=1015cm-3,tox=50nm) and curve (b) with more recent data reported in [7] (NA=5.1017cm-3,tox=5nm). All these calculations have been performed using a constant mobility: μn=550cm2V-1s-1.Idrift  versus drain voltage VD from (32) in linear scale compared with [4] and [7]. (a) NA=1015cm-3,tox=50 nm and W/L=10, (∘) corresponds to [4]. (b) NA=5.1017cm-3 and tox=5 nm, (∘) corresponds to [7].
(a)(b)


## 5.1. Solution from Surface Potentialu(s)

ξ(y) is a function of uS(y) from (14). In (30), dξ(y) is replaced by(31)dξ(y)=dξ(y)duS(y)duS(y),
in which dξ(y)/duS(y) from (15), is function of uS(y).The Pao-Sah double integral can be expressed in terms of two singleiterated integrals [19], and the drift current is given by(32)Idrift=qniWLUT2×∫uS(0)uS(L){∫0uS(y)μneffeu-ξ(y)F(u,ξ(y))du}dξ(y)duS(y)duS(y).The integral in braces is integrated first and gives a function ofuS(y) and the final integral is integrated with respect to uS(y) from uS(0) to uS(L).In (32):(i)
uS(y) is the surface potential along the channel.(ii)
ξ(y) is a function of uS(y) from  (14).(iii)
dξ(y)/duS(y) is given by  (15).(iv)
uS(0) is the surface potential in y=0. At a constant Vg bias, it corresponds to uLow,k solution of (20) with k=Vg/δ.(v)
uS(L) is the surface potential in y=L. At a constant Vg, it corresponds to uS,m solution of (24) with m=VD/hUT.(vi)
μneff is the effective mobility which has extensively been studied previously [20, 21]. The correction over a constant mobility (μn=550cm2·V-1s-1) used in this paper, principally depends on F[u,ξ(y)] and is easy to implement in the integral.


## 5.2. Solution from Channel Potentialξ(y)=V(y)/UT

In the drift current (30), it is possible to define the surface potential uS  versus ξ from (24). In this case, the channel potential is sampled according to Taylor expansion, and the integral on ξ=mh must be replaced by a discrete summation(33)Idrift=qniUT2WL∑m=0VD/hUT{∫0uS,mμneffeu-mhF(u,ξ(y))du}h.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.Figure6 shows the simulation results in log scale with NA=1017cm-3 and tox=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%.Figure 6
Idrift  versus drain voltage VD from (32) in log scale. NA=1017cm-3,tox=10 nm and W/L=10. Markers (+) corresponds to (33).Figures7(a) and 7(b) show the simulation results in linear scale. Curve (a) is compared to data reported in [4] (NA=1015cm-3,tox=50nm) and curve (b) with more recent data reported in [7] (NA=5.1017cm-3,tox=5nm). All these calculations have been performed using a constant mobility: μn=550cm2V-1s-1.Idrift  versus drain voltage VD from (32) in linear scale compared with [4] and [7]. (a) NA=1015cm-3,tox=50 nm and W/L=10, (∘) corresponds to [4]. (b) NA=5.1017cm-3 and tox=5 nm, (∘) corresponds to [7].
(a)(b)


## 6. Transconductance

The transconductance of the drift current is defined from the derivative [3](34)gm=dIdriftdVg|VD=dIdriftduS|VDduSdVg|VD.The mathematical definition of the integral operator [22] gives the condition(35)Ifg(x)=∫0xf(u)du,then⇒dg(x)dx=f(x),
and, consequently(36)dduS{∫0uSeu-ξF(u,ξ)du}=euS-ξF(uS,ξ).The transconductance is calculated at a constantVD bias. The derivative of the drain current is only on uS and (36) in (32) results in(37)dIdriftduS=qniUT2WL∫0VD/UT{μneffeuS-ξF(uS,ξ)}dξ.The transconductance is(38)gm=qniUT2WL∫0VD/UTμneffeuS-ξF(uS,ξ)duSdVgdξ,
which from (31) gives a single integral of the surface potential uS(y)(39)gm=qniUT2WL×∫uS(0)uS(L)μneffeuS-ξF(uS,ξ)(duS(y)dVg)(dξ(y)duS(y))duS(y),
with(40)duS(y)dVg=1UT[1+γ0UTe-ξdH/du+dG/du2e-ξH+G|u=uS(y)]-1.Figure8(a) shows [gm,Vg] plots in log scale when the source is grounded (VS=0,VD=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 gm=μnC0VD(W/L) is a linear function  versus VD.The transconductance from (39). (a) gm=f(Vg)(VS=0 and VD=5 V), (b) gm=f(VD)(VS=0 and Vg=5 V) (o) linear charge sheet model,. (c) normalized transconductance [G(iF),iF] plots, parameters are NA=1017cm-3,tox=10 nm and W/L=10.In the E.K.V. model [23], the authors introduce the normalized transconductance: G(iF)=gmnUT/iF with iF=IF/2nβUT2; n is the slope factor and β=μnCoxW/L. Figure 8(c) shows [G(iF),iF] plots using (39) and (32) (IF=Idrift). The function G(iF), the asymtotes G(iF)=1 and G(iF)=(iF)-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 limitsuS(0) and uS(L) given by iterative equations (20) and (24). Figure 9 shows the variations of these integration limits  versus VD, respectively, in strong (a) and weak inversion (b). The difference between uS(0) and uS(L) 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 limits of surface potentialsuS(0) and uS(L)  versus VD. (a) strong inversion, (b) weak inversion. NA=1017cm-3,tox=10 nm,W/L=10.
(a)(b)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 uLow,k which leads to Vg[uLow,k] equal to Vg with an absolute error less than 10-10.


## 8. The Diffusion Current

### 8.1. Integral of Current Diffusion

The diffusion current density of ann-MOSFET is given by(41)Jdiff(x,y)=-μnKTdn(x,y)dy.As shown in Figure1, 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(42)Idiff=W∫0xdJdiff(x,y)dx=-μnKTWddy{∫0xdn(x,y)dx}.Equation (28) gives the equivalence to potentials(43)∫0xdn(x,y)dx=niUT∫0u(s)eu-ξF(u,ξ)du,
and after integration on y, (42) becomes(44)Idiff=-qμnWLniUT2∫0Ld{∫0uS(y)eu-ξF(u,ξ)du},
and reduces to single integrals(45)Idiff=qμnWLniUT2×{∫0uS(0)eu-ξ(0)F[u,ξ(0)]du-∫0uS(L)eu-ξ(L)F[u,ξ(L)]du}.Equation (45) is equivalent to the well known expression(46)Idiff=μnWLUT[Qinv(0)-Qinv(L)].This expression is easily calculated from the exponential variations ofQinv  versus Vg in weak inversion. Equation (13) gives the exact relation Vg=f[uS(y),ξ(y)], and the inversion charge Qinv  versus uS(y) is given from (28).Figure10 shows [Qinv,Vg] plots according to (13) and (28) for different VD bias. Qinv(0) corresponds to ξ(y)=0 and Qinv(L) corresponds to the drain bias ξ(L)=VD/UT. Qinv(L) is rapidly negligible in (45). Figure 11 shows Idiff(Vg) plots in weak inversion calculated from (45) and compared with (46) from Qinv in (28).Figure 10
[Qinv,Vg] plots in log scale for different channel voltages. NA=1017cm-3 and tox=10 nm.Figure 11
[Idiff,Vg] plots in weak inversion (VD=0.3 V) (–-) from (46) and (■) from (45). NA=1017cm-3 and tox=10 nm.In weak inversion, the diffusion current is well captured by(47)Jdiff=J0eVg/NKT,
with N=1.38, which is in good agreement with the slope of [Qinv,Vg] plots (Figure 10).


### 8.2. Diffusion PlotsIdiff(VD)

A complete graphIdiff(VD) can be plotted by sampling uS(L) from (24). Figure 12 shows the contribution of the diffusion current (Idiff) compared with the drift current (Id) in strong and weak inversion (Id1>Idiff1 and Id2<Idiff2). The diffusion current is negligible in strong inversion and dominant in weak inversion.Figure 12
The driftId (–-) and diffusion Idiff (....) currents  versus VD in log scale. Id1 and Idiff1 strong inversion: Id2 and Idiff2 weak inversion. NA=1017cm-3,tox=10 nm, and W/L=10.


## 8.1. Integral of Current Diffusion

The diffusion current density of ann-MOSFET is given by(41)Jdiff(x,y)=-μnKTdn(x,y)dy.As shown in Figure1, 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(42)Idiff=W∫0xdJdiff(x,y)dx=-μnKTWddy{∫0xdn(x,y)dx}.Equation (28) gives the equivalence to potentials(43)∫0xdn(x,y)dx=niUT∫0u(s)eu-ξF(u,ξ)du,
and after integration on y, (42) becomes(44)Idiff=-qμnWLniUT2∫0Ld{∫0uS(y)eu-ξF(u,ξ)du},
and reduces to single integrals(45)Idiff=qμnWLniUT2×{∫0uS(0)eu-ξ(0)F[u,ξ(0)]du-∫0uS(L)eu-ξ(L)F[u,ξ(L)]du}.Equation (45) is equivalent to the well known expression(46)Idiff=μnWLUT[Qinv(0)-Qinv(L)].This expression is easily calculated from the exponential variations ofQinv  versus Vg in weak inversion. Equation (13) gives the exact relation Vg=f[uS(y),ξ(y)], and the inversion charge Qinv  versus uS(y) is given from (28).Figure10 shows [Qinv,Vg] plots according to (13) and (28) for different VD bias. Qinv(0) corresponds to ξ(y)=0 and Qinv(L) corresponds to the drain bias ξ(L)=VD/UT. Qinv(L) is rapidly negligible in (45). Figure 11 shows Idiff(Vg) plots in weak inversion calculated from (45) and compared with (46) from Qinv in (28).Figure 10
[Qinv,Vg] plots in log scale for different channel voltages. NA=1017cm-3 and tox=10 nm.Figure 11
[Idiff,Vg] plots in weak inversion (VD=0.3 V) (–-) from (46) and (■) from (45). NA=1017cm-3 and tox=10 nm.In weak inversion, the diffusion current is well captured by(47)Jdiff=J0eVg/NKT,
with N=1.38, which is in good agreement with the slope of [Qinv,Vg] plots (Figure 10).


## 8.2. Diffusion PlotsIdiff(VD)

A complete graphIdiff(VD) can be plotted by sampling uS(L) from (24). Figure 12 shows the contribution of the diffusion current (Idiff) compared with the drift current (Id) in strong and weak inversion (Id1>Idiff1 and Id2<Idiff2). The diffusion current is negligible in strong inversion and dominant in weak inversion.Figure 12
The driftId (–-) and diffusion Idiff (....) currents  versus VD in log scale. Id1 and Idiff1 strong inversion: Id2 and Idiff2 weak inversion. NA=1017cm-3,tox=10 nm, and W/L=10.


## 9. Effects of Source Bias

In Section5, the Pao-Sah integral is calculated with a zero bias source. If the source is biased with a voltage VSB  versus bulk silicon, the reduced drain-source potential ξ(y) along the channel varies from ξ(0)=VSB/UT to ξ(L)=VD/UT. Then, the surface potential uS,m starts at uS,mS given from (24) with mS=VSB/hUT and stops at m=M=VDB/hUT.The drift and the diffusion currents are always given by (32) and (45) by substituting uS,0 by uS,mS solution of (24) with mS=VSB/hUT.By setting(48)G(uS)=qniWLUT2∫0uSμneffeu-ξF(u,ξ)du.Equation (32) can be rewritten with m index in uS(49)Id=∫uS,0uS,MG[u(s)]dξduSduS-∫uS,0uS,msG[u(s)]dξduSduS,(50)Id=IF+IR.The integrals in (49) correspond, respectively, to the forward IF and reverse IR 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 Figure13 shows I=Idrift+Idiff=f(Vg) calculated from (32) and (45) for constant VDS and different VSB. The parameters of I=f(Vg) are VFB(flatband)=0, NA=1017cm-3,tox=10 nm, VDS=5 volt and VSB from 0 to 1 volt.Figure 13
The transfer characteristicsID=f(Vg) for constant VDB=5 volt. NA=1017cm-3,tox=10 nm and VSB=[0-1] volt.As a comparison, Figure14 shows the curves I=f(Vg) with the same parameters reported in [4]. The conditions are VFB(flatband)=0.92 for NA=1014cm-3and VFB=0.86 for NA=1015cm-3. The oxide thickness is tox=13 nm and the drain bias is VD=1 volt. Our simulations are in good agreement with these results.Figure 14
Transfer characteristicsID=f(Vg). (a) NA=1014cm-3,VFB(flat-band)=0.92 V. (b) NA=1015cm-3,VFB=0.86 V. Markers (■) and (⋄) are from [4].


## 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.Figure15 shows the user interface for current calculation. The parameters are doping NA, oxide thickness tox, and gate and drain voltages Vg and VD.Figure 15
The user interface of current and transconductance calculations.Vg*=Vg(uLow,k) calculated from uLow,k in (20).


## 12. Annex: Taylor Polynomials of Inverse Functions

If a functiony=f(x) has continuous derivatives up to (n)th order f(n)(x), then this function can be expanded in the following Taylor polynomials [25](51)f(x)=f(x0)+∑1n(x-x0)nn!f(n)(x)+R(n),
where R(n) is called the remainder after n+1 terms.This expansion converges over a certain range ofx, if limn→∞R(n)=0, and the expansion is called the Taylor Series of f(x) expanded about x0.In inverse function, ify=f(x) is a one-to-one function, then f-1 is continuous and, if f-1 has continuous derivatives up to (n)th order, f-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 functionuS(y)=f-1[ξ(y)], The first derivative variations given by (25) are shown in Figure 16 for different Vg bias. The first derivative is define in all the [0,VD] region and varies from 1 to ≈10-13 without discontinuity. The strong decrease observed for a VD value depending on Vg is due to the asymptotic value in the surface potential uS(y) when uS(y) reaches uUp, solution of E(u)-G(u)=0.Figure 16
The first derivativeduS(y)/dξ(y)  versus drain voltage VD for different Vg bias. NA=1017cm-3 and tox=10 nm.



---

*Source: 268431-2010-07-05.xml*