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

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.


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 V T [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][10][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.

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 div F = 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 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.
The presentation relates to a p-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 V D , V (y) ranges from 0 at the source side (y = 0) to V D at the drain side (y = L).U T = KT/q is the thermal voltage.In the following, we use the reduced potentials u b = ϕ b /U T , u(x, y) = ϕ(x, y)/U T and ξ(y) = V (y)/U T instead of ϕ b , ϕ(x, y) and V (y), respectively.u S (y) = u(x, y)| x=0 , F S (y) = F(x, y)| x=0 , n S (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 in n(x, y) with the expressions: After solving the 1-D Poisson equation dF(x, y)/dx = ρ(x, y), these expressions give the 1-D electric field by and its derivative dF x, y dξ y As e ub e u(x,y) in the channel region, this expression is equivalent to dF x, y dξ 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 → ∞) = N D , (4) must be corrected by   20) and (b) equation (22).N A = 10 17 cm −3 and t ox = 10 nm.Vg = V T when (a) crosses the line With this exact formulation of the donor density, the electric field gradient from the 1-D Poisson equation ( 5) now reads which gives the electric field where we have defined In this form, the derivative dF x, y dξ 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 F S (y) = F(x, y)| x=0 is given by substituting u and u(x, y) by u S (y) in ( 10)- (11).

The Surface Potential Dependence to Gate Bias
From electrostatic considerations [14], the gate voltage is expressed by The exponential expression in (13) allows to separate u(s) and ξ(y) with the derivative   1) and (14).N A = 10 17 cm −3 and t ox = 10 nm.
Drain voltage V D (V) in which we introduce the dimensionless quantity Several approximations to (13) have been reported.For instance, the equation is often used in practice in a large variety of surface potentialbased 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:

