A perturbation approach for low frequency noise in junction field effect transistors

ii NOl~it- HOLLAND

A Perturbation Approach for Low Frequency Noise in Junction Field Effect Transistors Riccardo Sacco

Dipartimento di Matematica Politecnico di Milano Via Bonardi 9, 20133 Milan, Italy Danilo Erricolo

Politecnico di Milano P.zza L. da Vinci 32 20133 Milan, Italy and Emilio Gatti

Dipartimento di Elettronica e Informazione Politeenico di Milano P.zza L. da Vinci 32 20133 Milan, Italy

Transmitted by S. Rinaldi

ABSTRACT A mathematical analysis of low frequency noise in junction field effect transistors (JFETs) is carried out by means of a perturbation approach on the stationary semiconductor device equations. Existence and uniqueness of the solution of the perturbed system is proved to hold. Discretization of the equations by the box method is then performed and some numerical results are discussed.

1.

INTRODUCTION

In the development of low-energy X spectroscopy detectors, the channel noise coming from the input junction field effect transistor (JFET) stage

APPLIED MATHEMATICS AND COMPUTATION 74:161-190 (1996) © Elsevier Science Inc., 1996 655 Avenue of the Americas, New York, NY 10010

0096-3003/96/$15.00 SSDI 0096-3003(95)00091-U

162

R. SACCO ET AL.

plays a relevant role in the spectral noise behavior of the devices. This noise component exhibits a power spectral density that grows up at the low frequencies according to a 1 / f ~ law. In recent years a series of measurements has been performed that show a random telegraph signal (RTS) of the drain current of the J F E T due to the mechanism of trapping and detrapping of electrons by traps inside the silicon [1]. In this paper a mathematical analysis of the physical problem is carried out through a perturbation approach on the stationary semiconductor device equations. Weak formulation for the system of perturbed equations is considered together with its relative finite element approximation. Existence and uniqueness of solution of the discretized equations is then proved to hold, this implying by means of a density argument the same property for the solution of the continuous problem. The resulting finite element stiffness matrix turns out to be a nonsingular block matrix, where two blocks are symmetric and positive definite, one is diagonal strictly positive, and the last one turns out to be a -(M-matrix) (cf. [2], p. 81). In the actual discretization of the differential problem the well-known box integration method (BIM) ([2-5]) has been employed instead of the finite element method introduced before, this nevertheless leading to a stiffness matrix that shares the same properties holding for the finite element matrix described above and that is proved to be nonsingular. Some numerical experiments are then reported, concerning the direct evaluation of the amplitude of the RTS as a function of the trap location in the J F E T and of the induced charge on the gate.

2.

THE DRIFT-DIFFUSION (DD) MODEL OF A SEMICONDUCTOR DEVICE In the following, we shall deal with the very simple J F E T prototype structure represented in Figure 1. In all the numerical examples the source electrode as well as the two p+ gates are kept fixed at ground potential, while the drain electrode is suitably biased in order to have the transistor working at saturated regime (the potential difference VDs between drain and source is taken equal to 2 Volt). Because the aim of this analysis is to evaluate the RTS amplitude, we have to focus on the basic principles of this noise mechanism. The latter is a result of the capture and release of an electron by a trap located somewhere in the device and having an energy level inside the silicon forbidden gap. Once an electron is trapped, the drain current decreases because the current field distribution is permanently modified by the Coulomb field of the electron trapped; this, in turn, allows the difference between the steady drain current corresponding to an empty trap and the one relative to the

163

Perturbation Approach in Transistors

drain

source

1.5gn

FIG. 1.

Prototype transistor.

situation of an occupied trap to be taken directly as the measure of the RTS amplitude. The electric behavior of a semiconductor device at steady state is commonly modeled by the following set of partial differential equations (see, e.g., [6, 7])