The Surface Potential Dependence to Drain Bias
The explicit relation ξ(y) = f [u S (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 [u S (y)] and to the inverse function u S (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 u S (y) compatible with the denominator E[u S (y)] − G[u S (y)].

Boundary Limits of Surface Potential at Constant Gate
Voltage.Equation ( 14) has only physical meaning when E[u S (y)]−G[u S (y)] > 0. At a constant Vg value, the reduced surface potential u S (y) ranges between u Low (lower value), solution of (13) with ξ(y) = 0  and between u Up (upper value), solution of (13) with ξ(y) → ∞ leading to e −ξ(y) → 0 u Low and u U are defined by implicit relations to the gate voltage Vg.Following a methodology previously described in [17], u Low and u Up 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/δ, u Low and u Up at iteration k are calculated values at previous iteration according to In equations ( 20) and ( 22), the initial value k = 0 corresponds to flat band (Vg = 0), and the starting conditions are u Low,0 = u Up,0 = u b .
In (20), u Low,k reaches the threshod voltage V T when: k = k T = int(V T /δ); with the function int = (number down to the nearest integer).As a result, V T = 1.2819 if N A = 10 17 cm −3 and t ox = 10 nm and k T = 12819.u Low,kT = 15.712= −u b when the step is δ = 10 −4 and the corresponding band bending is 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 u Low and u Up 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 u Low and u Up are used as integration limits in Pao-Sah integral (see Section 7).14) is an explicit relationship between u S (y) and ξ(y).As this has been done in previous section with the calculation of the lower (u Low ) and upper (u Up ) limits of the surface potential, the inversion of ( 14) can be obtained by a first order Taylor expansion of the inverse function u S (y   setting ξ(y, m) = m • h in which h is the sample step size and m an integer:

The Inversion of ξ(y)
In (24), initial value for u S is u S,0 = u Low,k ; u Low,k from (20) with k = Vg/U T , and iteration stops at m = M = V (y)/hU T .
Note that this representation is u S (y) versus V (y) and not u S (y) = g(y) which would need the knowledge of the variation of V (y) versus y.Such 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.
Figure 4 shows a complete set of curves u S (y) versus V (y) from Vg = 1.2 up to 5 V. u S (y) = u S,m is calculated from (24) and ξ(y) from u S,m in (14).This figure illustrates the pinch-off voltage in strong inversion which appears when u S (y) − ξ(y) becomes to decrease.
In terms of surface electron density n S (y) = n i e uS(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 dn S /dy gradient.
In strong inversion, the transition between n S (y) = cte and the n(y)αe −ξ(y) regimes happens for a voltage V t (Vg).The inset plots in Figure 5 shows that V t (Vg) = 0 when Vg is equal to the threshold voltage 28 V defined in the charge sheet model [18].

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: Or, in terms of potential J u x, y , ξ y = qμ n n i e [u(x,y)−ξ(y)] F u x, y , ξ y . (27)

Active and Passive Electronic Components
With the introduction of the inversion charge density where we recall F = −U T (du/dx) and x d the inversion length in x direction; the expression for the drain-source drift current follows in which a total channel width W is assumed.The Pao-Sah double integral then reads Since ξ(y) can be expressed in terms of u S (y) or alternatively u S (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.(14).In (30), dξ(y) is replaced by

Solution from Surface Potential u(s). ξ(y) is a function of u S (y) from
in which dξ(y)/du S (y) from (15), is function of u S (y).The Pao-Sah double integral can be expressed in terms of two single iterated integrals [19], and the drift current is given by μ neff e u−ξ(y) F u, ξ y du dξ y du S y du S y . (32) The integral in braces is integrated first and gives a function of u S (y) and the final integral is integrated with respect to u S (y) from u S (0) to u S (L).
In (32): (i) u S (y) is the surface potential along the channel.
(iii) dξ(y)/du S (y) is given by ( 15).(iv) u S (0) is the surface potential in y = 0.At a constant Vg bias, it corresponds to u Low,k solution of (20) with Vg, it corresponds to u S,m solution of ( 24) with m = V D /hU T .(vi) μ neff is the effective mobility which has extensively been studied previously [20,21].The correction over a constant mobility (μ n = 550 cm 2 • V −1 s −1 ) used in this paper, principally depends on F[u, ξ(y)] and is easy to implement in the integral.

Solution from Channel Potential ξ(y)
= V (y)/U T .In the drift current (30), it is possible to define the surface potential u S 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 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.

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 N A = 10 17 cm −3 and t 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] (N A = 10 15 cm −3 , t ox = 50 nm) and curve (b) with more recent data reported in [7] (N A = 5.10 17 cm −3 , t ox = 5 nm).All these calculations have been performed using a constant mobility: μ n = 550 cm 2 V −1 s −1 .

Transconductance
The transconductance of the drift current is defined from the derivative [3] The mathematical definition of the integral operator [22] gives the condition and, consequently The transconductance is calculated at a constant V D bias.The derivative of the drain current is only on u S and (36) in (32) results in The transconductance is which from (31) gives a single integral of the surface potential u S (y) with Figure 8(a) shows [g m , Vg] plots in log scale when the source is grounded (V S = 0, V D = 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 In the E.K.V. model [23], the authors introduce the normalized transconductance: respectively, in weak and strong inversion are in good agreement with this model.

The Accuracy of Pao-Sah Integral
We have integrated the Pao-Sah equation between the limits u S (0) and u S (L) given by iterative equations (20) and (24).Figure 9 shows the variations of these integration limits versus V D , respectively, in strong (a) and weak inversion (b).The difference between u S (0) and u S (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 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 u Low,k which leads to Vg[u Low,k ] equal to Vg with an absolute error less than 10 −10 .

Integral of Current Diffusion. The diffusion current density of an n-MOSFET is given by
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 Equation ( 28) gives the equivalence to potentials and after integration on y, (42) becomes and reduces to single integrals du .
(45) Equation ( 45) is equivalent to the well known expression This expression is easily calculated from the exponential variations of Q inv versus Vg in weak inversion.Equation (13) gives the exact relation Vg = f [u S (y), ξ(y)], and the inversion charge Q inv versus u S (y) is given from (28).
In weak inversion, the diffusion current is well captured by with N = 1.38, which is in good agreement with the slope of [Q inv , Vg] plots (Figure 10).

Diffusion Plots I diff (V D
).A complete graph I diff (V D ) can be plotted by sampling u S (L) from (24). Figure 12 shows the contribution of the diffusion current (I diff ) compared with the drift current (I d ) in strong and weak inversion (I d1 > I diff1 and I d2 < I diff2 ).The diffusion current is negligible in strong inversion and dominant in weak inversion.

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 V SB versus bulk silicon, the reduced drain-source potential ξ(y) along the channel varies from ξ(0 Then, the surface potential u S,m starts at u S,mS given from (24) with m S = V SB /hUT and stops at m = M = V DB /hUT.The drift and the diffusion currents are always given by (32) and (45) by substituting u S,0 by u S,mS solution of (24) with m S = V SB /hU T .
By setting The integrals in (49) correspond, respectively, to the forward I F and reverse I R currents introduced in the E.K.V model [23].The same expressions can also be defined in the diffusion current in (45).

Transfer Characteristics
The transfer characteristics represented in Figure 13 shows I = I drift + I diff = f (Vg) calculated from (32) and (45) for constant V DS and different V SB .The parameters of I = f (Vg) are V FB (flatband) = 0, N A = 10 17 cm −3 , t ox = 10 nm, V DS = 5 volt and V SB from 0 to 1 volt.
As a comparison, Figure 14 shows the curves I = f (Vg) with the same parameters reported in [4].The conditions are V FB (flatband) = 0.92 for N A = 10 14 cm −3 and V FB = 0.86 for N A = 10 15 cm −3 .The oxide thickness is t ox = 13 nm and the drain bias is V D = 1 volt.Our simulations are in good agreement with these results.

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 N A , oxide thickness t ox , and gate and drain voltages Vg and V D .

Annex: Taylor Polynomials of Inverse Functions
If a function y = 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] where R(n) is called the remainder after n + 1 terms.This expansion converges over a certain range of x, if lim n → ∞ R(n) = 0, and the expansion is called the Taylor Series of f (x) expanded about x 0 .
In inverse function, if y = 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 function u S (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, V D ] region and varies from 1 to ≈10 −13 without discontinuity.The strong decrease observed for a V D value depending on Vg is due to the asymptotic value in the surface potential u S (y) when u S (y) reaches u Up , solution of E(u) − G(u) = 0.

Nomenclature
Vg: Voltage gate with bulk silicon grounded ε S and ε ox : Silicon and silicon oxide permittivity Numerical applications use SI units, excepted dopant concentration in cm −3 , ε S and ε ox in farads.cm−1 .

Figure 1 :
Figure 1: The channel scheme in the Gradual Channel Approximation.

Figure 2 :
Figure 2: Upper u Up and lower u Low limits of the reduced surface potential versus gate voltage Vg.(a) equation (20) and (b) equation(22).N A = 10 17 cm −3 and t ox = 10 nm.Vg = V T when (a) crosses the line u Low = −u b .

Figure 15 :
Figure 15: The user interface of current and transconductance calculations.Vg * = Vg(u Low,k ) calculated from u Low,k in (20).

Figure 16 :
Figure 16: The first derivative du S (y)/dξ(y) versus drain voltage V D for different Vg bias.N A = 10 17 cm −3 and t ox = 10 nm.

t
ox : Oxide thickness L: Channel length W: Channel width n i : Intrinsic carrier concentration in cm −3 N A and N D : Dopant concentrations in cm −3 U T = KT/q: Thermal voltage ϕ(b) = ϕ b = −U T ln(N A /n i ): Bulk potential of pdoped silicon u(x) = ϕ(x)/U T : Reduced potential V T = −2ϕ b + γ |2ϕ b |/U T : Threshold voltage in strong inversion γ 0 = (1/C ox ) 2KTε S n i : Intrinsic body factor with C ox = t ox /ε ox γ = (1/C ox ) 2KTε S N A : p-type semiconductor body factor.