{ -div(tc(V~0)) = q( p - n + N1 - NA) in ~ c

divJ n = qR div J p = - q R

R2

(1)

where Jn = - q ~ n n V O + qD,~Vn J p = - q t % pV~b -

qDpVp"

(2)

In our applications the open bounded set II is a rectangle (or the union of rectangles) and mixed boundary conditions are supplied on the boundary a l l for the unknowns ~O, n, and p. These latter denote, respectively, the electric potential and the electron and hole concentrations. In (1) and (2) N 1 is the donor concentration, NA represents the acceptor concentration, R = R ( 0 , n, p) is the net recombination rate, and J~, Jp are the electron and hole current densities; q is the electron charge, ~ is the silicon dielectric constant, /x,~ and D~ are, respectively, the electron mobility and diffusion coefficient for electrons, with an analogous meaning for/xp and Dp. Einstein's relation Dv = /xv Vth (V = n, p) is assumed to hold, being Vth the thermal voltage. Moreover, the carrier mobilities /x~ are nonlinear functions of IV01, bounded and strictly positive according to

m = m ( x , lVCL), o < m < ~ , < + ~ .

164

R. SACCO ET AL.

The set (1)-(2) is usually addressed as the drift-diffusion model for charge transport in a semiconductor device. In particular, the first equation in (1) is a classical Poisson's equation for electric potential, while (2) are two continuity equations for electron and hole current densities. From a mathematical point of view, (2) are, for a given convective field b = - V ~ , two convection-diffusion problems, with a possibly highly dominating convection in a small subset of the device domain ~ . This in turn demands a stable discretization scheme, as we shall point out later on. 3.

L I N E A R I Z A T I O N OF T H E S E M I C O N D U C T O R DEVICE E Q U A T I O N S

The physical problem that motivates our numerical investigations consists in the evaluation of the drain current I D corresponding to two distinct situations differing from each other because of the presence of an electron trapped somewhere in the device. Indeed, it should be noted that the capture of an electron by a single trap produces a current density field J ' , which does not coincide with the quantity J evaluated when no electrons are trapped. With reference to Figure 1, this gives rise to two different values of the current I v at the drain electrode, defined as

D

where S D is the electrode surface and v its related outward normal. Because we may represent a single trap as a pointwise perturbation of the doping profile, this suggests that a perturbation approach may be profitably used as a mathematical model of the physical situation. Accordingly, we impress a pointwise variation PS(x - x 0) to the doping profile in the first of (1), where x = ( x l , x2) t is the generic coordinate vector (t being the transposition operator) and x 0 is the spatial position of the perturbation P, and then perform a linearization of (1) in the neighborhood of the working point of the device. This latter is defined by the potentials impressed at each electrode, which are the Dirichlet boundary conditions for the unknown ~b in problem (1). Before considering the linearization of (1) it should be stressed that a J F E T is a unipolar device, i.e., charge transport is only due one type of carrier. We are interested in an n-channel J F E T so that set (1) may be reduced, as far as the perturbation approach is concerned, to the first two equations in the system. In order to be consistent with the former assump-

165

Perturbation Approach in Transistors

tions, we must switch the recombination-generation mechanism of (i.e., R = 0), so that the set of partial differential equations to be solved is - d i v ( e V q J ) = q( N~ - n) divJ n = 0

in lI c R 2"

(4)

where Jn is as in (2) and which must be supplemented by suitable boundary conditions for ¢ and n (usually of Dirichlet-Neumann type). Let ¢ -- ~0 + ~, n = n I + ~, N 1 = No + PS(.), where ~, ~ represent the perturbation unknowns, Pt~(.) is the impressed perturbation and q~, n 1, NO refer to the nonperturbed situation. We point out that the perturbation ti in the mobility tG (in the following denoted b y / x ) has been neglected in the linearization procedure, although the latter has been considered field-dependent in the solution of the unperturbed equations (4). The inclusion of and of the hole transport equation in the perturbation model do not present any additional difficulties and will be used in further engineering applications of the method. After some algebra the linearized equations turn out to be given by

{

-div(6V~)

=

q(P3(x

-

Xo) - ~) i n / 2 c R 2,

div( q/z( VthVfi - ~Vq~))

(5)

div(qlxnlV~) = 0

which can be written out in a more concise notation as ~(s) = f

(6)

where s = (~, h) t, = "J - d i v ( 6 V ~ )

+ qh

~,(s) - d i v ( q l x n , V ~ ) + div( q/x( VthV~ -- ~Vq~))

(7)

and

(8)

166

R. SACCO ET AL.

Homogeneous Dirichlet boundary conditions will be assumed for ~ and ~ at the electrodes (F0), while homogeneous Neumann conditions hold on the remainder of the device boundary (F 1) (see Figure 2). A perturbation analysis on equations (1) has been proposed for instance by Gnudi et al. [8]. Their approach consists of applying the linearization procedure to the discretized semiconductor device equations, while we first act at a continuous level on the first two equations in (1) and after that we consider their discretization. This leads to the same algebraic problem for Poisson's equation, which is indeed linear in {~, n, p}, while a different linear algebraic system is obtained by the discretization of current continuity equation for electrons due to its nonlinear dependence on the convective term b = -V~o. Existence and uniqueness of both weak solution s and its finite element approximation s h will be stated in the next subsection. 3.1. A Result of Existence and Uniqueness Let L2(f/) be the space of square integrable functions over the device domain f~ c R2; we introduce a Sobolev space, denoted in the following by V, which is relevant to our analysis V = H)o,D = {~b: ~b ~ L2(a),V~b ~ (L2(1~))2 and ¢ = 0on r0} Let ~ h = U L=~Ti be a weakly acute triangulation of the domain f~ ([9], p. 117), i.e, for any triangular element T ~ ~ h each angle is less than or equal to ¢r/2. Also let M = NO + N be the total number of mesh points, where NO are the nodes lying on the Dirichlet part of 3f~ and N are the nodes

r1

r1

ro

r0

D

S F

1

G2

F

1

FIG. 2. Device domain considered in the numerical simulation. The electrodes are named: S (source), D (drain), G1 (gate) and G2 (gate).

Perturbation Approach in Transistors

167

lying in the interior of ~ and on the Neumann part of the boundary. We define on the triangulation ~ h a suitable finite element space for the approximation of problem (6) Wh, 0 = {~b: ¢~ ~ C ° ( K ) , c~1T ~ P I ( T ) for all T ~ ~ h and ~blr° = 0} where Pk(T) denotes the space of polynomials having degree ~< k and defined over the element T. We shall denote in the following the couple s = (~, ~)t by (u, v) t and we introduce the change of variable v = p e ~/v,h,

(9)

which formally extends to the perturbation problem the exponential relation between free carriers and electric potential (in semiconductor device literature p is commonly referred to as the Slotboom variable [10]). Substituting (9) into equations (5), the following system in the unknowns w = (u, p)t is obtained

-div(~Vu) + a(w)

=

qpe ~/ % =

qP~(x

-

x0)

-div(qtLnlVu) + div( q ~ Vth e ~/ % V p ) = 0

(10)

supported by the boundary conditions

ulr ° . O, .

Pra.

O,.

0u Ov Ir]

O,

Op 0vlr~ - 0 .

The equivalent formulation of problem (5) in terms of the unknowns (u, p)t will turn out to be very useful in the proofs on Theorems 1 and 2. The following result holds:

THEOREM 3.1. Let be s = (~, ~t)t ~ V × V and s h = (~h, nh) t ~ Wh,0 × Wh, 0, respectively, the weak solution of ~(s) = f and its finite element approximation. Then s and s h exist and are unique. PROOF. See Appendix 1. Let us now get a step further in the discretization of problem (5). To this aim it is useful to start from problem (10), where three symmetric and definite differential operators in divergence form may be clearly identified.

168

R. SACCO ET AL.

Use of the finite element method (FEM) to find an approximation of (10) in the product space Wh, 0 × Wh, o leads to the solution of the following linear system written out in block matrix form as Kw = f

(11)

where

K =

An A21

A12)

R2N×2N

A22

in (11) w = (u, p)t ~ R2N is the vector of nodal unknowns, f = (fl, fe) t ~ R 2N where fl = ( q 6 ~ ) t ~ RN, f2 = 0 ~ R N, 6ir is the Kronecker symbol defined as r {1 6i = 0

i=r i¢ r

and N is (are) the pointer(s) to the mesh node(s) where the perturbation charge(s) is (are) located. The square block matrices A u and A21, respectively corresponding to the discretization of operators --~A(.) and -div(q/xnlV(-)) , turn out to be symmetric and positive definite; matrix A22 , relative to operator div(qVa/~ e ~/v'W(.)), is symmetric and negative definite, and matrix A~2 is diagonal with strictly positive entries. It must be stressed that the choice of variable p defined in (9) gives rise to possible overflow/underflow problems because of the presence of the term e e/v~h We can go back to the original unknown fi and then discretize problem (5) by applying relation (9) at each node of the mesh, this being equivalent to multiplying matrices A12 and A22 by the diagonal matrix D = ( e~,l v% e~21 v,,,..., e~ul y,,,), where {q0j}N=l are the nonperturbed nodal potentials. We remark that such a procedure changes the nature of matrix A22, which, although no longer still symmetric, is however a -(M-matrix) as it has been also pointed out by Brezzi et al. [11]. To minimize the programming effort, we have chosen to implement the linearized model (5) within the framework of the two-dimensional device simulator HFIELDS [5], which is a general purpose code for the numerical solution of semiconductor device equations employing the well-known box integration method (BIM). Details of the discretization procedure will be extensively described in Section 5; in particular, we shall state that the linear system arising from use of the BIM admits a unique solution, thus extending to the box method approach the result stated in Theorem 3.1 for the FEM approximation.

Perturbation Approach in Transistors 4.

169

V A L I D A T I O N OF T H E T W O - D I M E N S I O N A L SOLUTION OF T H E LINEARIZED SYSTEM

In this section we shall discuss the conditions to extend to the case of a point change excitation in the real three-dimensional space the results of a two-dimensional simulation of a device excited by a uniformly charged line. The key points of the argumentation are the following two steps: (1) The perturbation equations are linear, therefore superposition holds; and (2) the three-dimensional device in Figure 3 is supposed to be indefinite in the Z-direction. From the theoretical point of view, the effect of a point charge on the amplitude of the RTS, as computed by a three-dimensional simulator, is the same obtained with the same amount of charge distributed along an infinite line of charge and calculated by two-dimensional simulator. Indeed, let us call a . ~/ the variation of the current per unit length at the drain electrode due to the infinite line charge having linear density equal to % The latter may be thought of as composed by an infinite number of equal elements of infinitesimal length d W and charge q = ~. d W. Each of these elements (which may be considered a point charge) will contribute a variation of current [3. q = ft. ~. dW, where /8 is a proportionality factor still to be determined. The total current variation caused by the union of all these elements as W--* o¢ will then be [3.31. W, so that the current per unit length resulting from the so reconstructed infinite wire will be lim w-.~( [3~ W ) / W = [3~. Necessarily, in order not to have an absurd, a must be equal to /3. We point out that by applying the theorem above, which is rigorously valid for W ~ 0% to the case of a device of finite width W, the errors

0

L

11

FIG. 3.

Ideal three-dimensional device with exciting charge Q.

170

R. SACCO ET AL.

inherent to the approximation introduced rapidly vanish as W increases. These errors in fact arise from the component of the field in the Z-direction caused by the point charge Q = T. W, which are fields that are negligible when the charge distance from the electrode planes z = 0 or z = W is larger than K . s, where K amounts to some units. The Z-component of the fields, as is well-known, decay exponentially like exp ( - K-~), as happens for instance in the case of Coulombian fields in a waveguide of width s. It is worth noting that in most electronic semiconductor devices the width W is much larger than the depth s. Before going on with the discretization of (5), we should take a deeper insight into the linearized transport equation, namely div( qtt( VthVfi - 5Vq~)) - div( qlznlV(b ) = O.

(12)

Let us introduce the following current densities: Jl = qt~( VthVa -- fiVq~) and J2 = -(q/zn2V(~) so that the total current density of the perturbed problem may be defined by the following solenoidal vector field

•1 = J1 + "]2.

(13)

The first term of equation (12)

div(,]l)

= div( q~( VthV~ - ~Vqo))

(14)

has the structure of a classical drift-diffusion term. It may be interpreted as the divergence of a current density due to the perturbed electron distribution ~ drifted by the electric field -V~o and relative to the nonperturbed situation. Concerning the numerical approximation of equation (12), we may employ the well-known Scharfetter-Gummel [12] scheme (SG) to discretize the convection-diffusion term (14), while the second term in (12) represents a new feature of this problem and will be discussed in more detail later. 5.

D I S C R E T I Z A T I O N OF T H E LINEARIZED SYSTEM

As mentioned earlier, the device simulator HFIELDS [5] makes use of the box integration method which may be briefly explained as follows: Each node in the finite element mesh is surrounded by a subdomain, called a box, in such a way that the union of all the boxes covers Ft without overlappings (see Figure 4).

Perturbation Approach in Transistors

171

FIG. 4.

Box.

The equations are integrated over all the boxes (an consequently over 11) by a local application of the divergence theorem. In Figure 4 we show a box, which is defined by the area limited by the axis of the sides converging into the same point. We first employ the box method to discretize the linearized Poisson's equation in (5), where a piecewise linear approximation is assumed for over each element T; we obtain, denoting by II~il the box area

-sE

(oj - b----~'dij= qP~y - q~Jl~J j#=i

forall i = 1,...,

N. (15)

Sij

In (15) the sum runs over all the nodes j surrounding node i, sij and dij a r e drawn in Figure 5, and ni is the perturbed electron concentration at node i. It may be easily shown that the BIM approximation for the Laplace's operator is exactly equivalent to employing a piecewise linear finite element scheme (see, e.g., [13]). In particular, the two discretization approaches give rise to the same elemental stiffness matrix, while in general different right-hand sides are obtained unless a suitable quadrature rule and choice of the triangulation are considered. The next step is now the discretization of the linearized transport equation (12); because this latter is a convection-diffusion problem with a highly dominating convective term b = -V¢p, some care has to be taken in its numerical approximation. The first term div(J 1) = div(q/~(VthV - ~V~)) may be treated by the Scharfetter-Gummel method, which assumes the current density flowing along the edges of each element to be constant. Integration of the term

172

R. SACCO ET AL.

j

k

Fro. 5. Localtrimlgle. div(J 1) by the BIM combined with the previous assumption on the current density field gives

f div(.]l)d['~ i = f J l . V dE i ~ E J~Jdij ~i

Fi

(16)

j-~ i

where d i j is the cross-section relative to edge eij (see Figure 5) and v is the

outward oriented normal. The constant values J1'J of the current densities deriving from the SG approximation are given by

(17) where D = IZ Vth, s i j denotes the length of edge eo, B(x) = x / ( e ~ - 1) is the Bernoulli function, and A ~i~ = ~Pj - ~i. The approximation of div(Jl~ according to the SG-BIM scheme is thus given by f a d i v ( q/z( VthAfi - fiA~))

dlli

f o r all i = 1, . . . , N

(18)

Perturbation Approach in Transistors

173

where the sum runs again over the neighboring nodes j of node i. It may be proved [14] that the following accuracy estimate holds for the ScharfetterGummel finite element approximation in the one dimensional case max

In(%) - fih( xj){ < Ch2

j e [ 1 . . . . . N]

where 5h(x) denotes the approximate solution, C is a positive constant and N+ 1

= U [ i=1

being h the maximum interval in the decomposition. It must be noted that the BIM applied to the convection-diffusion equation (12), together with the Scharfetter-Gummel approximation of the current density field, may not be directly interpreted as a classical finite element scheme. Indeed, the electron concentration ~t is properly defined at each mesh node, but no shape function is available to represent it over each element T. The approximation of div(J 2) = div( - q~nlV(O) is performed as follows. We first apply the divergence theorem, obtaining

f, d i v ( - q/.enlV~) d i l l :

f, J2.vdF~ =

-frqp~ny(O.vdF~. (19)

In order to compute the contribution of J2 = -qttnlV(~ to the total flux evaluation we proceed exactly in the same way as done for the term J r Precisely, the possibly rapidly varying function n 1 is numerically interpolated along the element edge e~j by the usual Scharfetter-Gummel exponentially fitted shape functions, while the perturbation electric field - V ~ , slow-varying function is treated as in the unperturbed problem by a piecewise constant approximation over each element. Therefore, J2 turns out to be calculated with the same accuracy as "]1 and is taken into account for a proper current flux evaluation by extrapolating on all the pertinent portion of the box boundary (denoted by AB in Figure 5) its value assumed at the mid-side point Pij of the triangle edge eij as is done for J1 in the classical box method. This is equivalent in particular to assume an electron density nl(s) constant along AB and equal to its value at the mid-side point Pij of the edge e/j (see Fig. 5). If a nonvanishing electric field is present, a stable and proper evaluation of nl(Pij) comes

174

R. SACCO ET AL.

directly again from the Scharfetter-Gummel approximation, giving

nl( P#)

1 "}" e A~°~j/2Vth

+ 1

"-I- e - A ~ ' j / 2 V t h

where A qhj = ~ ; -

~z.

The expression above degenerates in the usual average of the value hi(Pj) and nl(Pi) as A ~ j - . 0, while an upwinding weighting of the two nodal electron concentrations is automatically obained if IA ~ijl >> 1. The second term in (12) is thus discretized as

f div(- q;x niV~) fi ~

qtx n 1

da i = - f F,

V~. v dFi = - E qP~, j~ i

"~-----~dijltl( ~:~ij) Sij

for all i = 1 , . . . , N so that the complete approximation of (5) is given by the following linear system, for all i -- 1 , . . . , N

j~ i

Sik

d`j B j¢- i

8ij

j¢ i

G

(20)

8ij

We may write out (20) in a more compact form as

/~4" = f

(21)

where 4, = (~, fi)t ~ R2N is the vector of nodal unknowns and the righthand side has already been defined in (11). The stiffness matrix /( may be

Perturbation conveniently

Approach in Transistors represented

175

in block form as -&A(.)

‘=

div(J,(*))

where the differential operators A(*), div(J,(.)), di&,(.)) stand for the I E RNX N is the identity matrix and corresponding stiffness matrices, A = diag(la,], ]&I,. . . , Ifi,l> E RNxN. Because the triangulation Xh employed in the numerical experiments is of weakly acute type the matrix for - div(J,(.)) turns out to be a -(M-matrix) as already pointed out in Section 3.1, while the symmetric stiffness matrix for - CA(.) is positive definite (cf., e.g., [13]). The same is true also for the stiffness matrix coming out from the approximation of div(J,(*)); this is an immediate consequence of the fact that the coefficient n, is assumed to be in the continuous problem (5) a bounded, strictly positive known function (the electron concentration). Such a property is inherited by the discretization scheme because the electron concentration computed by solving the unperturbed equations (4) by the classical box method satisfies a maximum principle and is thus _strictly positive. This ensures in particular that the differential opera+ div(J,(.)) is V-elliptic [9] and that the corresponding block matrix div(J,) is positive definite, exactly as happens for the Laplace operator considered above. From this statement it turns out that the four block matrices corresponding to the box method approximation of problem (5) share exactly the same properties of the corresponding A,, (i, j = 1,2) introduced in (11) in the case of a finite element approximation. Moreover, a similar result of existence and uniqueness for system (20) may be proved, thus extending Theorem 3.1 to the BIM discretization of problem (5). This is stated in the following

THEOREM 5.1.

The solution

PROOF. See Appendix 6.

NUMERICAL

i% = (4, fdt of (20) exists and is unique.

2.

RESULTS

In this concluding section we present some results relative to the numerical solution by the perturbation method of two physical problems of particular interest. All the computations have been carried out by suitably modifying the two-dimensional device simulator HFIELDS [5] quoted earlier, which employs the direct sparse matrices package SPARSPAK [15] to

176

R. SACCO ET AL.

solve the linear system arising from the BIM discretization of the perturbed equation (5). In the first numerical example we study a one-dimensional test problem to check the accuracy of the method by a direct comparison with the available closed form solution. The other example is more closely related to the physical problem that motivated this paper, namely, the numerical study of the RTS of the drain current in a J F E T due to the mechanism of trapping and detrapping of conducting electrons by traps located somewhere throughout the device. Results of the simulations performed on the prototype transistor in Figure 1 at both thermal equilibrium and biased situations will be discussed. Let us now describe the geometry and the main physical parameters relative to the first numerical example. The model for the device to be studied, shown in Figure 6, is a threedimensional uniformly doped resistor; the electron traps are simulated by a plane of fixed charge in the middle of the device, denoted by the shaded area in the figure. Invariance of the solution along each direction parallel to the plane of charge makes the problem purely one dimensional. According to the frame of reference as in Figure 6 the closed form solution of equations (5) is given by qPL O sinh(( x + L ) / L D ) =

sA

cosh(L/LD)

qPL D sinh(( x -

A P ~( x) =

-L
L)/LD)

cosh(L/LD)

O
sinh((x + L)/LD)

AL--D

cosh(L/LD)

P

sinh((x - L)/LD)

AL D

cosh(L/LD)

- L <.
(23)

0 < x <~ L

where A is the device cross-section and the values of the symbols introduced are listed below L = 5 • 10 -4 cm D = 10 -4 cm W = 10 -4 cm A = 10 -8 cm 2 n 1 = 1015 cm -3.

(24)

Perturbation Approach in Transistors

177

hz

0

j

2t.

FIc. 6.

x~

/

Sample device used to check the numerical routines.

The Debye length L D introduced in (22)-(23) is the characteristic device length defined as

LD=(' 1/2~ ~Vth) qnl

and physically represents the average screening length of an electrostatic field in an electron gas of density equal to n r With the data in (24) L D assumes the value 1.28.10 -5 cm, while the peak values of the solution at x = 0 /zm are, respectively, ~(0) = 1.94- 10 -4 V and ~(0) = 7.8" 1012 c m - 3.

We show in Figures 7 and 8 the discrete solutions of system (20), where a close agreement with the values above as well as the exponential decay of (22)-(23) may be clearly appreciated. It must be noted that, despite the apparent simplicity, this test problem is a real challenge because of the strong coupling between the two perturbed equations (5). This absolutely demands a full coupled solution of system (20), whereas a block iterative procedure would result in unstable oscillating perturbed potential and electron distributions, as we have observed in several simulations that have not been reported here for the sake of brevity. To assess the sensitivity of the discrete scheme to a proper mesh refinement in the neighbourhood of the perturbing line of charge, we have simulated the above considered uniformly doped resistor by placing the excitating lines of charge at four noninteracting positions throughout the

178

R. SACCO E T AL.

I

~..

4, 6. D/stance (urn)

8.

I0

FIC. 7. Perturbed potential distribution for the sample device.

T

1

A~

l: o

4.

6.

Distance (urn~

10

FIG. 8. Perturbed electron concentration for the sample device.

Perturbation Approach in Transistors

179

device. This latter has been divided in four adjacent subregions 1~ (i = 1 , . . . , 4) of length A x = 2.5 /zm each, where properly increasing values of the mesh size h~, respectively equal to h 1 = A x/128, h 2 = A x/32, h 3 = A x/8, and h 4 = A x/4, have been taken. The grid spacings, expressed in terms of the Debye length LD, are h 1 = 0.148L D, h 2 = 0.592LD, h 3 = 2.368LD, and h 4 = 4.736L D. The resulting perturbation potential ~ at thermal equilibrium is plotted in Figure 9; as a check of the accuracy of the four discretizations, we have extrapolated from the corresponding curves the values of the Debye length at distances of 2 L D away from the perturbating charge. The so computed values of L D have exhibited excellent agreement in the case of the mesh sizes h 1 and h2, while the estimated percentage errors for grids h 3 and h 4 turned out to be 10% and more than 100% respectively. Similar results have been obtained simulating the considered device with single pointwise excitating charges. In this case the analytical solutions for both charge density and potential may be shown to be proportional to the expression K 0 ( r / L D) where L D is the Debye length previously defined, r is the radial coordinate, and K 0 ( ~ ) is the modified Bessel function of order zero. We have successfully verified that also in this case the solutions for more refined grids h 1 and h 2 have been in close agreement with the expression above. To conclude this analysis, we report that numerical experiments have demonstrated that in our two-dimensional model pointwise and line charge perturbations

it1 .... IE i .... ....

I

o.

4.

FIG. 9.

6.

6.

10

P e r t u r b e d p o t e n t i a l as a f u n c t i o n of t h e m e s h size.

180

R. SACCO ET AL.

(respectively corresponding to line and plane charge perturbations in a reM three-dimensional structure) produce equal variations of the current at the terminals, provided that the total amount of charge in the two cases is the same. This expected result may be regarded as the counterpart of the theorem stated in section 4, as concerns the extension from the case of a line charge excitation ot the case of a point charge excitation. In the following we are going to discuss some numerical results relative to the simulation of the prototype transistor shown in Figure 1. We first want to investigate the effect produced by a single electron trapped somewhere throughout the device on the perturbed solution (~, ~). To this aim, let us suppose to put one perturbation charge (i.e., one trap) in the middle of the device at thermal equilibirum. The perturbed potential and electron distributions are shown in Figures 10 and 11, respectively, which exhibit a clear symmetry according to the geometry of the device, the position of the trap, and the equilibrium state of the unperturbed problem. Because the purpose of this work is to apply the perturbation analysis to evaluate the amplitude of the RTS, the following mapping procedure will be carried out. A single trap is located on a particular mesh node at x = .~, and the numerical solution of system (20) is performed. The resulting (~, g) are

< ¢o

ae

D

FIc. 10. Perturbed potential resulting from an electron trapped in the middle of the device at thermal equilibrium.

Perturbation Approach in Transistors

Z_:

z

181

h

!

s

FIG. 11. Perturbed electron concentration resulting from an electron trapped in the middle of the device at thermal equilibrium.

employed to evaluate the current density ,], which in turn allows us to determine the drain current vaiation according to (3). The procedure above is then repeated for each node of the mesh, except the ones lying on the Dirichlet part of the boundary F = F 013 F 1. The drain current variation I D as a function of the position of the trap is the main result of our numerical investigation of the starting physical problem. A plot of I D = ID(X) is drawn in Figure 12 and it actually represents what we mean by the word "mapping." The graphic shows that the most sensitive location to determine the amplitude of the RTS is the so called pinch-off zone, i.e., the region of the device where the conducting channel is thinner and the current density is higher. The obtained solution for the perturbation distributions ~ and ~ allows an easy calculation of the functional

fs~, eV(b. v dS, which gives the induced charge on the gate resulting from a trapped charge at ~ (the point in the device domain where the trap is located). The map of

182

R. S A C C O E T AI_

c~

CL

FIc. 12.

Mapping of the prototype transistor normalized to

gin~Cgs (bandwidth).

0.

FIG. 13.

Induction coefficient versus the position of the trapped charge.

183

Perturbation Approach in Transistors

the ratio of the induced charge on one of the gates with respect to the trapped charge at K is given in Figure 13. This quantity is significant in order to calculate the drain current change if the considered gate is left floating in voltage or is connected to a known load. We can observe that the induction is complete for a trap near the gate under investigation and equal to zero for a trap near the opposite gate. With the above discussed results the problem of determining the amplitude of the RTS in a prototype J F E T is completely solved and the algorithms presented in this article can be easily employed to study doping profiles, geometries, and Dirichlet-Neumann boundary conditions representing more realistic actual devices. ACKNOWLEDGMENTS We are indebted to Prof. Laura Gotusso, Department of Mathematics, Polytechnic of Milan, Italy, for several helpful and enlightening discussions during the preparation of this manuscript and also to Prof. Antonio Longoni, Department of Electronics, Polytechnic of Milan, Italy, for many useful discussions on the physical problem. We are also grateful to the referee for his suggestions and comments that have been of great help in preparing this revised version of the paper. This research has been supported by M.U.R.S.T. 40% and 60% research contracts and by I.N.F.N. APPENDIX 1 In this section we report the complete proof of Theorem 3.1. Because of the linearity of (10), to prove existence and uniqueness of the solution it is enough to examine the behavior of the kernel of .~(w). If the solution of the homogeneous problem is w = 0, the linear operator ~ ( w ) turns out to be nonsingular, which means that ~ ( w ) = f admits the unique solution w = - 1f. We first observe that w = 0 (i.e., v = p = 0) is a solution of ~ ( w ) = 0. Indeed (10) reduces to

-Au = 0

in

c~u ~lFo = O,

-~0 o~P IF1

(25)

which admits the solution u = 0. Conversely, if we set u = 0, we find v - 0. This guarantees the existence of a solution of (10).

184

R. SACCO ET AL.

Let us now consider a generic p ~ V; from the first equation in (10) with P -= 0 it follows A 't/, =

q

--V g

q e,/v,h

=

(26)

,ff

Under the assumption a(x) = q/xn 1 = constant = a, substitution of (26) into the second equation of (10) gives q

- a A u + div( q/x Vthe ~/ v'~Vp) = - a - p e */v~h + div( q/x Vthe~/ v,~Vp) = 0 . (27) Let us now define the perturbation current density

J = q/xVthe*/v'hVp = qix( VthVv- vV~p)

(28)

so that (27) becomes - d i v J + q p e ~/vth = 0

in (29)

PPr0 = O,

Op

O-'-~lrl

= 3. Vlr ~ = O.

Under the assumption that a(x) is constant, the standard weak form of (29) is: Find p ~ V such that fn-divJxdf~+

f q e 'p/v'hpX d~ = O for a l l x E

V.

(30)

Integrating by parts, taking into account the boundary conditions for .] and X, (30) becomes

faq/xVthe¢/V~hVp. Vxdl~ + f qe~/V~hpxd~ = 0 for all X ~ V (31) Letting X - P in (31) gives

fa( q/xVthe~/VthlVpl2 + qe~/Vthp2)6 d ~ = O

forall X ~ V

(32)

Perturbation Approach in Transistors

185

which implies p - 0 in ~ . This agrees with the observation t h a t the bilinear form

B( u, v) - fa(J " vv + qe~/ V'~pv) dft = ff~( qtzVthe*/Y~Vu" vv + qe*/V'~uv) da

(33)

is coercive over V, i.e., B(u, u) > / 0 and B(u, u) = 0 if and only if u = 0, for all u ~ V. Letting p = 0 in the first of (10), where again we consider P -~ 0, it results u = 0, from which we m a y conclude t h a t in the case a(x) is a constant w = (u, v) t = 0 or, equivalently, because of (9): u = v = 0 in ft. T h e general situation of a(x) nonconstant over ~ m a y be easily reconducted to an equivalent problem where a(x) is constant over each element T ~ ~h. The advantage is twofold, namely, we shall be able to prove uniqueness of the discrete solution w h and, by a density argument, t h a t operator ~ ( w ) is nonsingular. In fact, let us consider the finite element approximation of problem (10): Find (uh, Ph) E ( W h , 0 X Wh,o)such t h a t

[ ff~(~VUh" VXh + qe~/VthDhXh)

d~ =

0

fa ( qttnlVuh " VXh - qtX Vthe~/ V'hVph" VXh)

- 0

for allxh

~ Wh,o (34)

We obtain from the first equation in (34), proceeding exactly in the same way as in the case a(x) = constant,

fasVuh" Vxh dl2 = - faqe ~/ V~h, I~h __ Ah dl~ or, equivalently,

~_~f EVu h • VXh dT T "T

:

-

~., fTqe ~/Vt~phxh dT.

(35)

T

The finite element test function Xh is a linear function of x, equal to 1 at each internal node Pi and to 0 at each other node Pk surrounding Pi, so

186

R. SACCO ET AL.

that for each element T ~ ~: h it holds, denoting by 0 the angle between VXh and Vu~, C T ~- V X h ' V U

h = I V X h l IVUh]COS 0 =

constant

From (35) and the second equation in (34) in follows

fj,n~V~.Vx~da= E f q,nqVu~'Vx~dT=-F,c~ffl~,~dT. T

Let us denote a ( x ) =/xnz, and

T "~T =

fTI-tnl dT/I T]. Then we obtain

/q-(x)Vuh. vx~ da = 2

q~Tf_Vu~,vxh dT

(36)

T

and reminding (35)

E q~TfVUh T

aT

" VXh d T = - Y'~ q2aT [ e ~/V%OhX h dT. T ~ "T

Substituting this latter equation into the second one in (34), we have

fnq.n,V~" Vx~ da = faqa(x)Vu~ "Vx~ da = Y'. q-STfAVu a • VXh dT T

"11

= - - E q 2 ~ T f e'P/VthPhXh dT T ~ f~ The resulting weak form for (29) in the case of a(x) nonconstaat over then reads: Find Ph ~ Wh,0 such that

B(ph, xh) = E T

for all Xh E Wh,o-

qvthP'e*/V'hVph'Vxa + q2°tTe~/V'hp~xa '~

)

dT=O

Perturbation Approach in Transistors

187

If we choose again Xh = Ph then B( Ph, Ph) = 0 if and only if ph(x) -- 0 for all x ~ ~ . This implies also uh(x) = 0 for all x ~ 12, because substituting ph(x) into the first equation of (34), we obtain

for all Xh ~ Wh, o

A( uh, Xh) = f s V u h " VXh d12 = 0

and letting Xh = uh, we have

A(u h,uh) =fa~LVuhl2 d 1 2 = 0 .

¢* u ~ = 0

The concluding result p = 0 and, consequently, u = 0 for all x E ~ immediately follows from the density argument W~, 0 "-) V. • APPENDIX 2 In this section we give a complete proof of Theorem 5.1. We shall follow the same lines of the proof of Theorem 3.1, namely, problem (5) will be considered in the case f = 0, in order to show that the solution of the homogeneous problem is unique and coincides with the null vector. Using the same notations introduced in Section 3.1, the system of perturbed equations reads { - d i v ( ~ V u ) + qpe ~/v,~ = 0 div( qtznlVu ) + div( qtt Vth e ~/v'hVp) = O"

(37)

Integrating (37) over the device domain 12 and applying the divergence theorem gives

- ~Vu" y i d'yi + qfa p e ~/v~h d ~ i = 0 (38)

(-q,-,lVu

+ J) .,,,

= o

where the sum above extends over all the boxes {12 i}~ 1..... N and the current density J = qtzVthe~/v~hVp has been already introduced in (28). The box

188

R. S A C C O E T AL.

method approximation of the first equation in (38) gives

~ , ~ u, - uj d~j + qpie ~'/ V'~lll,[} = O. lI,

j¢ i

sii

We remark that the term within the brackets vanishes for each box of the decomposition, so that it results

E 8 j4: i

u i -- uj

d,j = - qpie ~ ' ' / " h l a J

(la ) + la ))

- qpie ~ ' / % E

=

8iJ

j• i

for all i = 1 , . . . , N where the meaning of 1 ~ ( p = 1,2, i = 1 , . . . , N) is clearly shown in Figure 14. With reference to the same figure above, we define the quantities yiPj

( p = 1,2, i = 1 , . . . ,

. . . . ~/t~lT, n 1 i j d~pi j

\\ x

N)

/ i/ \ k.

Ii

\

\

\

\

\

~

i

i

/

I

J FIG. 14. A particular of the box decomposition.

Perturbation Approach in Transistors

189

where n[ j - nl(Po). The box method approximation of the first flux term in the second equation in (38) gives

{

1 1

2 2)}

where the constants ¢1i are strictly positive. Using thepreviously introduced SG-BIM approximation of the current density flux J • v,, the second equation in (38) becomes

~, ~ ql~Vthe,P,/Vt~BI Aq~iJ} p i - pj vthf~i} 0. (39) a~ j,~ ~ Vth siJ d~j + qp~e ~'1 = We can write out this latter relation in matrix form as

(Mo + Vo)p-- M *p = 0

(40)

where the symmetric and positive definite matrix Mp corresponds to the SG-BIM approximation of the differential operator - d i v J and the matrix Do is a positive diagonal perturbation of Mp, so that the resulting sum M* turns out to be a symmetric positive definite and strictly diagonally dominant matrix. This, of course, implies that the solution of (40) is p - 0. Substituting p -= 0 into the first equation of (38) immediately gives u =- 0, which completes the proof. • REFERENCES 1 K. Kandiah, M. O. Deighton, and F. B. Whiting, A physical model for random telegraph signal currents in semiconductor devices, J. Appl. Phys. 66:2 (1989). 2 R. S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1962. 3 R.E. Bankl D. J. Rose and W. Fichtner, Numerical methods for semiconductor device simulation, IEEE Trans. Electr. Dev. ED-30 9:1031-1041 (1983). 4 E.M. Buturla, P. E. Cottrell, B. M. Grossmann, and K. A. Salzburg, Finite element analysis of semiconductor devices: The FIELDAY Program, IBM J. Res. Dev. 25:218-231 (1981).

190 5

6 7 8 9 10 11

12 13 14

15

R. SACCO ET AL. G. Baccarani, R. Guerrieri, P. Ciampolini, and M. Rudan, HFIELDS: A Highly Flexible 2D Semiconductor-Device Analysis Program, in Proc. Fourth Int. Conf. on the Numerical Analysis of Semiconductor Devices and Integrated Circuits (J. J. Miller, Ed.), Boole Press, 1985, pp. 3-12. P. Markowich, The Stationary Semiconductor Device Equations, SpringerVerlag, New York, 1986. S. Selberherr, Analysis and Simulation of Semiconductor Devices, SpringerVerlag, Vienna 1984. A. Gnudi, P. Ciampolini, R. Guerrieri, M. Rudan, and G. Baccarani. Sensitivity analysis for device design, IEEE Trans. Comp. Des. CAD=6 5:879-885 (1987). P . G . Ciarlet and J. L. Lions (Eds.), Handbook of Numerical Analysis. Finite Element Methods (Part 1), North-Holland, Amsterdam, 1991. J . W . Slotboom, Computer-aided two-dimensional analysis of bipolar transistors, IEEE Trans. Electr. Dev. ED-20:669-679 (1973). F. Brezzi, L. D. Marini, and P. Pietra, Two-dimensional exponential fitting and applications to drift-diffusion models, SIAM g. Numer. Anal. 26:1342-1355 (1989). D. L. Scharfetter and H. K. Gummel, Large-signal analysis of a silicon read diode oscillator, IEEE Trans. Electr. Dev. ED-16:64-77 (1969). R.E. Bank and D. J. Rose, Some error estimates for the box method, SIAM J. Numer. Anal., 24:777-787 (1987). R. Sacco, Convergence of a second-order accurate Petrov-Galerkin scheme for convection=diffusion problems in semiconductors, Appl. Num. Math., 11:517-528 (1993). A. George, J. Liu, and E. Ng, User guide for SPARSPAK, Waterloo Sparse Linear Equations Package, Report No. CS=78-30, University of Waterloo, 1980.