h-..
'l .mu}:

, .
.rvd.’

L1.»

m
a.“ _

u .3)
DI‘ "
nfiwnuwm.
.J...h¢¢ .
.d . .Lvh “v; 1
“Sr“
1'

.. Av..-
...:;~v. .Yt’u
v .1 II.;\

 

I... if}... I.
‘I: g ._ 4:!!3. vv
.3 $2311.30}.

. .
'.‘. Eli?! .
3.5:}. 515%:

. , 2

‘\
\. :53»):

L

25
tech.»

23‘
)
{

-E
1"
-‘

.5‘“?

.
U

{I'MUHQWL .3. r
“firm?

Pie... .
. €5.25... ‘1!!! .5} .
yllv,!-I - ."on-unit..i‘l ‘3-"
av ' ’ furl-ID {.(‘t I
ovlplr

\vat.. 1!;

 

f . .
. .11? o p .9.
1 lkvhui...»z.u v. . (r.

 

 

 

 

 

 

 

....u..
..L ufiJ» r . n . .. ,
5.331....7 C.‘ , . . . . . .
. . _ . . . . . . . 931.2252”?
. . ‘ .H. .. u.. . .. ...‘ sCIII... Y‘i...u
. .
v I . . ‘n . ' . ‘ . - .. . Ln birth
«‘L‘.I.og;w.nl.k I L. ..I|I.I.'O it. . 1'1! .1..1.Al...d..l“1 .4 I 4 .‘d‘¢.)ifiyt rule .141 25!.- I.D.Ilu. 3940. 1 l. o: .‘ .I. .4, 7, 4-113 2, x. . .nulA. ‘ . II
I | | | Ill ill. III'II .
‘ t ufl

‘ at Smurawflgfihfinkw? . .

HIGAN STATE UNIVERSITY LIBRARIES

mflimmwmuu:I:;:I1:n\i.minl

5 ‘\ 3 1293 00594 8_42

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ruamv
Michigan State
i University

 

 

This is to certify that the

dissertation entitled

AC SIMULATION OF FIELD EFFECT TRANSISTORS
WITH A HYDRODYNAMIC TRANSPORT MODEL

presented by

Yao-Tsung Tsai

has been accepted towards fulfillment
of the requirements for

Ph . D. degree in Electrical Engineering .

 

 

‘ , ' I
Ma' professor

Date 5/3/?0

MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771

 

PLACE IN RETURN .BOX to remove We checkout from your record.
TO AVOID FINES return on or before due due.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative AOtIONEMRpOflufl lnetItution

AC SIMULATION OF FIELD EFFECT TRANSISTORS
WITH A HYDRODYNAMIC TRANSPORT MODEL

by

Yao-Tsung Tsai

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of '

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering
1990

ABSTRACT

AC SIMULATION OF FIELD EFFECT TRANSISTORS
WITH A HYDRODYNAMIC TRANSPORT MODEL

by

Yao-Tsung Tsai

Improvements in compound semiconductor materials and processes are producing
small- geometry, high-frequency field-effect transistors (FETs) for digital, microwave
and millimeter-wave applications. The design of the FETs requires accurate and
efficient computer-based modeling tools which include the appropriate physical
phenomenon. A computer-based simulation tool using the semiconductor hydro-
dynamic transport equations has been developed, evaluated and applied for studying
the DC and AC behavior of III-V FETs including MESFETS and MODFETs.

The hydrodynamic transport model which solves the continuity equation, the con-
servation of energy equation and the Poisson equation is the basis of this simulator.
The model does the DC solution by numerically solving these equations using the
box-integration method for discretization and the Newton method for the discrete equa-
tion solution. The DC solution was verified by comparison to a Monte Carlo particle
solution of the FET structure. The AC solution was accomplished by applying the
sinusoidal small-signal analysis (S 3A) technique to the hydrodynamic transport equa-
tions. This is the first time that the hydrodynamic transport equations have been
solved for the AC solution by this S 3A technique. The AC solution was compared to
the AC solution found using the Monte Carlo particle simulator with a Fourier decom-
position solution. FET structures of various geometries were studied with respect to

their AC performance. The AC parameters simulated include y-parameters, current

Yao-Tsung Tsai

gain, unilateral power gain, f T and f m. In particular, the performance of FETs in the

millimeter-wave frequency range was simulated and studied.

ACKNOWLEDGMENTS

I would like to thank my adviser, Professor Timothy A. Grotjohn, for his constant
guidance, encouragement, valuable suggestion and inspiring discussions throughout the

course of this study.

Thanks are also given to guidance committee members, Professors D. K.
Reinhard, G. M. Wierzba and C. Y. Wang for taking time to serve on the examining
committee.

Finally, I am also grateful to my wife Kitty for her support during this research.

iv

Table of Contents

 

 

 

 

LIST OF TABLES ..................................................................................................... viii
LIST OF FIGURES ................................................................................................... ix
Chapter 1. INTRODUCTION .................................................................................... 1
l. 1 Semiconductor Device Modeling .................................................................... 1
1.2 Statement of Purpose - - _- - ....... -- ..... 3
1.3 Thesis Preview -- - - - .............. - - - - .................... 3
Chapter 2. HYDRODYNAMIC SEMICONDUCTOR EQUATIONS ..................... 7
2.1 Model History .............................................................. - .....
2.2 Model Development .........................................................................................
2.3 Relationship Between HTM and DDM .......................................................... 15
2.4 GaAs MESFET Simulation ............................................................................. 19
2.5 Energy-Dependent Parameter Calculation. .................................................... 23

Chapter 3. NUMERICAL SOLUTION OF THE SEMICONDUCTOR

EQUATIONS ........................................................................................... 28

3.1 Finite-Difference Method ................................................................................ 28
3.2 Finite-Element Method .................................................................................... 30
3.3 Box-Integration Method ................................................................................... 32
3.4 Scaling ........... ‘ ................................................................................................... 35
3.5 Modified Scharfetter-Gummel Technique for J and S. ................................ 36
3.6 Matrix Equation Solution ................................................................................ 42
3.6.1 Review of Matrix Equation Method. ....................................................... 42
3.6.2 Matrix Coefficient Allocation and Calculation. ....................................... 44
3.7 Boundary Conditions ....................................................................................... 49

vi

Chapter 4. DC SIMULATION OF MESFETS USING THE
SEMICONDUCTOR HYDRODYNAMIC TRANSPORT

EQUATIONS ........................................................................................... 51
4.1 Current-Voltage Characteristics Compared with Monte Carlo. .................... 51
4.2 Source and Drain Resistance Studies of Short Channel MESFET’s .

Using Two-Dimensional Device Simulators .................................................. 65
4.2.1 Source and Drain Resistance Models ........................................................ 66
4.2.2 Bias Dependence of Source and Drain Resistance ................................... 70
4.2.3 Source and Drain Resistance in Submicron GaAs MESFETS ................. 74

Chapter 5. AC SIMULATION OF MESFETS USING THE
SEMICONDUCTOR HYDRODYNAMIC TRANSPORT

EQUATIONS ........................................................................................... 86
5.1 AC Simulation Methods .................................................................................. 86
5.2 Sinusoidal Steady-State Analysis ( 53A ) using
Hydrodynamic Transport Model ..................................................................... 89
5.2.1 Model Development .................................................................................. 89
5.2.2 Boundary Condition .................................................................................. 94
5.2.3 Small Signal Current Calculation ............................................................. 94
5.2.4 Y-Parameters Calculation ......................................................................... 96
5.3 H'I'M Y-Parameters Compared with Monte Carlo and
Drift Diffusion Models .................................................................................. 97
5.4 AC Performance of GaAs MESFETS ............................................................ 102
5.4.1 Drain Bias Dependence ............................................................................ 105
5.4.2 Gate Bias Dependence .............................................................................. 106
5.4.3 Gate Length Dependence ........................................................................ 106
5.4.4 Gate-Source Spacing Dependence ............................................................ 106
5.4.5 Gate-Drain Spacing Dependence .............................................................. 106
5.4.6 Epilayer Thickness Dependence ............................................................... 106
5.4.7 Substrate Dependence ............................................................................... 107
Chapter 6. HYDRODYNAMIC TRANSPORT MODEL FOR MODFET ............. 116

6.1 Hydrodynamic Transport Equations for the MODFET .................................. 116

vii

6.2 Simulation Results and Discussion ................................................................. 122

Chapter 7. CONCLUSIONS AND RECOMMENDATIONS .................................. 131

Appendix A. FINlTE-ELEMENT DISCRETIZATION FOR
DRIFT-DIFFUSION MODEL ............................................................ 133

Appendix B. DERIVATION OF THE CURRENT DENSITY AND
ENERGY FLUX EXPRESSIONS USING THE
MODFIED SCHARFETI'ER-GUMMEL TECHNIQUE ................... 137

Appendix c. REGRID ALGORITHM ...................................................................... 142

LIST OF REFERENCES ........................................................................................... 145

LIST OF TABLES
Table Page

3.1 Scaling factor ................................................................................. 37

5.1 f m dependence on device parameters for ND = 2><1017cm‘3 115

viii

LIST OF FIGURES

Figure Page
1.1 Different device modeling levels .................................................. 2

1.2 (a) Average drift velocity of electrons in GaAs
versus time for different fields. (b) Drift velocity
relaxation time versus the electric field ........................................ 4

2.1 The displaced Maxwellian distribution can be uniquely

described by three quantities: n, V, and T ................................... 8
2.2 Computer flow chart for (a) drift diffusion model,

(b) hydrodynamic transport model ................................................ 18
2.3 GaAs band su'ucttu'e ...................................................................... 24

2.4 Results of Monte Carlo simulation for GaAs with
N0 =2x1017cm‘3 and T0=300°K :
(a) average electron velocity v , (b) total electron
energy fi, and (c) electron potential energy U p (fi)
versus electric field E .................................................................... 26

2.5 Results of Monte Carlo simulation for GaAs with
N0 =2x1017cm’3 and To=300°K :
(a) electron energy relaxation time t§(§), (b) electron
mobility p.(§), and (c) ratio of thermal energy to total

energy, r(§), versus electron total energy é ................................. 27
3.1 Finite-difference mesh point notation ........................................... 29
3.2 Non-uniform finite-difference mesh .............................................. 29
3.3 Finite-element mesh ....................................................................... 31

3.4 Box-integration method based on (a) point basis, (b) rectangular

element basis, and (c) triangular element basis ........................... 33
3.5 Notation for both I and S discretization .................................... 39
3.6 Grid number assignment for a rectangular mesh ......................... 45

Figure Page
4.1 Simplified two-dimensional MESFET geometry ......................... 52
4.2 The current-voltage characteristics for both HTM and

Monte Carlo .................................................................................. 54

4.3 Grid structure used in the two-dimensional Monte Carlo
simulation ........................ 55

4.4 Potential distribution (V) over the entire device for

VGS = 0.0 V and VDS = 1.5 V .................................................... 57
4.5 Electron density distribution (m ’3) over the entire

dCViCC for VGS = 0.0V and VDS = 1.5 V ................................... 58
4.6 Electron energy distribution (eV) over the entire

dCViCC for VGS = 0.0V 311d VDS = 1.5V .................................. 59
4.7 Potential distribution at y=0.l pm for

VGS = 0.0V and VDS = 1.5V ..................................................... 61
4.8 Electron density distribution at y=0.l um for

VGS = 0.0V and VDS = 1.5V ..................................................... 62
4.9 Electron energy distribution at y=0.1 pm for

VGS = 0.0V and VDS = 1.5V ..................................................... 63
4.10 Current voltage characteristics for both HTM and DDM .......... 64

4.11 MESFET structure showing the lateral gate depletion
regions ALGS and ALGD .............................................................. 68

4.12 Source (solid-line) and drain (dashed line) resistance
versus gate voltage for a silicon MESFET using the power
dissipation definition .................................................................... 72

4.13 Source resistance versus drain voltage with a gate voltage
of 0.3 volts ................................................................................... 73

4.14 Change in the lateral depletion width ALGD versus
drain voltage ................................................................................. 75

4.15 Drain resistance between cross section A (shown in Figure
4.11) and the drain contact versus drain voltage ........................ 76

xi

Figure Page
4.16 Drain resistance versus drain voltage for a silicon MESFET 77
4.17 Power dissipation (dashed line) and electron heating

(solid line) in the drain resistance region Rm versus

drain voltage for a LG = 0.4 pm GaAs MESFET ..................... 81

4.18 Electron temperature versus position between the source
and the drain as calculated using the energy transport
simulator ....................................................................................... 82

4.19 Equivalent drain resistance versus drain voltage for a
LG = 0.4 pm GaAs MESFET ..................................................... 83

4.20 Total electron energy above ti; I‘ valley minimum
(solid line) and the electron kinetic energy (dashed line)
versus position between the source and drain at a depth
of 0.08 pm .................................................................................. 85

5.1 Fourier decomposition method ..................................................... 87

5.2 Charge partitioning method for (a) a two terminal device,

and (b) a three terminal device .................................................... 90
5.3 Notation for small- signal current density calculation ................. 95
5.4 The three-terminal y-parameter equivalent circuit ...................... 95

5.5 A comparison of y-parameters versus frequency for
both the HTM method and the Monte Carlo method ................. 98

5.6 A comparison of y-parameters versus frequency for
both the HTM method and the DDM method ............................ 99

5.7 The three gains A1, Av and 0;, versus
frequency from the HTM ............................................................. 100

5.8 The three gains Al, AV and GU versus
frequency from the DDM ............................................................ 101

5.9 (a) Circuit model of a MESFET. (b) The extrapolation of
the low-frequency gain overestimates the unity power gain
frequency, f m ............................................................................ 104

xii

Figure Page
5.10 (a) GU versus frequency for various VDs-

(b) f m versus VDS ..................................................................... 108
5.11 (a) GU versus frequency for various Vas-

(b) f m versus VGS ..................................................................... 109
5.12 (a) 0,, versus frequency for various L6.

(b) f m versus LG ...................................................................... 110
5.13 (a) GU versus frequency for various L65 .

(b) fm versus L63 ..................................................................... 111
5.14 (a) 6,, versus frequency for various LCD.

(b) flmx versus LCD .................................................................... 112
5.15 (a) GU versus frequency for various epilayer thickness a .

(b) f“m versus a ......................................................................... 113
5.16 (a) 0,, versus frequency with or without substrate.

(b) ND versus depth ..................................................................... 114
6.1 Basic device topology of the MODFET ...................................... 117
6.2 Energy band diagram for a MODFET ......................................... 117
6.3 Equilibrium ohmic contact boundary condition model ............... 121

6.4 The current-voltage characteristics of the MODFET for .
VGs = 0.0V and VGS =0.2V ......................................................... 123

6.5 The conduction band distribution (eV) over the entire
device for VGs = 0.0V and VDs = 1.0V ................................... 124

6.6 Electron density distribution ( m‘3 ) over the entire
device for VGs = 0.0V and VDS = 1.0V ................................... 125

6.7 Electron energy distribution (eV) over the GaAs layer for
VGS = 0.0V and VDS = 1.0V ..................................................... 127

6.8 Longitudinal electron current density distribution (A/m 2)
over the entire device for VGs = 0.0V and Wm = 1.0V .......... 128

6.9 The y-parameters versus frequency for the MODFET
using HTM ................................................................................... 129

xiii

Figure Page
6.10 GU and A, versus frequency for the MODFET using HTM 130
Appendix

Figure Page

C.l Regrid algorithm for rectangular mesh ....................................... 143

CHAPTER 1
INTRODUCTION

1.1 Semiconductor Device Modeling

Modeling for semiconductor devices means to produce a representation or simula—
tion of a device, or to make a description or analogy which helps to visualize the dev-
ice characteristics that can not be directly observed. Modeling is often done by solving
the appropriate mathematical equations that describe the device operation. The sem-
iconductor equations consist of a set of partial differential equations which must be
solved subject to a pre-defined set of boundary condition over a specified domain.
There are two approaches to describe the device behavior: (1) closed-form analytical
solution, and (2) numerical solution. In many circumstances, it is possible to simplify
the device model to be a closed-form analytical expression. The solution can be
directly computed using this analytical expression with minimal computer time, hence,
the analytical solution is suited for circuit-level simulation. However, they are
severely limited in their range of application and accuracy because of the multi-
dimensional and non-ideal nature of most modern devices. The numerical approach
requires considerably more computer time than the analytical method, but usually pro-
duces more accurate results and provides greater flexibility. Figure 1.1 shows the

different device modeling levels.

In order to obtain higher speed and higher integration performance, the size of
semiconductor devices has been drastically reduced to submicrometer dimensions due
to progress in fabrication technology. Numerical simulation techniques are required
for submicrometer channel length devices in order to design and improve the devices
for various applications. The numerical models used to model semiconductor devices

may be classified into two types depending on the solution techniques. The two

 

 

 

Closed-f
method

 

 

 

 

Modeling

 

 

 

 

 

Numerical
method

 

 

 

 

 

 

 

 

 

 

__‘

Monte Carlo
particle model

Figure 1.1. Different device modeling levels.

 

 

 

3

technique are a continuum approach treating the macroscopic quantities of carrier tran-
sport directly and the particle approach treating the macroscopic quantities as averages
of many microscopic events. The continuum approach has a variety of methods
including the drift diffusion model which has received the most use in the past. The
drift-diffusion model (DDM) or classical model is based on solving the continuity
equation and Poisson equation. This makes two central assumptions: (1) there is a
steady-state thermal equilibrium between mobile carriers and the crystal lattice, and (2)
there is a stationary relationship between the electric field and the carrier velocity
(local field dependent mobility). In FET’s with submicrometer channels, the electrons
are not transported under the carrier-lattice equilibrium conditions, and nonstationary
effects such as velocity overshoot become significant. The non-stationary effects can
be understood with the aid of Figure 1.2 [1]. When a constant electric field is applied
in a homogeneous uniformly-doped material, it takes some time for elecu'ons to reach
steady state. Figure 1.2(a) shows that velocity overshoot occurs before steady state has
been reached. The velocity overshoot can improve the device performance for submi-
crometer devices. DDM makes the stationary assumption which is equivalent to using
the velocity values at the times greater than 5ps in Figure 1.2(a), and neglects the velo-
city overshoot effect. Figure 1.2(b) calculates the delay time for electrons to reach
within 95% of their steady state value. Therefore, the drift-diffusion model does not
accumme predict importrnent phenomena related to nonstationary effects which occur

in today’s small- geometry semiconductor devices.

For accurate modeling of nonstationary effects, two major approaches have been
developed. The first is the Monte Carlo particle simulation [1-4] which gives accurate
solutions based on detailed transport and band structure parameters but at a large
expense in computation time. This large computation time is unattractive for device
optimization and the extraction of circuit simulation model parameters. The second is

the hydrodynamic transport simulation [5-11] which solves the conservation equations

 

 

 

 

 

 

 

&
I

 

Relaxation time, ps
AIL

l

l l
0.5 1.0 1.5 2.0
Electric field, MVm’l

 

(b)

Figure 1.2. (a) Average drift velocity of electrons in GaAs versus time for different
fields. (b) Drift velocity relaxation time versus the electric field. This relaxation has
been defined in the text.[l]

5

for particles, momentum, and energy. This method greatly reduces the amount of CPU
time compared with the Monte Carlo simulation method.

1.2 Statement of Purpose

Due to nonstationary effects, the DC, AC, and transient performance of field
effect transistors must be carefully evaluated. Tire primary goal of this study is to
develop a two—dimensional device simulator using the hydrodynamic transport model
(HTM) for III-V FETs. The simulator should include the following features:

a) solution of the three conservation equations (particle conservation, momentum con-

servation, and energy conservation) coupled with Poisson’s equation,
b) DC and AC simulation capabilities,
c) transport parameters calculated from Monte Carlo simulations,
d) regrid capabilities based on potential, electron concentration and electron energy,
and
e) heterojunction capabilities for the simulation of heterostructure FETs.

The simulator will be used to model the high-frequency (1-100’s GHz ) behavior
of III-V FETs ( MESFETS and MODFETs ), and to compared to the results obtained
from the Monte Carlo method or to the results available in the literature. The results
of this study will not only provide information on device performance, but will also
contribute to the basic understanding of device modeling using hydrodynamic transport
model.

1.3 Thesis Preview
The thesis is divided into seven chapters. Chapter 2 develops the hydrodynamic
transport model (HTM), reviews previous works and presents the model for this

research. Chapter 3 describes the numerical technique for this model. In particular,

the discretizations for current density and energy flux are described in detail. The DC

6

simulation solution of the MESFET is discussed and verified using a Monte Carlo par-
ticle simulation in Chapter 4. The sinusoidal small-signal analysis (S 3A) technique
applied to the HTM is described and verified in Chapter 5. The AC performance
dependence on device parameters is also studied. The HTM model is extented to the
DC and AC simdations of the MODFET in Chapter 6. Conclusions are given in the

last chapter.

CHAPTER 2
HYDRODYNAMIC SEMICONDUCTOR EQUATIONS

This chapter starts with a review of the literature on the hydrodynamic transport
model followed by the development of the model. The relationship between the hydro-
dynamic transport model (HTM) and the drift-diffusion model (DDM) will be dis-
cussed. Also, the HTM used in previous work and in this study will be included in

this chapter.
2.1 Model History.

One of the earliest solutions of the Boltzmann traIrsport equation (BTE) for the
hot-electron problem, which was to influence subsequent developments in the field,
was given by Stratton[12]. Stratton expanded the nonequilibrium distribution function
in spherical harmonics, assuming the first two terms to be important. Applications to
device simulation using Stratton-based transport model may be found to date in Cu

and Tang [13], Cook[14], and Widiger et al. [15].

A further development that was more general in approach was given by Blo-
tekjaer[16]. Blotekjaer solved the transport equations for electrons in a two-valley sem-
iconductor. Instead of expanding the distribution function into harmonics as Stratton,
Blotekjaer derived the model by taking the first three moments of the BTE to give the
macroscopic quantities electron concentration n,- , drift velocity V} , and electron tem-
perature T ,- ,where i indicates the valley (see Figure 2.1). The basic assumption of
this theory is that the distribution of electrons within each valley can be described ade-
quately by three quantities, namely, electron density nL , drift velocity 17}, , and tempera-
ture T L in the lower valley, and n”, V”, and TU in the upper valley. It follows that
all quantities which depend on the distribution function are uniquely determined by the
parameters nL,ITL,TL ,nU,VU, and TU. This assumption is satisfied by displaced

7

 

 

 

 

Figure 2.1. The displaced Maxwellian distribution can be uniquely described by three
quantities: n, I7 and T. n is electron density and describes the area under the distribution
curve. 9' is the mean velocity and describes how far the distribution function is displaced
from its equilibrium value 50 = 0, and T is electron temperature and describes how
wide the particle velocity spreads away from its mean velocity r70.

9

Maxwellian distributions which can be expressed as[17]

‘ e 3/2 [ m.(V°-\T)2
_ . m I I
my") " "' [zit/car] ”H 2kBT ] ' (2'1)

 

 

Equation (2.1) is described by three quantifies : n,- , r7} , and T,- . However, the analysis
is not restricted to any particular distribution function. Applications to device simula-
tion using the Blotekjaer-based transport model may be found to date in Cook and
Frey [5]. and Snowden and Loret [7].

2.2 Model Development

This section will develop the first three moment equations based on the

Boltzmann transport equation (BTE) for a single, parabolic energy band.

Charged particles in semiconductors can be characterized in terms of their posi-
tion in space r and velocity v at time t. The density of particles It (r,t) may be
described by means of a distribution function f (r,v,t ), which is itself a function of

phase and velocity space as well as time. The density of particles is given as

n(r,t) = If (r,v,t)dv . (2.2)

The BTE can be written as [17]

El EL L.§L=§L
aIHar +m" av (a:)° ‘ (2'3)

where f is a distribution function f (r ,v ,t), F represents external forces, m' is parti-
cle effective mass. The right hand side includes the randomly-timed scattering events
that the particles experience. Since (2.3) does not have a closed form solution for the
devices being considered, one approach to solving the BTE consists of simulating the
motion of one or more carriers at a microscopic level with Monte Carlo methods [1-4].
However, this category of simulations is computationally intensive, and there fore, with

a few exceptions only, not suitable for engineering application.

10

An alternative approach is to solve the moment equations, which are obtained
from (2.3) through multiplying by various functions of velocity, <D(v ), and integrating
over velocity space. Essentially, <I>(v) takes the values 1, v( or vi), vv( or vivj ), ...,
thus giving rise to the zero-order, first-order, second-order, ..., moment equations,
respectively. The procedure is straightforward, though tedious for higher orders. It
replaces an equation for the distribution function, f (r ,v ,t ), by equations which are

function of r and I only.
The average over velocity space of an arbitrary function ¢(r ,v ,t) is defined by

= 1
n(r,t)

 

<¢> [pn- ,v ,:)f (r ,v ,t)dv (2.4)

where n (r ,t) is defined in (2.2). Hence, multiplying (2.3) by ¢(v) and integrating over

velocity space, the general moment equation may be written

"e E-<fl>

8(n <<D>) +
m' 3"

a:

 

1.010 ¢>) -
8r

C

=[%......]

where F is replaced by eE . e is the particle charge and E is the electric field.

Thus the BTE for the distribution function f is replaced by a set of equations
containing averaged quantities. Each moment equation introduces the next higher-order
velocity moment due to the second term in (2.5). The moment equations are then an
infinite set of equations unless some additional assumptions are used to break the chain
of equations and restrict the variables to a manageable number. These additional

assumptions will be considered below.

The zero-order moment equation is obtained from (2.5) by putting (D = 1

11

 

 

 

 

an _ Q.
where
v=<v>=ijvfdv . (2.7)
n

The zero-order moment equation (2.6) is the carrier continuity equation. According to
(2.6) the increase of particle density plus the divergence of particles equals the

increase in density due to collisions.
The first-order moment equation, putting (b = m ' v,- , gives for the ith component

of the equation

3(m ' m7”

at + V~(m'n<vv,->) - neE; = [1(m‘n17i)] . (2.8)

a: c

The general second-order moment equation uses d) = (m ' /2)v,- vj. For the development

of the energy conservation equation it is adequate to take the cases i = j and sum over

i = 1,2,3. Then, (I) = i—m’vz, and this gives

a 1 e 1 e _
at(n<2m v2>)+V(n<2m vzv>)-neEv

= [58:01 <-;-m' v5] (29)

where <v,-vj> and <vzvi> may be written

 

<v2vj> = Eat-Viv!) = 2 [—% + f V} + :1 vi + vivivj] (2.11)

 

 

12
with pij, the pressure tensors, being given by
pi, = m‘n<(v,. - 9})(vj - VI.» = m‘j(v,. — 17,-)(vj- 71-)de = pi,- (2.12)
and qijk, the heat tensors, being given by
q... = m‘n<<v.- — Vszj - 17,-va - a» . (2.13)

Substituting (2.10) and (2.11) into the first-order moment equation (2.8) and into the

second-order moment equation (2.9) gives
l(m'nV¢)+ 2331(1). +m'anVo)-neE- = i(m‘mT) (2.14)
at ‘ i=1 arj U ‘ j t at ‘ C

and

a 1 1 3 a 1 1 __ _ y. 1_ _
E [zpii + 2’" "Vi 2+] 1.2::1'3—[2 -qiij+ Zpiivj +pijvi +m "(Evizvj)]

C

- new; = [%(n<%m‘vfi>)] (2.15)

for the momentum and energy conservation equation, respectively.
From elementary gas kinetic theory, we define
"ks Tn = m‘ to. - vim,- - v7>fdv (2.16)
where k3 is Boltzmann’s constant and Ti]- is related to the pressure tensor as
P.) = "kaTij . (2.17)
Finally, the components of the heat flux vector are defined as

_ 1 I- 2 _ _ 3 1
Q.- -3". [(v ~97 (v.- -vi)de-J}_.31-2-qijj (2.18)

by (2.13).

13

Using (2.17), (2.14) reduces to

 

' '1 3 ankT
8m nv, +2 —-B—T-L+— a -—-,-'v(m nvj) -enEi
a; '-1 ar- 3"
1— 1 r1
= a —(m m?) (2.19)
a: c -
or in vector notation
9%": + V-Wm'nTr) + V-nkBI - enE = [%(m’nfi] (220)
where
3 3
T=',EET""

I is a tensor, and a,- and aj are unit vectors. Assuming 1' is scalar, i.e. I = T_I_, then

the moment conservation equation (2.20) becomes

 

am‘nV

a:

 

+ vv(m‘m = en E-V(nkBT) + [%(m'n7)] (2.21)

C

 

 

 

or using (2.6), (2.21) can be rewritten as

 

fl +v "477+ —V(nkB T) - —E= [#2:] (2.22)
m m a:

a:

 

 

 

where the collision term was written as

av _ 1 3 . “231
[-:-],-..[..mm], ..[..],-

In (2.21), m'n'v' is the momentum density. The left-hand side is the rate of change

plus the outflow of momentum density. The right-hand side represents the force

l4

exerted by the electric field and by the particle pressure nkB T, and the rate of momen-
tum density change due to collisions.

Using (2.17), (2.18), and E = nkpl' = nkp TL, the second-order moment equation
(2.15) reduces to

a; l ._, . a .1. .4— —
at(2nkBT+ 2m nv )+V [[anBT'I' 2m nv ]v+Q+nk3Tv]
-J-E= i(n<lm'v7'>) (224)
a: 2 c '

where J = en‘v' is the current density. Using

w = %nkpr + -;-m‘nr3= n§ (2.25)
-2 .1. ._,
g- 21.,“ 2m v (2.26)
s = (W + nkBTW+ Q (2.27)
and
aw _ 3 l .
[.Jwval

(2.24) can be written as

 

3W - W.
at +VS-JE— [311: (2-29)

 

 

 

Ol'

15

 

333% + V-(VW) = env-E - v-(v-nkB T) - V-Q + [59%] (2.30)

 

 

 

where g is average electron energy and S is the flux of energy flow. These two equa-
tions are equivalent expressions of the energy conservation equation. (2.30) contains on
the left-hand side the rate of change plus the outflow of kinetic energy density W . On
the right-hand side the first term is the energy supplied by the electric field, the second
term is the work performed by the particle pressure, the third term is the divergence of
the heat flow Q , and last term is the rate of change of kinetic energy density due to

collisions.

(2.6), (2.21), and (2.30) are the first three moment equations for the conservation
of particles, momentum, and energy (see [16 equations (1), (2), (3)]). The hydro-
dynamic transport model solves the three moment equations with some appropriate

assumptions.
2.3 Relationship Between HTM and DDM

The improved accuracy of the HTM as compared to the drift-diffusion model
(DDM) can be demonstrated by showing the assumptions necessary to derive the DDM
from the first two moment equations. The DDM model does not consider the energy

moment equation as it assumes the energy is always at its equilibrium value.

For any quantity X , the relaxation time approximation for the scattering term

reads

 

3X _ ano
[., L- .

where X 0 is the value X at equilibrium and tx is the relaxation time for X quantity.

For the momentum relaxation time approximation, V = 0 at equilibrium. For the energy

16

relaxation time approximation for scattering terms, e = -q for the case of electrons,

 

        

and
I. > IV. + -v-.Vv , (2.32)
1: 8
(2.22) reduces to
T +—g-E- = -

Rearranging terms gives

k k T
(1,, (§)E + i’L—JQVT + fiéflvn (2.33)

1:
where 11,,(é) = q mfg) is the electron mobility. Applying the two central approxima-
m

 

tions of DDM to (2.33) ,i.e., (1) carrier temperature - lattice temperature equilibrium
(T = T 0), where To is room temperature or lattice temperature, and (2) stationary rela-

tion between electric field E and electron energy E, (2.33) can be reduced to

7,, = - u. (E )E + 0,, (E )—:-Vn (2.34)

k T
where D,, (E) = u" (E )-iq—. D" is the electron diffusion coefficient. In (2.34), un(§)

and 0,,(5) have been replaced by 11,, (E) and D,, (E) respectively due to the stationary

assumption.

Using (2.33), the current density expression for HTM becomes

k T
em? = q nun(§)E + ##VT + i-Sflvn . (2.35)

G:
3

Using the velocity equation (2.34) the current density expression becomes

17
J" = enVn = q (n p." (E )E + D" (E )Vn) (2.36)
which is the well known current density expression for DDM.

Comparison of (2.35) and (2.36) shows the fundamental difference between these
two models. The DDM neglects the current due to temperature gradient which may be
comparable in size with the other terms in equation (2.35). Also, DDM assumes
locally field dependent mobility which may result in large errors for submicrometer
device. The errors occur because of the time it takes for the electrons to reach station-

ary status, and if the time is compatible with the transient time of the electron travel-

ling through the channel, then the non-stationary effect should be taken into account.

Besides the BTE or the moment equations, Poisson’s equation is solved to obtain
the potential distribution. The Poisson equation for unipolar (electron) semiconductors
is

v2“; = .. .‘BLWD _ n) (2.37)
where \v is electrostatic potential, 8 is the dielectric permittivity of the material, ND is

the doping, and n is electron concentration. The electric field E is obtained directly

from the potential using the relationship
= -V‘V . (2.38)

Therefore, the difference between solving DDM and HTM can be logically
described in Figure 2.2. The comparison between HTM and DDM for the velocity

component due to electric field may result in three cases:
(1) v = u(§)E < u(E )E , non-stationary and undershoot.
(2) v = u(§)E = NE )E , stationary.

(3) v = “(QB < NE )E , non-stationary and overshoot.

18

£9 a

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  

 

>——~ .
V Solve. Poisson’s
Solve. Poisson’s equatron
equanon ;
Solve energy
conservation
equation
Calculate field i
d d t aramt
epen en p er Calculate energy
l1 (E) ’ D (E) dependent paramter
u (a). r, (a)
Solve continui Solve. continuity
equation ty equation
No
Converge ?

 

    

(a) (b)

Figure 2.2. Computer flow chart for (a) drift diffusion model, (b) hydrodynamic transport
model.

19

If the electric field varys slowly in space and time as in the case of long channel
devices operating at low frequencies, then the non-stationary effect is not important. If
the electric field varys rapidly in space and time, then the non-stationary effect is
important and the electron energy has to be determined using the energy conservation

equation.
2.4 GaAs MESFET Simulation

The three moment conservation equations based on solving the Boltzmann’s tran-
sport equation have been developed in Section 2.2 for the single valley. In the case of
multivalley semiconductors, these equations must be written in each valley. However,
the generation of transport equations for electrons in multivalleys leads to a highly
complex non-linear model. The problem is too complicated to be of any practical use
and is simplified considerably by deriving a single equivalent electron gas, whose
parameters are described by weighted averages of the electron population in each val-
ley as obtained from steady-state Monte Carlo simulations. Even though GaAs is a
three valley material, some researchers have used a two valley model [18,19]. Using a

two-valley model, the single electron gas model uses the following quantities[20]:
n = nL + nu
V = (1'17 U (§))VL + F U (§)VU
t t
m’ = (l-FU(§))mL + FU(§)m U

5 = (1": U (§))§L + Fu(§)(A§UL + 50) (2.39)

where L denotes the lower conduction band valley, U denotes the equivalent upper
valleys, FU (fi) = nU/(nL + nu) is upper valley fraction, and A501. is the energy
separation between minima of the lower valley and upper valleys. The energy refer-
ence for EL and fig are the respective valleys. The conservation equations of the single

electron gas model with the relaxation times approximation are [20,21]

20

 

i + V-(n v) = 0 (2-40)
3:
av 51E 1 Vgéz
._ .V V nk T = — .
at + v v + m'(§) + m' (§)n ( B (§)) Tm(§) (2 41)
a §—}2'k3 To
J? + V?! v(§ + kBT(§)) + V'Q = -qn V'E - "W (2-42)
with
= 'é-M' (5)"2 + %kBT(§) + F U (§)A§UL - (2'43)

The model equations (2.40)—(2.43) still need further simplification in practical use
for compound semiconductors. There is a diversity in these models caused by making

different assumptions as reviewed below.

Feng and Hintz [8] neglected heat flow Q and employed a relationship between E

and T given by
a: é-m‘vz + ‘3'](37‘ . (2.26)

The F U (§)A§UL , though comparable in magnitude with 5 when electrons become hot,
has also been neglected. This approximation is often made by authors working on the

modeling of compound semiconductor devices [7,8].
The next simplification made by several authors [6.9-11] is to neglect the time

. . dv . . .
derivative term (7) and the convective term v-Vv 1n the momentum conservation

equation. This gives

 

v = - u(§)E +

k k T
3:“) VT + ill-(31v): (2.33)

qn

The bar over v has been dropped for simplicity. Hereafter, the bar is always dropped

out.

21

Curtice and Yun [9] developed a carrier temperature model which used

1 kBT
V =-’ - [H.(T)E + :(V—q—[KTNfl] . (2.44)

They showed the simulation results using (2.44) are similar to those generated using

the simpler form

1 kBT
V = - [H.(T)E + :(TLL(T)VII)] . (2.45)

This is based on the assumption that the diffusion current in (2.44) is generally a
sufficiently small part of the total conduction current.
A further simplification of the momentum and energy conservation equations

(2.41) and (2.42) frequently used in semiconductor device models is to assume that

spatial variations are small. This reduces these equations to [5]

d(m;t(§)v) = _qE iii); (2.46)

and

 

(2.47)

These equations neglect any diffusion contribution due to Vn or VT in the momentum
conservation equation, and they neglect any spatial variations in the energy conserva-

tion equation. Hence, this model is not suited to modeling Schottky barrier devices.

The models which use FU (§)A§UL in (2.43) can be found in Cook and Frey [6],
and Curtice and Yun [9]. Cook and Frey use the thermal energy according to the equa-
tion

3

aka T = E ' F U ($135111. (2-48)

where g is the total energy and A§UL is the energy separation between the lower and

22

upper valley. This relation can be used to modify the momentum conservation equa-
tion (2.41) to obtain a model for multi-valley semiconductors, the new velocity is

given by

T»: 2 3170 (g) 2
v = -7 {41: + .5 [145m 8: 1v: + a [g - F0 (@13ng ]Vn}. (2.49)

 

In this research, (2.48) has been replaced by

3
31,7- : g _ up (2.50)

where U p is the average potential energy due to the fraction of particles in the L and

X valleys, so (2.49) will be replaced by

m

a
v = -13, {qE +-§- [1%]Vg + 3527 [g — Up(§)]Vn}. (2.51)

If we define the ratio of thermal kinetic energy to total energy r by

 

 

3
—kBT
r(§) = 2 = ELUL’Q ’ (252)
§ P.
then (2.51) can be rewritten as
Tm 2 2
v = - . {qE + —V(r (§)§) + —r(§)§Vn}. (2.53)
m 3 3n

(2.53) provides an alternative for numerical implementation by switching U p (é) to
r (§).
The heat flux term in the energy conservation equation is expressed by assuming

[15]
Q = -KVT (2.54)

with the Wiedeman-Franz relation used for thermal conductivity K[12,22]

23

2
_£§L_
cm #0:) 4 Ag” (255)

A(&) is a dimensionless number of order unity given by [12]

lags-[1:45]]

(k3 T )2

 

A(§)= (2.56)

where t(§) is the energy-dependent relaxation time of a particular type of collision
[12]. Many authors have neglected the heat flux term while other authors have kept

this term, in this work the heat flux is retained in order to study its effect.
2.5 Energy-Dependent Parameter Calculation.

The energy transport model used here requires the stationary average values of
electron velocity vs, , average total energy g, and average potential energy Up as func-
tions of the electric field. A steady state three-valley Monte Carlo simulation for
homogeneous doping was performed to get their values. The simulator uses a three
valley model for the conduction band of GaAs. Only the electron transport is con-
sidered in the simulator, hole transport is not included. The three valley model includes
the 1", L and X valleys as shown in Figure 2.3. Each individual band is treated as a
non-parabolic band with an energy-crystal momentum relationship given by

11k“

El+aE =—— 2.57
( ) 2m, ( )

where E is the energy, a is the non-parabolic factor, k is the crystal momentum and
m. is the effective mass. The values of a and m' are different for each of the three

valleys.

The simulator uses the particle simulation techniques so that each of scattering
mechanisms can be included individually to more accurately represent the physics of

the electron transport process. The scattering scattering mechanisms included are

24

 

3
5
..D
B
o
'8
s L X
C
a F

é §X

L gr
‘ 111 l i (.100)
< >
g <_ Heavy hole band

m 0
E V: Light hole band
0 Split-off band
§
3"

Figure 2.3. GaAs band structure. gr , Q and 5x represent energies of the extremum
of the I‘, L, X and split-off bands, respectively.

25

acoustic phonon scattering, optical phonon scattering, intervalley scattering, and ion-
ized impurity scattering. Each of the phonon scattering processes may occur by pho-
non absorption and phonon emission. The particular electron scattering process which
occurs is determined by a statistical process using random numbers. The average
scattering rate (# collisions/unit time) is calculated and used to determine the probabil-
ity of each scattering process. The probabilities are then used in combination with a

random number generator to determine which scattering mechanism occurs[2].

The stationary results, v” versus E, i versus E , and Up versus E , for a daping
of ND = 2x10” cat-3, and a lattice temperature of T0 = 300° K ate depicted in Figure
2.4.

The steady-state relationship of equations (2.41) and (2.42) for homogeneous case

can be written as [8]

 

 

 

 

_ In' (§)vn(§)

THE) - 45 (g) (2.58)

and
§ - £0

“(a ' qE(§)vn(§) (2'59)

Note that v,,(§) and m’ (a) can be calculated since :09) is known.
. tm(§)q . . .
Equation (2.58) can be rearranged to treat , i) as a variable u(§) glvmg [6]
m
v...(§)

The mobility u(§), the energy relaxation time t§(§), and the ratio of thermal kinetic
energy to total energy r (é) as function of energy are depicted in Figure 2.5.

(a)

(b)

(c)

v
(105mm)

3:
(eV)

p
(eV)

3.0

26'

 

2.5 —
2.0 -
1.5 —
1.0 -
0.5 -
0.0

 

 

 

 

0

l
5

I I l

I ll
10152025303540

 

0.6
0.5 -
0.4 —
0.3 -
0.2 ..
0.1 —
0.0

 

 

0

0.6

l
5

 

llllll
10152025303540

 

0.5 1
0.4 -
0.3 —
0.2 a
0.1 _.
0.0

 

 

0

l
5

 

llllll
10152025303540

Electric field E (105V/m)

Figure 2.4. Results of Monte Carlo simulation for GaAs with ND =2x1017cm'3 and
T0=300°K : (a) average electron velocity v , (b) total electron energy §, and (c) elec-
tron potential energy Up (6) versus electric field E.

(a)

(b)

(C)

t§(§)
(p8)

“(5)

(m 2V'ls'l)

r(§)

27

 

2.5

2.0 —

1.0..

0.5 -

0.0

 

 

0.6

 

 

0.5 —
0.4 —
0.3 —
0.2 -
0.1 -
0.0

 

 

1.0

 

 

0.8 -

0.6 -

0.4 _

0.2 -—

0.0

 

 

 

l I l I
l .2 .3 4 .5

Total electron energy § (eV)

Figure 2.5. Results of Monte Carlo simulation for GaAs with ND =2x1017cm-3 and
T o=300°K : (a) electron energy relaxation time t§(§), (b) electron mobility u(§), and
(c) ratio of thermal energy to total energy, r(§), versus electron total energy é.

CHAPTER 3
NUMERICAL SOLUTION OF THE SEMICONDUCTOR EQUATIONS

In general it is not possible to obtain closed-form expressions which describe
satisfactorily the operation of modern semiconductor devices. Therefore, numerical
techniques are used to solve the full set of semiconductor equations over a specified
domain.

The most common numerical techniques used to solve the set of partial
differential equations which constitute the semiconductor equation are finite-difference,

finite-element, and box-integration techniques.
3.1 Finite-Difference Method

The continuous derivatives of the semiconductor equations are replaced by discre-
tized finite-difference approximation derived from truncated Taylor series. The Poisson
equation in two-dimensions is usually discretized using a ’five-point’ difference
approximation which gives (see Figure 3.1)

Vi+l:z'-‘l’i,z _ Wi=z'-Vi-1.£‘ Vi:!’+1_‘Vi,£ __ Vi:l""l’i:g—l

 

 

a, a;_1 + b, bJ-l
a;-+-a,--1 bj+bj-l
2 2
= -% [Now—nu] (3.1)

where ND , n, and e are defined in (2.37).

The finite-difference discretization of the current continuity equation is more cru-
cial due to the electron concentration exponentially depending on the potential. For
the DDM, the Scharfetter and Gummel scheme [23] provides the electron current den-

sity as

28

29

 

 

 

i,j+lo
. . b- . .
l-1,J 1,] J 1+1,J
__9 - :r\ - <?___
ai-l ai
bj_1
0 node
- - . half-point
1,J-1C|)

I
Figure 3.1. Finite-difference mesh point notation

 

Figure 3.2. Non-uniform finite-difference mesh

30

 

ILL-+3- l

__an .1. _._._q_.-.
J l - a- {W443 [kBT(\Vi+l VJ] "t3 [kBTWt Wl+l)]} (3.2)

where a,- is the mesh spacing between it =i and x = i +1, and B(x) is the Bernoulli

function

X

ex—l .

B(x):

 

(3.3)

Equation (3.2) allows the electron density to follow an exponential variation between

mesh points. The DC electron continuity equation reads

 

 

 

V'J = qU (3.4a)
Jx,i+-;- — Jar—é- Jy,j+-% - Jy,j-;-
+ = qU (3.4b)
ai+ai_1 bj+bj-'l
2 2

where U is the net recombination rate. The subscript n for J has been neglected for
unipolar devices without any ambiguity. It is straightforward to construct the discre-

tized continuity equation (3.4b) using the discretized current density equation (3.2).

Finite-difference techniques applied to semiconductor modeling are well esta-
blished and there is considerable information available on the stability and convergence
properties of these schemes. A major disadvantage of the finite difference technique
when compared with the finite-element method, is that the requirement for a finer
mesh in one specific area of the device necessitates the presence of a finer mesh in
other parts of the device (see Figure 3.2). This leads to a surplus of nodes in regions

which could have a relatively course mesh.
3.2 Finite-Element Method
The finite-element method is fundamentally quite different from the finite-

difference method. The domain requiring analysis is sub-divided using piecewise

approximations to produce a mesh with m nodes at intersections of the mesh as shown

31

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.3. Finite-element mesh.

 

 

 

 

 

32

in Figure 3.3. The basis of the finite-element method is to approximate the solution in
each element. This method was not generally applied to semiconductor device model-
ing until the early 1970’s when Barnes and Lomax [24] reported one of the first two-
dimensional finite-element simulations which was used to model a short gate-length
GaAs metal-semiconductor field-effect transistor (MESFET). A finite—element discreti-

zation for the drift-diffusion model is discussed in Appendix A.

An important advantage of the finite-element method over the finite-difference
schemes is that there is no additional complication introduced by using elements of
different size. Finite-element methods provide a flexible means of investigating sem-
iconductor devices with non-planar geometries and regions of highly nonlinear field
and carrier distributions. The principal disadvantages of the finite-element method
compared with the finite-difference approach are that initially more effort is required to
implement the technique and that the stability and convergence criteria are not as well
understood for the semiconductor equations. Additionally it is difficult to develop spe-
cial finite-element schemes which are competitive with the exponentially fitted finite-
difference method.

3.3 Box-integration method

Box-integration methods [25,26] link exponentially fitted carrier concentrations to
the finite-element method which divides the simulated domain into elements. This is

the technique used in this investigation and it is described in this section.

The general semiconductor equation can be expressed as
VP (u ,ti) = c (u ,ll) (3.5)

where F represents any physical flow quantity like electric flux density D , current den-
sity J, or energy flux S, and u represents the unknown state variables which are being
solved, e.g. (um ,§). The divergence theorem states for a generalized function F that

the volume integral transforms into the surface integral given by the relationship,

33

 

 

 

 

 

 

P2
led}. >1
#F23
(a) P1 a; 1:13 P3 7:31—35
I"
P4

 

(b)

 

 

 

 

(C)

 

 

: = >——
PI 1:12 P2
Figure 3.4 Box-integration method based on (a) point basis, (b) rectangular
element basis, and (c) triangular element basis.

34
IV-de = I F'W (3.6)
vol surf

where 21’ is the outgoing unit vector normal to surface dA. Applying (3.6) to (3.5)

gives

[IV-2121A = [adv . (3.7)
vol

surf

In the case of two dimensions, (3.7) becomes
jF-wl = jcdxdy . (3.8)
I A

The numerical application can be understood with the following examples.

Example 1: By ad0pting nomenclature shown in Figure 3.4(a) and by treating (3.8) on
a point basis, the discretized form for (3.8) is

(F35 " F13”? + (F34 " F2947 = 63(53):). (3.9)

Example 2: By treating (3.8) on a rectangular element basis as shown in Figure 3.4(b),

the discretized form for (3.8) is
F12d12+F14d14=c1A1 for point P1,
F23d23-F12d12=c2A2 for point P2,
—F23d23 — F43d34 = 03.43 for point P3,

F43d34 "’ F14d14 = C4144 for point P4. (3.10)

35

Example 3: By treating (3.8) on a triangular element basis as shown in Figure 3.4(c),
the discretized form for (3.8) is

F12d12_F3ld13=clAl for POintPlt
FudB‘F12d12=C2A2 for pOint P2,
17316131 - F23d23 = C3A3 for pOint P3. (3.11)

Usually the F’s are either electric flux density D , current density J, or the energy
flux S. The next task to use this method for computer simulation is discretizing F

efficiently and stably. First, it is simple for D ,

‘Vl‘Wj
Dr'j = EEij = 93—h]. (3.11b)

where hij is the spacing between nodes i and j. For the DDM, the Scharfetter-Gummel
(S-G) technique is applied for J according to equation (3.2). For the case of HTM, the
modified S-G technique for both I and S will be discussed in Section 3.5.

This hybrid method, which can be understood as a finite-difference method on a
triangular element or rectangular elerrrent, has been proven to work satisfactorily for
many applications [27,28]. The box-integration method allows exponentially fitted car-
rier concentrations similar to the finite-difference method. The finite-box method is a
good choice to discretize semiconductor equation for devices with non-planar
geometries and regions of highly nonlinear field and carrier distributions. It has been

implemented in the FIELDAY [27] and PISCES [28] programs.
3.4 Scaling

Since the unknown variables ( v.2: ,§ ) in the basic equations (2.40)-(2.43) are of
greatly different orders of magnitude, it is appropriate to scale the unknown variables

and parameters to 0(1) for numerical considerations. Some scaling factors for DDM

are discussed in Selberherr[26]. In this research, the scaling factors used are

summaried in Table 3.1.

The Poisson equation (2.37) and the moment equations (2.40)-(2.43) are

transformed into the following equations for numerical simulation in scaled form. The

scaled equations are

V-(Ww) = -<Nn - n)

 

 

339$;- - V-J = 0
and
34? + v-s =J-E-nfg-(2;
with
J = l1(§)(nE + V(nT)) = -nv
s = nv(§ + T) - ATn u(§)VT
= —:-T + Up(§)
and
r2 = 8
(xozqco’Wo)

(3.12)

(3.13)

(3.14)

(3.15)

(3.16)

(3.17)

(3.18)

where the assumption (2.32) is used, and (3.17) neglects (l/2)m* v2 [6.9]. Notice that

there is no kg or q parameters shown in (3.12)-(3.17), and that the magnitudes of all

the scaled unknown variables are closed to 0(1). The discretization for J, S, and J-E

will be the topic of the next section.

3.5 Modified Scharfetter-Gummel technique for J and S

The current density J and energy flux S are expressed in (3.15) and (3.16)

respectively. This section will discuss the discretization of (i) J, (ii) S , and (iii) JtE

37

Table 3.1 Scaling factors.

 

 

quantity

x,y (position)
\ll (voltage)

n,ND (doping)
MP.) (mobility)

I, Tim, 13: (1211116)

T (temperature)

i (energy)

I (current density)

S,Q (energy flux)

v (velocity)

symbol

x0

W0

V0

value

max(x,y)
1‘3 T0
4

 

maX(Np (x .y ))

maX(u(§))

XOZ

110%
To = 3WK

 

1‘3 T0 = 4W0
quoCo‘Vo

x0

quoC W
10% = (———0'—0)\Vo

10

x0

‘0

 

 

38

for the box-integration technique. Note these discretization techniques can apply to

both rectangular and triangular meshes.
(i) Current density J discretization:

Following the work in Tang[22], J between points i and i +1 in Figure 3.5 can be
expressed as

J =u.[nE +V(nT)]
=u(nE +nVT +TVn)

= )1 [ME + VT) + TVn]

= ttT [(5 + VT)71,-n + Va] . (3.19)

Rearranging terms gives

J 1
[II (E + VT)Tn + n
= VT (a + 1)%n + Vn (3.20)
where
E
= _ . .21
a (3 )
(it

If E: -%' $15, and J are all assumed to be constant in the interval [1, i+1], (3.20)

can be integrated to yield [see Appendix B ]
1
Ti

Ti+l

 

J... = u.<%) ntB <—x.> — natB an] (3.22)

 

where

39

 

"I
W: J “’1' + 1
n, m "H 1
it a1' + 1
3 II 3 p X
xi xm xi+l

Figure 3.5. Notation for both J and S discretization. The two neighboring mesh points are
xi and xi.” where w, n and i are evaluated. xm, the midpoint between xi and xi”, represents
the location where E, J and S are evaluated.

i

Ti+l

 

xb = (ct + 1)1n (3.23)

The mesh point m is midway between points i and HI. B (—x,,) and B (xb) are Ber-

noulli functions defined in equation (3.3).
(ii) Energy flux S discretization:

S canbediscretizedin asimilarmannerasJ. Defining
§= 3T . (3.24)
2 5

the percentage of thermal energy relative to total energy is given as
,(g)=__=___‘ P =-—. (3.25)

r (é) is an energy-dependent parameter which can be calculated from three-valley
Monte Carlo simulation as described in Figure 2.5. Substituting (3.24) and (3.25) into
(3.16) gives

S = -J(1.5T,é + T) — ATn ttVT

= —tt(nE + V(nT))(-l"_iT + T) - ATn uVT
= 1&1}; + l)T(nE + V(nT)) — ATn ttVT
= -116 [T (nE + V(nT)) + Afi‘lTn VT]

= arts [Tn (E + Ao-IVT) + TV(nT)]

= -1.15T[-71-:(Tn)(E + Ao-IVT) + V(nT)J (3.26)

where

41
5 = _-_ + 1 , (3.27)

Rearranging terms gives

_5_ _. 1 .1.
$81, — (E + A8“ VT)T(Tn) + V(nT)

= VT (11 + A5‘1)-11-:(Tn) + V(nT) . (3.28)

Again, if we make the approximation that E, VT, S , 8, and AS"1 all change little in
the interval [1' , i+1], (4.31) can be integrated to yield [see Appendix B]

 

 

S... = —u..&..% i“.- [(nT)tB(-xb)-(nT).-+1B(xa)] (3.29)
In
Ti-l-l
or
= 3322.4... . _. _ .
s", —tt,,,5,,, dx (rTgh [(ang),B( xb) (ang),+,a (14)] (3.30)
In—
(rT§)i+l
where

i

Ti+l

 

x, = (a + A5'1)ln (3.31)

(iii) J -E discretization:
A major discretization problem to be solved in HTM regards the inner product

J 'E, which represents the forcing term for carrier heating. The proposed discretization

scheme for J -E is based on the simple vector relationship [29]
J.E = _V.(\VJ) + \Vv.J . (3.32)

By remembering that in the DC case, V-J = qU with U being the net recombination
rate, the second term on the right side of (3.32) can be neglected for unipolar devices

to further simplify (3.32). However, in the AC case, the second term on right side of

42

(3.32) has to be retained because it is not zero in equation (3.13). Hence, both terms
have to be included to allow DC and AC simulation even for unipolar devices. Apply-
ing the box-integration method to (3.32) is similar to the examples in Section 3.3. It is
worth mentioning that (3.32) does not involve the problem of computing the current
density J on the node; rather all the physical parameters appearing in (3.32) are either
nodal values of scalar quantities or the current density through the sides emanating

from the node.
3.6 Matrix Equation Solution
3.6.1 Review of Matrix Equation Methods

To simplify notation, the HTM equations can be expressed as

F VOW! .g)
F(\v.n.§) = F,,(\v.n.§) =-- o (3.33)

F §0IUI .5)
where F v denotes the Poisson equation, F n denotes the continuity equation and F g
denotes the energy conservation equation. Various numerical techniques may be used
to solve for (3.33) equal zero. The two principal approaches to solving the HTM equa-
tions are the coupled method (Newton’s method) and the decoupled method

(Gummel’s method) [30].

The Newton’s method is used to linearize the partial differential equations. Given
an initial guess, the solution of the non-linear equations is obtained by iterating the
matrix equation

Val-“v 3F“, an,“
a? a? 83;; AW" Fvwkrnkrgk)
n n n k =_ k k
An Fn(\v*.n .§ ) (3.34)

Agk F§(Wk,flk,§k)

 

 

 

 

 

 

.3)!

where k denotes the iteration count. The correction vector for the k -th iteration is

But an 8
33.1223
an 8§J

43

given by
AW” = Wan-1 _ wk , (3.35a)
Auk = n"+1 - n" (3.35b)
and
Afi“ = gm — g" . (3.35c)

All of the variables in the problem are allowed to change during each iteration, and all
of the coupling between variables is taken into account. Due to this tight coupling, the
Newton’s algorithm has a fast convergence. The matrix size will be 3N x3N (for 3-

coupled equations), where N is the number of grid point.

In Gummel’s method, the equations are decoupled such that each one can be
regarded as an independent equation for each iteration cycle. The equations are solved

and updated sequentially. At the k-th step, Gummel’s method can be formally written

 

 

 

as
[aFW(\.;,In",§k) J'A‘Vk = _Fv(\vk.nk’§k) (3.36a)
[317" (¢;;’nk’§k) ]-An" ___ _Fn(wk+l’nk’§k) (3.36b)
[8174 Wkgénhl’g") ]-A§" = -F§(\|I*+l.n"+l.§k) . (3.360)

Note that the most recent variables are used in equations (3.36a)-(3.36c) by making
use of equations (3.353)-(3.35c). At each stage only one equation is being linearized
and solved by Newton’s method, so the matrix A has N rows x N columns regardless
of the number of coupled equations being solved, where N is the number of grid

points.

44

Both the Gummel method and the Newton method produce a linear matrix equa-
tion of the form AX = B . This matrix equation can be solved by direct or iterative
techniques. Two direct methods are Gaussian elimination or LU-decomposition. Three
common iterative methods are Jacobi, Gauss-Seidel, and successive-over-relaxation
(SOR) methods [25]. In general the matrix size A is large due to a large number of
mesh points, hence a banded matrix LU decomposition solution technique is used to

reduce the computation time needed to solve AX = B .
3.6.2 Matrix coefficient allocation and calculation

Equations (3.12)-(3.14) govern the device transport behavior. Integrating these

equations gives

 

jwootsz) + (ND - n)]dxdy =F,, ‘ (3.37)
J-giidxdy — IV-Jndxdy = I%dxdy + F, = o (3.38)
M (7.8 _ ‘ .E _ Ho
I a: duty + II (J. n1§(§))1dxdy

= a—gtéaxay + F, = o (3.39)

where F w, F" , F g denote respectively the DC part of the integrated Poisson equation,
continuity equation, and energy conservation equation. In particular, F" and F g are

expressed as

F, = — jV-Jndxdy (3.40)

fi-fio
1§(§)

 

F; = [NS - (Jn-E - n )]dxdy . (3.41)

In this study, a rectangular mesh with variable mesh spacing was used as shown
in Figure 3.1. Since the number of grid points in the x-direction is larger than that in
the y-direction for the MESFET, the grid points are numbered as shown in Figure 3.6.

45

 

 

 

 

 

 

 

 

 

 

 

 

 

‘ >X
Y
1— .— _
. . . o . . . o o o
jp-l
jp-ny jp jp+ny
jp+1
O I O O O O O
O O O
my Zny nx-ny

Figure 3.6 . Grid number assignment for a rectangular mesh. The arrangement gives a re-
duced matrix of bandwidth (2ny+ l) for Gummel’s method and (6ny+5) for Newton’s meth-
od.

45

The rectangular mesh can be refined using a regrid generator as described in Appendix
C. Applying the box-integration method [see example 1 in Section 3.3], equations
(3.37)-(3.39) can be rewritten as follow.

[(XZV‘V)» Jp+ny " (”‘7‘")ij Jp Li?

+ [(X2VW)jp,jp-+l -' (12Vw)jp_up]d_x 4' (ND - [IL-pa; = O (3.42)

3 — — — _
[3?- ]jpdx dy- [U ij Jp+ny “lip-row)” +(ij.jp+l_ij-l,jp)dx]=0 (3.43)

an —— _
["an—té' ] dxdy+ ijp+ny Sip-npr)dy +(S'Jp1r+1 “Sip-1.19%] .

 

jp
-l(S'E- d! dy = 0 3.44

or they may be expressed as:

Fm = (3.45)
an — —
3;]me dy + F”, = o (3.45)
33‘; . Err—d? + pa.» = (3.46)
JP

where ijp, Fan’ F“? denote the DC part at grid point jp in (3.42),(3.43) and
(3.44), respectively. The at"? and 47 were defined earlier in Figure 3.4(a). The matrix
equation solution of (3.42)-(3.44) is done by using either (i) Gummel’s method, or (ii)

Newton’s method.

(i) Gummel’s method

Gummel’s method solves equations (3.45)-(3.46) sequentially as shown in Figure

2.2. The matrix equation for the DC solution of the Poisson equation can be expressed

47.

 

 

 

 

 

 

as
’31:“,l an,1 817,",1 ] Am F v.1
a‘llr 3W2 a\I’N AW Fm
' = - ' , (3.4721)
3”va ava ava . .
L 3W1 W2 3%] JLA‘I’NJ {7wa
or
ANXNXle = Ble - (3-47b)

Similarly, the matrix equation for both continuity equation and energy equation can be
easily formulated. The matrix A is a band diagonal matrix with five or less non-zero
elements in each row. In equations (3.42)-(3.44) for grid point jp , the highest grid
number is jp+ny ,and the lowest is jp -ny , therefore the bandwidth is (Zny +1). By tak-
ing advantage of this banded property, the size of matrix A can be reduced from N xN
to (2ny+l)xN. For example, if nx=100 and ny=20, then N=2000, and matrix A is
reduced from 2000x2000 to 41>QOOO. This means only 1/48 of the original memory

space is enough, and that the computation time is reduced substantially.

One thing worth mentioning is that Gummel’s method requires the addition of
one more term to the main diagonal of matrix A for the Poisson equation, otherwise
the iteration will not converge. This can be understood by rewritting the Poisson equa-

tion (3.12) as

V‘XZVW‘Q'H + (ND _ nk-t-l)

+1 _
= v-xsz‘“ + (ND — nkexp(-‘1’:—U-—‘l’:) = o (3.48)
t

k T
where the superscript k indicates the k-—th iteration and U, = Ji—lulo for the scaled

Poisson equation. The added diagonal term came from the exponential dependence of

48
n on potential.
The coefficients in matrix A are calculated by numerical evaluation. For example,

. . 3F .- F ,- (v- +Aw- )—F ,- (41-)
A019,”): 8;”: = WP JP AJP WP JP
in VIP

 

(3.49)

and

B Up) = -F(.,,,-,, (vjp) . (3.50)

The numerical evaluation for A causes longer computer time than analytical expres-
sion, but it provides great flexibility for changing device models. After A and B have
been calculated, the matrix equation is solved by the LU-decomposition method in this
study.

(ii) Newton’s method

The second approach to solving the semiconductor equations is Newton’s method
which solves all three equations simultaneously. This method provides the capability
for the sinusoidal steady state analysis which will be presented in Chapter 5, and it
allows the simulation of high current densities. However, this method requires 3N x3N
for the. size of matrix A . In order to reduce the size of matrix A , the unknown vari-

ables are arranged in the following order

 

 

 

 

 

 

 

 

 

P317,” 3F,“ BFWJ 3F,“ 1 21,1 im1
aw. an. 3&1 3§N Ag: pg
: = - , (3.51a)
a}; . a); . a); . . . . 3):" . AWN va
_5.~_ w _a_lz. ,,, 5~ A,” FM
_ 3V1 3"1 3§r 351v 4 Ag” Féflj
b a L

01'

A3Nx3NX3N><l = 33er - (3-51b)

49

The matrix A is also a band diagonal matrix with fifteen or less non-zero elements in
each row. Instead of X =[Aw1,- - - ,AWNAnl. - - - ,AnNAfil, - . . ,A§N]T which is
commonly used in Newton’s method [31], the rearranged X shown in (3.51a) has been
used. By doing so, the matrix A is more compact in bandwidth. After a careful
arrangement for matrix A to match X, the size of A can be reduced from 3N x3N to
(Gay +5)x3N . The reduced matrix A can save substantial memory space. For example,
nx=100 and ny=20 gives N =2000, with matrix A reduced from 6000x6000 to
125x6000. This means only l/48 of the original memory space is used. In order to
reduce the grid point number and the size of matrix A , ny should be kept as small as

possible for this arrangement.

There are fifteen or less non-zero terms in matrix A for each row or each column.
The non-zero terms are also evaluated numerically. The matrix equation obtained is
solved by LU-decomposition as indicated for the case of Gummel’s method.
3.7 Boundary Conditions

The above set of coupled equations are solved numerically using a box-
integration formulation on a variable spacing mesh. The solution is subject to the fol-

lowing boundary conditions.

(i) Neumann boundary conditions are applied to the free surfaces. The boundary con-

ditions are
11:0, i=0, _a§____0 (3.52)
Bun 8a,, Bun

where u” is the normal direction. This is equivalent to

E,“ = o, J,_ = o, s,_ = o. (3.53)

(ii) Dirichlet boundary conditions are applied to the elecu'odes. The boundary condi-

tions are

50
w = constant, as = nd = Nb, and §= £0 (3.54)

where n, , nd and N D are respectively carrier concentrations under the source and
the drain and the doping level.
(iii)A thermionic emission-diffusion boundary condition is applied on the Schottky

contact at the gate [32]. The boundary conditions are

 

W = _¢bn + Vgs (3.55)
J = (Mr (n we exp( :43." )) (3.56)
B 0
and
g = g0 . (3.57)

The thermionic emission-diffusion model is used to calculate the electron concen-
tration at the Schottky contact. 44),," is the barrier height, v, is the effective recom-

bination velocity, and NC is the effective density of states in the conduction band.

CHAPTER 4
DC SIMULATION OF MESFETS USING THE SEMICONDUCTOR
HYDRODYNAMIC TRANSPORT EQUATIONS

4.1 Current-Voltage Characteristics Compared with Monte Carlo.

The study of small geometry GaAs MESFETs is a important application of the
HTM simulator described in the previous chapters. The DC simulation of the MES-
FET is presented in this chapter along with the verification of the accuracy of the
HTM simulator. The DC simulation of the MESFET using the HTM method has been
compared to calculations using the Monte Carlo particle method. This comparison
serves to verify the accuracy of the hydrodynamic transport model because the Monte
Carlo method is known to include the important physics of electron transport in MES-
FETS. In particular, comparisons of Monte Carlo results. to experimental data have
shown good agreement [62]. The structure used for the simulations had L65 =0.4pnr ,
LG = 0.4m, and LCD = 0.6].tm as shown in Figure 4.1. The device doping was a
constant n-type doping of 2.0><1017cm'3 for the active layer. The depth of the active

layer was 0.1um and the Schottky barrier potential was 0.75 V.

The HTM model gives the current voltage characteristics shown in Figure 4.2.
The HTM model used had one additional simplification to the model presented earlier
in Section 3.5. This simplification was implemented to improve the convergence for
higher V05 values. The source of the convergence problem was the strong dependence
of r(§) on energy for high drain voltages as shown in Figure 2.5. The difficulty occurs

in the discretization of S given as

 

d (rTg) 1
Sm = —p,,,, 8”, dx —-("—T§)i_- [(anQiB (-xb) — (”’T§)i+13 (xb )] (3.30)
('T§)i+1

51

52

 

 

 

L LGS .L LG J. LGD J
Source l T Gate T 1 Drain
T — _ , -
K /

 

Active layer

 

 

Semi-insulating Substrate

 

 

 

Figure 4.1 Simplifid two—dimensional MESFET geometry.

53

where the ang terms are the unstable terms. The simplification made is that the r in
the ang terms is set to a constant value of rs =0.4. This value of r, =0.4 is an approxi-
mate value for the r shown earlier in Figure 2.5. The simplified energy flux expres-

sion becomes

 

d(fT§) 1
5,, = r. —u..6.. dx We): [(nT§)iB (—xb> - (am-.18 (24)] (4.1)
(r T§)i+l

The importance of the S j term to the solution of the FET currents is in the gra-
dient of the electron energy. The validity of the r, approximation can be seen in Fig-
ure 4.2 where the HTM and the Monte Carlo results are compared. The other material
parameters including r (E), u and 'cg where all determined from using the one particle
Monte Carlo simulator. The one modification made to the Monte Carlo data is an
adjustment of the low field mobility to u0=0.45m2V‘1s'l. This was done to match
low field mobility values found in the literature [7]. The one particle Monte Carlo
simulator used does not accurately simulate the low field region due to an over-
simplified acoustic phonon scattering model. The acoustic phonon model used treats
the acoustic phonon scattering process as an elastic process which is a valid approxi-

mation only at higher fields.

The grid structure used in the two-dimensional Monte Carlo simulator consists of
small rectangles as shown in Figure 4.3. The simulator allows recess gate structures
with a insulating region on each side of the gate. The simulator solves the Poisson
equation using the finite difference method and it solves the Boltzmann transport equa-
tion using the Monte Carlo particle method. The boundary conditions used in the Pois-

son equation solution are:

1) Potential equals source voltage

m-j

54

 

 

 

 

 

 

3
—-—— :HTM
2.5“ eeeefi ..... :MC
2— VGS=O.0V
...... oooe0.0....coo...00.000000000000000...-1|
..a ............
ID (A ’07”) 1.5..
Keri-0.2V
"Ha .................. a oooooooooooooooooo .0
1.1 ".4 ..............
9""
.5—
0 I I I
0 .5 l 1.5 2
VDs(Volt)

Figure 4.2. The current-voltage characteristics for both I-lTM and Monte Carlo.

55

l nyssp lnysg ”ydgl nydsp ngy

_ nxress

 

- ngx

E Source-drain grids

Insulator grids

Figure 4.3. Grid structure used in the two—dimensional Monte Carlo simulation.

56

2) Potential equals drain voltage
n-c
3) Potential equals gate voltage
e-f, g-f, h- g
4) Normal component of electric field is zero
(a-b), c-d, e-d, i-h, i-j, a-m, b-n
5) Normal component of electric flux density is constant
d-l, f-l, g-k, i-k
6) Potential equals substrate voltage
(a-b)
Note that the boundary (a-b) has two possible boundary conditions (4 or 6).
The boundary conditions for the Boltzmann transport equation solution are:'
1) Reflective boundary to particles
i-j, i-k, g-k, f-g, f-l, d-l, c-d, a-b
2) Source-drain grids
a-j, b-c.
At the source and drain grids, particles are added to these grid cells each time step to
maintain a charge neutrality.
The computation times for the Monte Carlo particle and the HTM methods were
compared. For each I-V curve with 10 points as shown in Figure 4.2, it took about 4
to 6 hours CPU time for the Monte Carlo simulation on the CONVEX C-220 com-

puter, and it took about 5 minutes CPU time on the same computer with a grid size of

60x5 for HTM.

The typical two-dimensional distributions for the potential, the electron density,

and the electron energy obtained from the HTM model are shown in Figure 4.4, Figure

57

 

 

 

 

 

                            

a w .
see
i \ \ m
assesses“
\\\\\\ w... e0
\\\\\ m. are; r

   

 

d _ _

WV)
1.515.
.7601
.0050

0.0V and
1.5V. The source is the front left region and the drain is the fi'ont right region.

The gate extends fi'om 0.6 um to 1.0 pm .

over the entire device for V5;

Figure 4.4. Potential distribution (V)

VDs

58

 

 

 

 

 

 

 

 

 

 

 

 

 

/0 \ ° .000 W")

Figure 4.5. Electron density distribution (m ‘3) over the entire device for V55 = 0.0 V and
V05 = 1.5 V. The source is the front left region and the drain is the from right region.
The gate extends from 0.6 pm to 1.0 um .

59

§(eV)

 

 

.5162 r

.3570 r

.1979 r

0.0 V and

1.5 V. The source is the front left region and the drain is the front right region.

The gate extends from 0.6 pm to 1.0 pm .

Figure 4.6. Electron energy distribution (eV) over the entire device for V65

V03 =

60

4.5 the region directly under the gate is depleted with electrons present at the bottom
of the active channel. The electron concentration plot shows a dipole domain towards
the drain end of the channel. This dipole domain occurs for the GaAs MESFET
because of the multivalley nature of GaAs [33]. The electron energy distribution

shown in Figure 4.6 has a peak energy of approximate 0.52eV. ,

A more detailed picture of the potential, electron density and electron energy
along the lower portion of channel between the source and the drain is shown in Fig-
ure 4.7-4.9. In Figure 4.8, the dipole domain can be clearly seen. For these figures the
gate extends fi'om 0.6 um to 1.0 pm .

For comparison with the field dependent mobility model, the current-voltage
characteristics for DDM and HTM is shown in Figure 4.10. The DDM and HTM are
compared using the same relationship for velocity versus electric field as shown earlier
in Figure 2.4. Figure 4.10 shows that the HTM predicts a larger current than the DDM.
The DDM calculates the mobility from the electric field resulting in a stationary rela-
tionship between velocity and electric field. The HTM calculates the mobility as a
function of electron energy where the electron energy is calculated by solving the
energy conservation equation. Since the electron can not gain energy instantly from
the electric field, the electron energy for HTM is smaller than the energy correspond-
ing to the local electric field model used for the DDM. Hence, the HTM has higher
mobility than the DDM in the device channel. This results in the electron velocity for
the HTM being higher than the stationary values used in the DDM. This higher velo-

city is called a non-stationary effect and it leads to higher current [6,7].

61

 

 

 

 

 

2
1.5 —
Voltage1 _
(V)
.5 -
0 " l I I
0 .45 .9 1.35 1.8
X-position (pm)

Figure 4.7. Potential distribution at y=0.1 pm for VGS =0.0V and VDs = 1.5 V. The
gate extends from 0.6 um to 1.0 um .

62

 

 

 

 

 

 

. 3
2.5 -
2
Electron
Density 1.5—
(1017cm ‘3)
1 z
.5 —-
O I l I
0 .45 9 1.35 1.8

X-position (11171)
Figure 4.8. Electron density distribution at y=0.1 um for VGS = 0.0V and VDs = 1.5 V.
The gate extends from 0.6 pm to 1.0 um .

63

 

 

 

 

 

.6
.4 —
Electron
Energy
(eV)
.2 —
0 I I I
0 .45 .9 1.35 1.8
X-position (um)

Figure 4.9. Electron energy distribution at y=0.l pm for VGs = 0.0V and VDS = 1.5 V.
The gate extends from 0.6 pm to 1.0 pm .

 

3
—— :HTM
2.5 a ----- :DDM
2 —‘ VGS = 0.0V

 

In (A lcm) 1.5.4

 

 

 

 

 

VGS = -0.2V
1 _
5 _____ K (Effigy .......
/’,’ ______Y_G£-_--O.2V ______
0 I I I
0 5 l 1 5

VDS (Volt )
Figure 4.10. Current voltage characteristics for both HTM and DDM.

65

4.2 Source and Drain Resistance Studies of Short Channel MESFET’s using

Two-Dimensional Device Simulators1

One of the phenomena which is predicted by the HTM is that the electrons are at
high energies as they enter the drain region shown in Figure 4.9. The high electron
energy strongly influences the behavior of the parasitic drain resistance region which
exists between the gate and the drain contact. A study was conducted to model the

parasitic resistance regions using two dimensional device simulators.

The parasitic MESFET source and drain resistances strongly influence the
MESFET’s performance as the channel length of the MESFET is reduced. The
increased influence of the resistances, RS and RD , is due to the source-to-gate length,
L65, and the drain-to-gate length, LCD , not decreasing proportionally as the gate
length is decreased. Previous source and drain resistance studies[34-40] have modeled
the source resistance and drain resistance as constants which are independent of the
applied bias. One exception to using a bias independent model is the gate voltage
dependent model developed by Byun and coworkers[41]. These previous models are .
used for circuit simulator models with resistance values extracted from measured
current-voltage characteristics. The use of a bias independent source and drain resis-
tance model neglects three effects which become important in short-channel MES-
FETs. The first effect is the increase of the electric field along the current flow path in
the parasitic resistor regions as the channel length decreases. The electric field can
become large enough to yield field'dependent mobility effects. The second effect is
the increased impact of the lateral gate depletion region as the device geometries are
reduced. The third effect is the energy relaxation of the carriers as they leave the

channel region and move into the drain resistance region. These three effects cause

 

1 This section contains a paper published in IEEE Trans. on Electron Devices, vol. ED-37,
pp. 775-780. 1990.

66

changes in the parasitic resistance of the transistor as a function of the transistor
operating bias.

This section will examine the source and drain resistance in short channel silicon
and GaAs MESFETS using two-dimensional simulator tools including a drift-diffusion
simulator, an energy transport simulator and a Monte Carlo particle simulator. This
approach allows a careful study of the transistor’s internal potential, electric field, elec-
tron concentration and currents which yields an improved understanding of the source
and drain resistance for both design improvement of the transistor and for development
of models for circuit simulation. The models used for circuit simulation typically
describe the transistor as a parasitic source resistance, a parasitic drain resistance and.
an intrinsic MESFET model. Section 4.2.1 describes the use of a device simulator to
study source and drain resistance and compares the resistor values extracted using
other source and drain resistance definitions. Section 4.2.2 discusses the bias depen-
dence of the source and drain resistance values in a silicon MESFET. Section 4.2.3
presents simulations for a submicron GaAs MESFET using an energy transport simula-

tor and a Monte Carlo particle simulator.
4.2.1 Source and Drain Resistance Models

The source and drain resistance models are divided into three types according to
the definition of the resistance. The first definition is a constant value definition of the
source and drain resistance. These models use the I-V characteristics of the transistor
to extract bias independent values for R5 and RD . The constant value definition
models include the Hower and Bechtel model[34], the Fukui model[35] and the end
resistance model [36-39]. The Hower model extracts the sum R5 + RD by measuring
the drain-source resistance at VDs —> 0 as a function of the gate voltage. The Fukui
method extracts the value of Rs - RD by measuring the current flow through the gate

to either the source or the drain. This method finds RS — RD as

67

W GS dV GD
RS " RD = .276— lDrain open - ‘21;— lSource open° (42)
The end resistance method drives a current through the gate and source terminals. The
voltage at the open drain acts as a probe which gives a voltage which is related to the

source resistance as
VD
Rs = I—' - “Rd! (4.3)
G

where Rd, is the channel resistance. Similarly, a current is flowed through the gate

and drain terminals with source open circuited to give
Vs
RD = T- - ochh. (4.4)
G

The (1 depends on the bias condition and measurement method as described in the
literature[39].

The second definition[4l] is the geometric definition where the resistance is
defined for the region from the source to the edge of the gate for R s and from the
drain to the edge of the gate for RD as seen in Figure 4.11. The sum of Rs +RD is
found by measuring the value of Rs + RD + Rd, for transistors with various channel
lengths. Rd, is the active channel resistance under the gate. The sum R5 + RD + Rd,
is then plotted versus channel length and the extrapolation of the data points to LG = 0

gives R3 + RD . This method gives a gate bias dependent model.

The third definition, which is used for this research, is the power dissipation
definition of source and drain resistance. For this definition, the resistance is extracted
from a device simulator using the simulated potential, electric field, and current data.
The source and drain resistance is found by equating the parasitic resistance power loss
in a lumped model, ngm, to the power loss throughout the parasitic resistance region.

The equation for this parasitic equivalent resistance is given by

68

 

 

 

 

 

 

L LGS LG LGD
Source I AL GSI‘_' Gate “‘TALGD Drain
_ —
VI
K / IA
n+ —vw——|rw—| l'W-l—I—W— 11+
R81 R52 n RDzl Rm

 

 

 

l

l——"Y

x Semi-insulating Substrate

 

l
|
L

 

Figure 4.11. MESFET structure showing the lateral gate depletion regions ALGS and
ALGD. The thickness of the active layer is a.

 

69

j H: M
R = R regton (4.5)

1.3...
J -E is the power dissipation per unit volume as described by Navon[42] and Adler[43].
By selecting various cross-sectional regions, this method allows the resistivity versus
position to be determined The value of RS is found using the power dissipation in the
region extending from the source to the edge of the gate. This source resistance is
further divided as the resistance from the source to the edge of the lateral gate deple-
tion region, R s 1, and the resistance from the edge of the lateral depletion region to the
edge of the gate, R32. The edge of the lateral depletion region is selected as the y
position in Figure 4.11 where the surface concentration starts to decrease because of
the gate depletion effect. Typically, a decrease in the surface concentration to 90% of
the non-depleted value is used to indicate the edge of the gate’s lateral depletion
region. Knowing the value of R51 is important for understanding how changes in the

source-to-gate spacing changes the source resistance.

The drain resistance can also be divided into two sections which extend from the
edge of the gate to the edge of the lateral depletion region, RD 2, and from the edge of
the lateral depletion region to the drain, RD 1. For small drain voltages with the
transistor operating in the linear region both R01 and R02 have validity as being part
of the parasitic drain resistance. However, when the transistor is operating in the
saturation region, the use of R02 as part of the parasitic drain resistance is not valid as
the channel under the lateral depletion region is in a current saturation condition. This
region R02 becomes a part of the effective channel length that accounts for additional
canier transit time delay. Since the R02 portion of the drain resistance is strongly
dependent on the active channel behavior, the R02 portion of the transistor will be
considered as part of the intrinsic transistor device model. The behavior of R01 will

be of primary interest in this research.

70

These three methods give slightly different RS and RD values because of the
different definitions for the source and drain resistance used in each case. To compare
the three definitions, the structure shown in Figure 4.11 has been simulated and the
three methods have been used to calculate the source and drain resistance values. The
Si MESFET structtu'e simulated has LG =1.0 um, L63 = 0.4 pm and LCD = 0.6 pm.
The active channel doping was ND = 2 x 1017 cm‘3 and the channel thickness was
0.12 pm. The simulator used was the PISCES-II two-dimensional semiconductor dev-

ice simulator[28] which solves the Poisson equation and the continuity equation.

The Hower and Bechtel method was applied to the simulated current-voltage data
to get R, +RD = 2530 n-nm (resistance for a unit width of 1 pm). The Fukui
method gave a value of RD — R5 = 695 Q—um . Putting these two results together
gives R3 = 920 O—um and R0 = 1610 Q—um. The second method of finding the
parasitic resistance is by using transistors with different channel lengths. At a gate
voltage of 0.0 volts, this method yielded a value of RD + R5 = 3450 Q—ttm. This
value is greater than the Hower and Bechtel value because the resistance of the lateral
depletion regions is included in this second method. The new method proposed in this
paper gives values at V63 = 0.0 volts and VDs = 0.01 (linear region) of
Rs =R31+R32= 1452 Q—um and RD =RDl +R02=2005 Q—urn for a sum of
R3 + RD = 3457 Q—um which agrees closely with the second method. The values of
R51: 1130 Q—um and 1201:1635 Q—ttm gives a sum ofRSI +R01= 2765 Q—wn.
This sum of R51 + Rm gives a value close to that found using the Hower and Bechtel

method and the Fukui method.
4.2.2 Bias Dependence of the Source and Drain Resistance

The values of Rs and RD using the power dissipation definition are dependent on
the gate voltage, VGS . and on the drain voltage, VDs . The gate voltage dependence

was investigated by Byun and coworkers[41]. They demonstrated using the geometric

71

definition that the parasitic resistance decreases as the gate voltage of an n-channel
MESFET is increased. This occurs because the width of the lateral gate depletion
region decreases as the gate voltage increases. The power dissipation definition of the
parasitic resistances shows a similar gate voltage dependence. Using the PISCES-H
device simulator, the gate voltage dependence was calculate for the source and drain
resistance at VDs = 0.01 volts. The results are shown in Figure 4.12 where the drain
and source resistances, R5 and RD , decrease as the gate voltage increases. Also shown
in Figure 4.12 are the values of R31 and R01. Both of these resistances show an
increase as the gate voltage increases. The increase in the resistance occurs because
the lateral gate depletion region decreases in size so that the length of the source and
drain resistance regions, R31 and RD 1, increase. For a gate voltage of -0.9 volts the
lateral depletion width was 0.088 um and for a gate voltage of 0.3 volts the lateral
depletion width was 0.052 um . The use of this variable length resistor model for R51
and R 01 is appropriate when the gate’s lateral depletion width is modeled in the active
device portion of the MESFET model.

The drain voltage dependence of the source and drain resistance can also be
determined with the device simulator extraction of the parasitic resistors. The drain
voltage influences both the drain and the source resistances. Considering the source
resistance first, the source resistance changes as a function of the drain bias if the elec-
tric field in the source resistance region is large enough that the carrier mobility begins
to decrease. The carrier mobility decrease occurs in silicon due to carrier velocity
saturation. The lateral depletion width on the source end of the gate remains
unchanged with respect to the drain voltage. The source resistance values versus drain
voltage for a silicon MESFET with LGS = 0.4 um, L6 = 0.5 pm and LCD = 0.6 pm
are shown in Figure 4.13. The mobility versus electric field expression used to give

the mobility in the simulation[28] was

72

 

 

 

 

 

 

2500
2000 —
Resistance
(II-um )
1500 -
-1.0 -0.5 0 0.5
Gate Voltage (Volts)

Figure 4.12. Source (solid-line) and drain (dashed line) resistance versus gate voltage for
a silicon MESFET using the power dissipation definition. The triangles indicate the
resistances calculated using the gate edge definition and the squares indicate the
resistances calculated using the gate depletion edge definition.

73

1500

 

 

 

 

 

 

Source
Resistance 1300 -~
(El-um )
1200 - n a
an M1 _
1 100 I I I
0 1 2 3 4
Drain Voltage (Volts)

Figure 4.13. Source resistance versus drain voltage with a gate voltage of 0.3 volts. The
triangles indicate the resistance calculated using the gate edge definition and the squares
indicate the resistances calculated using the gate depletion edge definition. The barrier
potential height for the gate was 0.7 volts.

74

 

1 1/2 46
E = .
11() [1+<s.E/v...>2] Ila , ( )

where E is the local electric field, [.10 is the zero field mobility and v”, is the satura-

tion velocity.

The drain resistance, RD 1, is influenced by both the change in the lateral deple-
tion width of the gate and by the decrease in the mobility at high electric fields. The
change in the lateral depletion width versus drain voltage is plotted in Figure 4.14
where the width is seen to change by 0.135 pm for VDS changing fi'om 0.0 to 4.0
volts. The influence of the non-constant mobility as seen in Figure 4.15 is assessed by
plotting the value of the resistance from cross-section A (as shown in Figure 4.11) to
the drain contact using the resistance definition given in equation (4.5). This plot is
constructed using cross-section A so that the lateral depletion change is absent from
the resistance calculation. The effect occurring is that the drain region resistivity
increases as the drain voltage increases. In opposition, the lateral gate depletion region
increases in size as the drain voltage increases. The two effects are opposing each
other with one increasing and one decreasing the drain resistance value RD 1- The
inclusion of both effects is seen in Figure 4.16 where the drain resistance first

decreases slightly then increases.
4.2.3 Source and Drain Resistance in Subrrricron GaAs MESFETS

The power dissipation / elecu'on heating definition of source and drain resistance
is applied to short channel GaAs MESFETS in this section. The structure used for the
simulation had Les = 0.4 um, L6 = 0.4 pm and LCD = 0.6 um. The device doping
was a constant n-type doping of 2.0x1017 cm‘3 for the active layer. The depth of the

active layer was 0.1 micrometers.

One device simulator used was a two-dimensional simulator which self-

consistently solves the Boltzmann transport equation and the Poisson equation. The

75

 

 

 

 

0.15
3
Lateral 0'10 7
Depletion
Width
Change
IALG I
(“m 0.05 __
0 l T l
0 1 2 3 4
Drain Voltage (Volts)

Figure 4.14. Change in the lateral depletion width ALGD versus drain voltage.

76

1700

 

1600-

Resistance
9‘1”" ) 1500 —

 

 

 

1400 -4
1300 F I I
0 1 2 3 4
Drain Voltage (Volts)

Figure 4.15. Drain resistance between cross secdon A (shown in Figure 4.11) and the
drain contact versus drain voltage.

77

 

1900

1800-4

Drain }
Resistance 1700-—
(Q-IIM) ‘

1600—

 

1500

 

 

 

I I I
O 1 2 3 4

Drain Voltage (Volts)

Figure 4.16. Drain resistance versus drain voltage for a silicon MESFET. The dashed
line is for V3; = 0.0 volts and the solid line is for V65 = 0.3 volts.

78

Boltzmann transport equation is solved using the Monte Carlo particle method and the
Poisson equation is solved using the finite difference method[4,44,6l]. The simulator
includes the nonparabolicity of the I", L and X valleys, ionized impurity scattering,
intervalley scattering, acoustic phonon scattering and optical phonon scattering. The
simulator calculates the solution of the Poisson equation every 5 femtoseconds and the
Monte Carlo method moves the particles for the 5 femtoseconds of electron movement.
This process is repeated with an advancing time until an accurate steady-state solution
is reached. This simulator was used so that the non-stationary velocity overshoot and
the multivalley transport effects of GaAs would be included in the resistance calcula-
tion.

Another simulator used was a two—dimensional energy transport simulator which
solves the Poisson equation, the continuity equation and the energy transport equa-
tion[15,16]. The multi-valley nature of GaAs was treated using a single continuity
equation and a single transport equation to represent all the valleys. This was done by
using an energy dependent effective mass, momentum relaxation time and energy
relaxation time. The energy dependence of these values were extracted from a one-

particle Monte Carlo simulator[45-47].

The source and drain resistance behavior versus gate voltage for small drain vol-
tages was found to have a behavior for the GaAs MESFET similar to the Si MESFET
as discussed in Section 4.2.2. The source and drain resistance variation was due to the

change in the lateral gate depletion width with respect to the gate voltage.

The drain resistance at moderate and large drain voltages requires further exami-
nation because of non-stationary effects. In the energy transport solution, the steady-
state energy transport equation given by[7]

W - W0
J-E + V-a[p.WE + V(DW )1 — — = 0 (4.7)
T§(§)

79

is used. W is the energy density, W0 is the equilibrium energy density, D is the
diffusion coefficient, )1 is the electron mobility, a is the energy transport coefficient, “cg
is the energy relaxation time and i is the average electron energy. The first term is
the energy gained by the electrons from the electric field, the second term is the tran-
sport of the energy and the third term is the dissipation of the electron energy to the
lattice. The drift diffusion model used earlier for the silicon simulations assumes that
the energy is dissipated in the same region as the electrons gain energy from the elec-
tric field. The MESFET operating with a moderate drain voltage has a large transport
of energy from the active channel region into the drain region. The total power dissi-
pated in the drain region can be atuibuted to the energy transported to the region and
the energy dissipated in the region which was created by the J -E electron heating in

the region. The total power dissipation is written as
Pm = P“ + Pd, (4.8)

where I", is the power dissipation in the drain region due to energy transported from
the active channel region and P», is power dissipated in the drain resistance region
due to electron heating in the drain resistance region. An assumption is made that the
drain contract regions are ohmic contacts with an equilibrium electron density and an
equilibrium electron energy density. The design of the drain resistance region should

look at minimizing the value of Pd, .

The calculation of the total power dissipation in this region for the energy tran-

sport simulator is done using

W - w
Pm, = j ———°dr3. (4.9)

Rm region 14;)
The calculation of the electron heating in this region is given by

Pg: [ J-E M. (4.10)

Rm region

80

The method for calculating the total power dissipation in the resistance region using
the Monte Carlo simulator is to calculate the net phonon emission in this region per
unit time. The total power dissipation is given by

Eph [np‘h -' ugh]
t

 

(4.11)

Ptot=z

pk
where Eph is the phonon energy, n5), is the number of phonons emitted, ugh is the

number of phonons absorbed and t is the simulated time. The summation is performed

over the phonons of different energies produced by the various scattering mechanisms.

The drain resistance R01 is best understood by considering the power dissipation
and the electron heating. Figure 4.17 plots both of these quantifies for drain voltages
from 0.1 to 2.0 volts. The main feature to note is that because energy is transported
from the channel region to the drain resistance region, the power dissipation is greater
than the electron heating in the drain resistance region. The smaller value for electron
heating is due to the electrons moving across the drain resistance region as the result
of a gradient in the electron energy or temperature. The influence of the energy or
temperature gradient is demonstrated by considering the velocity expression for elec-

trons as given by[7]

v‘= - u(§)E + 5%?in + fl(i‘félvn (2.33)

where é is the electron energy, T is the electron temperature, It is the mobility, and n
is the electron concentration. The first term is the drift transport term, the second term
is the temperature gradient transport term, and the last term is the diffusion transport
term. Figure 4.18 demonstrates this large carrier temperature gradient between the
gate and the drain. The drain resistances calculated using R010“) = Pm, ”02 and
RD 1(eh) = Fe}, ”Dz are shown in Figure 4.19. The drain resistance due to electron heat-

ing shows a slight decrease initially as the drain voltage increases because of the

81

300

 

Power

(“W/W")

 

 

 

 

I l
0. l 0.5 1.0 1.5 2.0
Drain Voltage (Volts)

Figure 4.17. Power dissipation (dashed line) and electron heating (solid line) in the drain
resistance region Rm versus drain voltage for a L5 = 0.4 um GaAs MESFET.

82

 

2700

2100-

Temp.
(Kelvin) 150° ‘

900-

 

 

300

 

I T
0 0.5 1.0 1.5 1.8

Position (11m)

Figure 4.18. Electron temperature versus position between the source and the drain as
calculated using the energy transport simulator. The values are plotted at a depth of 0.08
pm. The gate extends from 0.6 pm to 1.0 pm. The drain voltage was 1.0 volts and the
gate voltage was 0.3 volts.

83

 

3200

2400 -
Drain
Resistance 1600 —
Q—ttm )

 

 

 

 

I I
0.1 0.5 1.0 1.5 2.0
Drain Voltage (Volts)

Figure 4.19. Equivalent drain resistance versus drain voltage for a LG = 0.4 pm GaAs
MESFET. The resistance values for the dashed curve are calculated using the total
power dissipation in the drain resistance region. The resistance values for the solid curve
are calculated using the electron heating in the drain resistance region.

84

movement of the lateral gate depletion region. However, for larger drain voltages the
electron heating resistance increases because the electric field increases resulting in a
reduced mobility. The RD w“) resistance which includes both the electron heating

energy and the transported energy increases monotonically with the drain voltage.

The two-dimensional Monte Carlo particle simulator provides a more detailed
description of the drain resistance region since a multivalley band structure is included.
The electrons are moving from the some to the drain at moderate or high drain vol-
tages scattering into the L and X valleys. When these electrons enter the drain resis-
tance region the potential energy associated with being in a higher valley must be dis-
sipated. The dissipation of the energy can be observed in Figure 4.20 where the total
energy of the electrons above the I‘ valley minimum and the kinetic energy of the
electrons are both plotted versus position between the source and the drain. This plot
shows the loss of the higher valley potential energy as the electrons move across the

drain resistance region.

 

85

 

0.24

0.18-

Electron
Energy 0. 12 —
(eV)

0.06 —(

 

 

 

 

I I I
0 0.4 0.8 1.2 1.6

Position (11m )

Figure 4.20. Total electron energy above the F valley minimum (solid line) and the
electron kinetic energy (dashed line) versus position between the source and drain at a
depth of 0.08 pm . The gate extended from 0.5 urn to 0.9 pm . The drain voltage was 0.5
volts and the gate voltage was 0.3 volts. Calculations were done using a Monte Carlo

particle simulator. 1

CHAPTER 5
AC SIMULATION OF MESFETS USING THE SEMICONDUCTOR
HYDRODYNAMIC TRANSPORT EQUATIONS

5.1 AC Simulation Methods

Simulation programs providing DC and transient solutions to the hydrodynamic
transport equations in two dimensions can be found in the 1iterature[7,9]. The third
device operating mode, namely small-signal AC operation, has received less attention
to date in the context of numerical device simulation. This is due partly to the need to
have an existing DC device solution upon which to build the small-signal AC analysis
and partly to the nature of the AC computation. This can be understood by knowing
that HTM is much more complex than DDM, even for the DC solution. Additionally,
the HTM becomes more complex for compound semiconductor devices. This chapter

will develop an AC model using the HTM for the GaAs MESFET.

The techniques for AC analysis of semiconductor devices using the DDM can be
found in the work of Laux[48]. Three standard approaches are Fourier Decomposition
(FD), Incremental Charge Partition (CP), and Sinusoidal Steady-State Analysis (534
). This section will briefly review the first two methods, the third method will be dis-
cussed in the next section.

(i) Fourier Decomposition of Transient Excitations (FD)

The FD method applies a step perturbation AVj to terminal j about a DC operat-
ing point as shown in Figure 5.1. The FD method gives the small-signal admittance

matrix component 1751-(03) as

86

87

 

 

 

 

 

 

 

 

 

V1“) 1'10)

“VJ(0)+AVJ- I A I,(ee)

0 t: 10— Device .5 0 >
av (0, Av 0).
Aoév 0) —-l——_ 2 *4” (1)2

Figure 5.1. Fourier decomposition method.

88

F{i,~(:) ‘ 11(0)}
17' (0))

i1 = ‘
F{vj(t) — mm}

_I,-(°°)-Ii(0) fl ,
- AVJ- + AVjF{lt(t)-1t(°°)} (5.1)

 

where F is Fourier transform operator, 1,- (0) is the steady-state current, V] (0) is the
steady-state voltage, v-(t) = VJ-(O) + AVju(t), and ii(t) is the transient response at ter-
nrinal i. Upon separating (5.1) into real and imaginary parts, the conductance and

capacitance matrix entries become

 

It‘(°°) "' I;(0) a) a ,
G," = + "'—"' [1i (t) - Ii (°°)] smwtdt (5.2)
1 AV]- AVj {
and
1 °° .
Ci} = TV]. g [tr-(t) - I.-(°°)] cosmtdt . (5.3)

The FD method needs DC and transient device simulation capabilities. It requires that
limitations be placed on the step At used in the transient solution in order to reduce
the error for high frequency admittance, and it requires that upper and lower limita-
tions be placed on the exciting voltage AV} in order to avoid harmonic generation and

to dominate numerical noise, respectively.

89

(ii) Incremental Charge Partitioning (CP)

Figure 5.2 shows the incremental charge partitioning method. The CP method

finds the capacitance and conductance matrix components as

A].
GU = 27‘? v, = const, k a: j (5.4)
and
AQ-
c,,- = 37311,, = const., k a: j (5.5)

I

where AQ; is the incremental charge associated with terminal 1', AI,- is the incremental
current at terminal i, and AV; is the incremental voltage applied at terminal 1' . The
CP method requires a DC solution only. The total charge Q is calculated by integrat-
in g electron density n over the P-region or N-region for a diode, or over the whole
simulated region for a MESFET. The AQ is obtained by subtracting two Q ’5 found
from two-successive V03 biases. The CP method provides only quasi-static (low fre-
quency) admittance. A disadvantage of this method is that the CP method depends on
insight into the physics of device operation in order to partition the incremental charge
to each contact. For an N-terminal device such physical insight is rarely available.

5.2 Sinusoidal Steady-State Analysis ( S 3A ) using Hydrodynamic Transport
Model

5.2.1 Model Development

The third technique for the AC analysis of semiconductor devices is the
Sinusoidal Steady-State Analysis ( s34 ) method. The 534 method works directly in
the frequency domain, and it requires DC and AC device simulation capabilities. The
admittance matrix obtained for the 53A method is rigorously correct as was the case
for the FD method. However, two characteristics of the S 3A make the technique

extremely accurate: (1) no time steps need be selected using the S3A and hence no

 

 

     

 

 

 

 

 

AQj
7

. P a)
J r

/@

(a)

Vj (0) +AVj
a"! I l I X ----

 

 

 

 

(b)

Figure 5.2. Charge partitioning method for (a) a two terminal device, and (b) a three
terminal device.

91

inaccuracies related to time discretization occur, and (2) the device is linearized about
the DC operating point so that harmonic generation within the device is precluded.

The S 3A method is performed after a DC solution has been achieved. Starting
from a DC bias condition, an input of given amplitude and frequency can be applied to
a device structure from which sinusoidal terminal currents are calculated. Then using

the relationship
+j0)C,-- = -..‘—, I7, = 0, k #j (5.6)
1'

the frequency dependent admittance matrix, and hence conductances and capacitances,
can be calculated. Note also that by varying the frequency and examining the various
device admittances, the current, voltage and power gains versus frequency can be
directly determined.

Special numerical techniques are often adopted for the S 3A method. The Poisson
equation, electron continuity equation, and energy conservation equation are expressed

as [see Section 3.6.2]

F (1(th ,fi) = 0

an

F,.(\II,n.§) + 3741 d)’ = 0

Frown .t) + 3375-4747 = 0 (5.7)

where (Ix— d? exists due to the integrated form for F W, Fmand F g.

The AC system is obtained by substituting time dependent functions of the form
C(t) = Co-l- Ce!“ into (5.7), where C: v, n, or §, and the 0 subscript denotes a
steady-state solution for the device. Performing a Taylor’s series expansion on (5.7),

and keeping only the linear terms, we have

an,
(Fw)dc + {TL} = 0

92

ar C warm-—
(Fn)dc+TC—dc +j dyn—

BF __
(F§)dc + {-35- }ch + jtodx dy(no§ + gait) = 0 (5.7b)

where (F Wu. , (F ,, )dc and (F 94‘. are expressions for the DC solution and therefore
should be zero when the DC solution has been reached. Hence, the AC solution

becomes

 

and

[381% ]dct+1wdxd)’(noa+§on)- (5.70)

For numerical implementation, the AC system (5.7c) at point 1' becomes

r a

 

 

 

 

3% 3"} 351' " .
N 3F“ 3F“ +'(od—x_ 3F“ Y] - 0 5 8
IE] 33ij aan, J dy a 351' gj - ( . )
F5. FE. . __ FE. . __ J
.310 3"} flmgodx‘ty 351' flamedx dyme

where F w,- denotes the DC portion of the Poisson equation solution in integrated form
at point 1', similar interpretations hold for F m- and F vi . In particular, all the terms
without any to frequency dependence form the Jacobian matrix which is available from

the DC solution done to get the operating point.

 

93
After assembling the global AC matrix, the AC system becomes
[J +jD] )2 =3 (5.10)

where J is the Jacobian matrix, D contains the contributions of the time derivative
terms in (5.7) to the matrix, B is a real vector dependent on AC driving voltages and
J? is the AC solution vector given by if = [itLfiJ-Ej ]T.

For the evaluation of (5.10), recall that the arrangement of unknown variables in
the DC case discussed in Section 3.6.2 gives the matrix J a size of (6ny+5)x3N . D
in (5.10) contains two non-zero diagonal bands, hence the reduced size of D is 2x3N.

Usually D is very small compared to that of J .

Equation (5.10) is a complex matrix equation. By splitting the AC solution vector
)2 into real and imaginary parts XR + 1X1, (5.10) can be written using only real arith-

metic as

8-?][§:J=[8]-

The system order doubles compared to the DC case. Computation time for the direct
factorization of the AC system can be excessive, but the computation time can be
reduced by using a block-SOR indirect solution techniques. The block-S OR (
successive-over-relaxation ) solution procedure alternates between the following two

equations until convergence is obtained:
x,y”1 = (1 - (DR )XR" + toy-Rm," + B) (5.12a)
X,"+1 = (1 — toR)X,“ + mRJ‘1(—DXR"+1) (5.1215)

where the superscript k denotes iteration number and (DR is the SOR parameter. Laux

[48] set (0R < 1 when using f 2 f7 / 10 for the DDM.

94

5.2.2 Boundary Conditions

The AC boundary conditions are based on the boundary conditions use for the
DC solution. The AC boundary used are as follows:

(i) Neumann boundary conditions are applied to all fi‘ee surfaces. This AC boundary
condition carries over directly from the DC simulation, hence this boundary condi-

tion has been already set in the Jacobian matrix J.
(ii) Dirichlet boundary conditions are applied to the source and drain contacts,
0:17”, ii,=;id=0, and§=0 (5.13)
where V” is the small signal directly applied at the contact. ii, and iid are respec-
tively the small signal carrier concentrations at the source and the drain.

(iii)The Schottky contact at the gate is described by the boundary condition
ii;=0and i=0. (5.14)

Here, because a thermionic emission diffusion model is used to calculate the elec-
tron concentration at the Schottky contact, the AC electron concentration at the

gate carries over directly from the DC simulation.
5.2.3 Small Signal Current Calculation

After solving equations (5.12a) and (5.12b), the small signal current can be calcu-
lated by using the small signal variables at each node. For a given branch between
nodes i and 1' +1 as shown in Figure 5.3, the current density is expressed in terms of

conduction current and displaced current as
i = im, + id“, (5.15)

where

~ 8]
Jcond = frag—tit C: = \Vt',\l’i+l."i,"i+l.§i 311d §i+l (5.16)

95

 

8’; Vi +1

51' j "t +1

gt l §i+ 1

3 it t p X
xr xln x1+1

Figure 5.3 Notation for small-signal current density calculation.

 

 

 

 

 

 

 

 

 

 

 

l . 2
C O
+ +
1L V2 ylz I71 yz1 '_1—
‘" — — t7
V‘ Yll o o Y22 2
' o o
3 3

Figure 5.4. The three-terminal y-parameter equivalent circuit.

96
and

id“? =jt0£§ =1.“ [Wt-:m]

 

 

[ (ill‘q’mhmg ] .[ (i’i-{Vi+1)real ]
=- we +1 are

5.17
h h ( )
where h is the spacing between nodes i and t' +1.
5.2.4 Y-parameter Calculation
The frequency dependent admittance matrix Y can be calculated as
l7..=—"—,I7,,=0,k¢j. (5.6)

Applying a small signal at one terminal, and only one, will provide the frequency
dependence of one column in the Y matrix. For an N -terminal device, N -1 small sig-
nal excitations are required to determine the Y matrix at a given DC bias point. For
example, two excitations is needed for a three terminal device. For this three terminal
device, the first step is to apply a small signal perturbation at terminal 1 to get y 11,
and 1’21- The second step is to apply a signal at terminal 2 to get Y12 and yzz. The y-

parameter equations can be expressed as

1'l _ 113’12 ‘71 .
1.1-1:21-10 »

The equivalent circuit is shown in Figure 5.4. Once the Y parameters have been
found, other equivalent parameters, e.g., hybrid parameters: H, scattering parameters:
S, impedance parameters: 2, etc. , can be calculated using existing parameter transfor-

mation tables[49].

97

5.3 HTM Y-parameters Compared with Monte Carlo and Drift Diffusion Models

The structure used for the Y-parameter simulations had LGS =0.6ttm ,
LG = 0.4um, and Lao = 0.6m [See Figure 4.1]. The device doping was a constant
n-type doping of 5.0x10160m ‘3 for the active layer. The depth of the active layer was
0.22pm . The Schottky barrier potential was 0.75 V.

A complete set of Y-parameters versus fiequency at VDS=1.5V and VG=-0.1V
has been calculated using the S 3A method. The seven frequencies are fi'om 24.4 GHz
to 170.8 GHz in steps of 24.4 GHz. Figure 5.5 presents a comparison of this data and
data obtained with the Monte Carlo particle simulator. In Figure 5.5, at low frequen-
cies all the imaginary parts of the y-parameters are small and at higher frequencies all
the y-parameters depart from the Imag(y)=0 axis. The Fourier-decomposition method
was used to generate the y-parameters versus frequency in the Monte Carlo simulation.
There is a general agreement in this comparison. All four HTM y-parameters have the
same behavior as the Monte Carlo results. The Y-parameters versus fiequency for the
HTM and the DDM are shown in Figure 5.6. This Figure shows that the HTM has a
larger Re(y 21) than does the DDM. Also, the ya and yzz for DDM go different direc-
tions with those for the HTM. This means that the DDM loses accuracy for short

channel length devices.
The complete y-parameter set permits many useful device attributes to be calcu-
lated such as cmrent gain A,, voltage gain AV[48], and unilateral power gain GU [49].

This is done as follows:

 

 

 

I I
A,(tn)= 5:] (5.19)
I I
AV((0)= ': :1 (5.20)
I - 12
GU ((0) = f 1’21 1’12 (521)
4 [RCCYll)RC()’22) - RcomReon)]

 

98

 

3 — —u-——:HTM(S 3A )
. . . . '0 ..... :MC

 

 

 

 

 

 

 

I I I I I
4 -3 -2 -1 o 1 2 3 4
Realm

~

Figure 5.5. A comparison of y-parameters versus frequency for both the HTM method

and the Monte Carlo method. The seven frequencies vary from 24.4 GHz to 170.8 GHz
in steps of 24.4 GHz.

99

 

 

 

 

 

 

 

4
11
I
3 - —-——:I-ITM
- - + - w’DDM
21
1-
Y22
2'22
11112180) 0— ylz
-1 _.
9’2
-2_
-3 _
'4 I I I I I r
-4 -3 -2 -I 0 l 2 3 4
R0810)

Figure 5.6. A comparison of y-parameters versus frequency for both the HTM method
and the DDM method. The seven frequencies vary from 24.4 GHz to 170.8 GHz in steps
of 24.4 GHz.

100

 

 

 

 

 

4o
30..
20—
Gain
10—
(dB)
v
0 ‘AGU
-10— AI
'20 I I IIIIIII FfiTIIWI
101 102 103

Frequency (GHz)

Figure 5.7. The three gains A], Av and GU versus frequency from the HTM. The eight
frequencies are 12.2, 24.4, 48.8, 73.2, 97.6, 122, 146.4 and 170.8 GHz.

101

40

 

30—1

20 -—

Gain
((13) 10 fl \

-10— A]

 

 

 

 

IIIIIIIII I IIIIITI
102 1

Frequency (GHz)

Figure 5.8. The three gains A], Av and GU versus frequency from the DDM. The eight
frequencies are 12.2, 24.4, 48.8, 73.2, 97.6, 122, 146.4 and 170.8 GHz.

102

where ReO'n) expresses the real part of y“. The unity-current-gain frequency f7 and
the maximum frequency of oscillation f m are both figures-of-merit that are correlated
to the microwave and millimeter-wave performance. f T and f mu are the frequencies

at which equations (5.19) and (5.21) are equal to unity, respectively.

In Figure 5.6, the HTM gives a transconductance gm = 120.28 mS/mm, and the
DDM gives gm = 41.44 mS lm for low frequences. These values came from the real
part of y21. Similarly, the gate capacitance c8 at low frequencies can be calculated
from the imaginary part of y“. Figure 5.7 and Figure 5.8 shows the three gain values
versus frequency for the HTM and DDM models, respectively. The eight frequencies
are 12.2, 24.4, 48.8, 73.2, 97.6, 122, 146.4 and 170.8 GHz. By setting current gain

|)’21'
I1’11I

 

A I (to): = 1 [48], which corresponds to OdB on the plot, the HTM gives a unity

current gain frequency of fT=45.98 GHz and the DDM gives fT=l7.1 0112. By set-
ting the unilateral power gain 0(0)) = 1, the HTM gives f mu = 154.95 GHz, and the
DDM gives f m = 40.49 GHz. This shows that the DDM underestimates fT and

f m, and that the nonstationary effects are important for submicrometer devices.
5.4 AC Performance of GaAs MESFETS

The AC performance of MESFETS depends on the geometric structure and the
bias. This section applies the HTM AC simulator to the study of these dependencies.
The gain and frequency performances of microwave and millimeter—wave transistors
are generally specified in terms of GU (5.21) and f max. f m , the frequency at which
GU is unity, is a particularly important figure of merit as it is the maximum frequency
of oscillation. It indicates the boundary between an active and passive device. Com-
mon practice has been to estimate f m using [extrapolation of the microwave fre-
quency gain measurements at — 6 dB per octave. However, some works [50-53] have

pointed out that parasitic resistance and capacitance cause the unilateral power gain of

103

FETs to roll off at a -12 dB per octave. Hence, the extrapolated f m value are con-
siderably greater than the actual f mu value.

The gain slope changes are explained in the work of Steer [51] as follows. The
circuit model of a common source MESFET is shown in Figure 5.9. The best possible
gain and frequency performance of the transistor will be obtained when the parasitics
are negligible, so Steer considered the performance of the intrinsic transistor alone
which has the unilateral power gain

8 :ORDS 1

G ___ (5.22)
U 4CGS R 1 (C03 - CDC 8»: ORDS ) (020'? 2‘02)

 

 

where

2 (R 3CosXCoc + Cos )2 + CDCngRDStzlz

 

P (5.23)
(CDCngRDS-CGS)
At frequencies much less than (2rtR1CGs )’1 and ignoring CDC, GU reduces to
2
8 as
GU = —”-w%— . (5.24)
4032C G SR I

The commonly used expression for GU, (5.24), rolls off at -6 dB/octave because of
1/002 term. However, with CDC in (5.22), there is an additional -6 dB/OCtave roll—off at
high frequencies due to the complex pole pair contained in the 1/(l-p 2(02) term. The
complex conjugate poles are at the frequency

1
21tlpl°

 

fp = (5.25)

For millimeter-wave transistors this pole frequency is typically below f m so that the
pole has a limiting effect on frequency performance [51].

104

 

 

 

 

 

 

 

 

 

 

INTRINSIC TRANSISTOR
R0 CGD RD
o—Wfi 4] IL new—-
0 _L+ D
Cos _ V1
Rpsi
RI CDC GMVl _]_CDs
T
(a) ,
i Rs
0
S

 

(b)

 

 

 

 

Frequency

Figure 5.9. (a) Circuit model of a MESFET. (b) The extrapolation of the low-frequency gain
overestimates the unity power gain frequency, fm.

105

One approach to estimating f m is to use an equivalent circuit model developed
from physical insight as well as 2-port measurements to determine the value for each
model component [51]. However, errors can be made due to the inaccuracy of the
equivalent circuit model, due to uncorrected parasitics in the measurement and due to
the extrapolation of low frequency measurements to high frequencies. Another more
accurate approach to estimate f m for short channel FETs is to extract the y-
parameters from a device simulator using the HTM directly in the frequency domain.
In this case, errors due to the equivalent circuit model and due to the quasi-static

assumption can be precluded.

Using the HTM model the frequency behavior of the MESFET with various
geometric and bias parameters have been simulated. In particular, investigations of the
G U versus frequency dependence on device parameters have been studied. The device
parameters studied include ( 1) drain bias dependence, (2) gate bias dependence, (3)
gate length dependence, (4) gate-source spacing dependence, (5) gate-drain spacing
dependence, (6) epilayer thickness dependence, and (7) substrate dependence. The
nominal device has LGS =0.4tlm, L6 = 0.4tlm, and LCD = 0.6m as shown in Figure
4.1. The device doping was a constant n-type doping of 2.0x1017cm’3 for the active
layer. The depth of the active layer was 0.1um and the Schottky barrier potential was
0.75 V.

5.4.1 Drain Bias [Dependence
Figure 5.10 shows GU versus frequency for VD = 1.0, 1.5 and 2.0V. The f max

determined from Figure 5.10(a) is shown in Figure 5.10(b). This maximum f m is

near VDS=1V°

106

5.4.2 Gate Bias Dependence

Figure 5.11 shows GU versus frequency for VG = 0.0,—0.2 and —0.4 V. The f M
determined fiom Figure 5.11(a) is shown in Figure 5.11(b).
5.4.3 Gate Length Dependence

Figure 5.12 shows GU versus frequency for LG = O.2,0.4 and 0.6m. The f m
determined from Figure 5.12(a) is shown in Figure 5.12(b). Some features can be
observed: (1) the f mu decreases with increasing L6, (2) the slope changes from -6
dB/ocrave to - 12 dB/octave particularly for the long channel device. The unity power
gain frequency using a -6dB/octave extrapolation over-estimates the f m and results in

a larger error for the longer channel device f m than for the short channel device.
5.4.4 Gate-Source Spacing Dependence

Figure 5.13 shows 6,, versus frequency for Log = 0.2.0.4 and 0.8m. The f m
determined from Figure 5.13(a) is shown in Figure 5.13(b). The f m decreases only
slightly with increasing LC. For this particular structure the parasitic source resistance

is seen to have only a small effect on f m.
5.4.5 Gate-Drain Spacing Dependence

Figure 5.14(a) .shows GU versus frequency for LCD = 0.2.0.4, and 0.8m. The
f m determined from 5.14(a) is shown in Figme 5.14(b). The longer LCD value gives

a smaller drain-to- gate feedback capacitance, hence, a longer LGD gives a higher f max.
S.4.6 Epilayer Thickness Dependence

Figure 5.15(a) shows GU versus frequency for a = 0.1 and 0.12pm. The smaller
epilayer thickness gives the larger f mu because of the fighter control by the gate vol-
tage VGS . Since a non-ideal substrate permits both current flow and a reduction in

fields at the substrate/epilayer interface, the sensitivity of performance to epilayer

thickness in a real device will be somewhat less than that presented here [11].

107

5.4.7 Substrate Dependence

Figure 5.l6(a) shows GU versus frequency with substrate and without substrate.
The doping profile with substrate is shown in Figure 5.l6(b). The bias point is

Vos = 1.0 V and VGS = -0.2 V. The device with substrate has a lower f mu.

Table 5.4.1 shows the various simulation with different dependence, the f max

corresponding to each case is also included.

108

 

 

 

 

 

 

 

 

 

40
——-u— VDS = 1.0V
——a-— V05 = 1.5V
30—? —e— V05 = 2.0V
20_
Gain f‘
((13)
10— ‘
\
o .‘ fir x
-10 l IIIIIIIT I 7 [Tllll
101 102 1
Frequency (GHz)
(8)
450
400..
350—
f max
(GHZ)
300-
250—
I
200 I I I
0 .5 l 1.5 2
VDS (Volt )
(b)

Figme 5.10. (a) GU versus frequency for various VDs . (b) f m versus V05.

109

 

 

 

 

 

 

 

 

 

40
—-u— V05 =-0.0V
—B— V05 =-0.2V
30 " —e-— V65 =-0.4V
20-4
Gain “
(dB)
10— “
O ’\E ‘
'10 I IIIIHII I 711111]
101 102
Frequency (GHz)
(8)
350
300-
f max
(6112)
250— (
200 I I F I
.05 0.4 03 -02 -0.1 0
V65 (Volt)
03)

Figure 5.11 (a) GU versus frequency for various VGS. (b) f m, versus V05.

110

 

 

 

 

 

 

 

 

 

40
—a— LG =0.4Il.m
30 —e— LG =0.6|,lm
20..
Gain
((13)
10—
° \
-10 j llllllll I lllllll
101 102 103
Frequency (GHz)
(8)
300
280..
260-4
f max
(GHz)
240-
220- I
20° F I I T
0 .2 .4 .6 .8 1
L603")
0!)

Figme 5.12 (a) GU versus frequency for various La. (b) f m versus Lg.

111

 

 

 

 

 

 

 

 

 

40
—~—— Los =0-7—llm
—9— L05 =04le
30 —e— L05 =0.8l.lm
20—1
Gain
((13)
10—
0
-10
I llllllll I lllllll
101 102 103
Frequency (GHz)
(a)
300
230— \
260..
f max
(GHz)
240.4
220-
200 I I I I
0 .2 .4 .6 .8 l
LGSOJM)
Cb)

Figure 5.13 (a) GU versus frequency for various L05. (b) f max versus L55.

112

40

 

Gain
((13)

 

 

 

 

 

 

 

 

 

'10 I IIIIIIII I IIITIFI
101 102 103
Frequency (GHz)
(a)
300 t
280—1
260—
f max
(GHZ)
240—
220-
200 I I I I
0 .2 .4 .6 .8 1
LCD (11"!)

(b)
Figure 5 . 14. (a) GU versus frequency for various Lap. (b) f max versus LCD.

113

 

 

 

 

 

 

 

 

 

4o
—-9— a=0.lum
__...__ a=0.12wn
3o.£
20-
Gain
((13)
10—
0 \r In
'10 I I IIIIIII I IIIIIII
101 102 1
Frequency (GHz)
(8)
300
280_
250_
fnwt
(GHz)
240-
2.20—
20° I I I I
0.06 0.08 .1 .12 .14 .16
awn!)
(b)

Figure 5.15. (a) GU versus frequency for various epilayer thickness a. (b) f max versus a.

114

 

 

 

 

 

 

 

 

 

 

40
—e— a =0.l},tm
—~— a = 0.1111): +substrate
30-4n
1
20—
Gain
(43)
10—
'10 I IIIIIIII I I IIIIII
101 102 1
Frequency (GHz)
(a)
4><10l7
3x1017—
(CZQQ) 2X1017
1x1017d
15
1.0x10 I I I I I I
0 0.05 .l .15 .2 .25 .3
Depth (Hm)
(b)

Figure 5.16. (a) GU versus frequency with or without substrate. (b) ND versus depth.

Table 5.1 fum dependence on device parameters for ND = 2x1017cm'

115

3

 

L003“)

145501-31)

1431304111)

mun)

VDs(V)

VosW)

fmu(GHz)

 

0.4

0.4

0.6

0.1

1.5

-O.2

278.8

 

0.5

207.2

 

1.0

416.5

 

2.0

223.7

 

-0.0

246.9

 

-0.4

346.6

 

0.2

292.6

 

0.6

243.5

 

1.0

217.5

 

0.2

280.4

 

0.8

265.0

 

0.4

254.0

 

1.0

296.8

 

0.12

226.9

 

 

 

 

 

0.1s

 

1.0

 

 

361.6

 

0.1s means device with substrate.

 

CHAPTER 6
HYDRODYNAMIC TRANSPORT MODEL FOR THE MODFET

6.1 Hydrodynanric Transport Equations for the MODFET

The basic device topology of the MODFET (modulation-doped field effect
transistor) is shown in Figure 6.1, and the energy band diagram along the y direction is
shown in Figure 6.2. Electrons accumulate at the interface between the two materials
as a result of band bending. The ionized donor impurities in the AlGaAs and the con-
duction band electrons in the quantum well are spatially separated. This results in a
substantial reduction in the ionized impurity scattering in the GaAs layer which leads
to enhanced electron mobilities, particularly at low temperatures. The potential well
formed at the heterojunction is normally narrow enough to have quantized energy lev-
els in the y-direction and behaves as a two-dimensional electron gas. Up to this date
there are relatively few numerical MODFET solutions using the hydrodynamic equa-
tions. Two reports are those of Widiger et al. [15] and Shawki et al. [54].

In Widiger’s model, electron transport in the MODFET can take place in b0th the
bulk GaAs and in the quantum well (assuming the AlGaAs is depleted). In the region
near the source where the fields are small the conduction process may be attributed
mainly to the lowest quantum sub-band. In contrast, at the pinchoff region towards the
drain end of the gate, the electron current will be principally in the bulk. Electron
transfer between the quantized subbands 2-D gas and the bulk 3-D gas is present
whenever the average electron energies are comparable to the subband energy spac-
ings. Widiger assumed only the lowest subband was a two-dimensional gas system.
The higher subbands were all treated as part of the three dimensional gas system. This
may be assumed since, when conditions allow significant transfer out of the lowest

subband, the higher subbands will be spaced sufficiently close so as to approach

116

117

Source Gate Drain

 

 

 

 

_ _
I x
AlXGa1_xAs :n
Y
Alea1_xAs : u
GaAs : u

 

 

 

 

 

Figure 6.1 Basic device topology of the MODFET.

AleaLxAs : u
l

Aleal-xAs=n one.

Ionized quantum well E
donors ///7 c
j] EF

 

 

 

 

 

 

Figure 6.2 Energy band diagram for a MODFET.

118

three-dimensional properties. The quantum well and bulk systems are then coupled by
allowing the electrons to scatter between the two-dimensional and three-dimensional
gas systems. The equation for transport in the quantum well are similar to those of a
three-dimensional gas ((2.44)-(2.46)) with the addition of coupling terms since elec—
trons can transfer to the bulk system. These coupling terms are analogous to

generation-recombination terms. Widiger assumed no current conduction in the

AlGaAs layer.

Shawki et al. [54] applied a hydrodynamic energy model that is valid in the
framework of a gradual variation in Al alloy composition to simulate the entire MOD-
FET region including the top AlGaAs layer. This model treats the electrons as a.
three-dimensional electron gas without considering two-dimensional quantum well
effects. This is justified from Monte Carlo simulations of submicron MODFET’s[55]
which demonstrated that electrons during the major part of their travel under the gate
are not confined in quantum subbands if the drain bias is high enough to heat the elec-
trons (for VDs= 0.5 V, the total subband population is found to be less than 15 percent
of the equilibrium one). Also, the literature results for the two-dimensional electron gas
(2DEG) behavior obtained from classical models based on either Boltzmann or Fermi-
Dirac statistics are generally comparable with those obtained from the exact self-
consistent solution of Schrodinger’s and Poisson’s equations [56]. This model without

2DEG treatment makes the model more compatible with that of MESFET’s.

This study will extend the MESFET HTM model in previous chapters to include
a heterojunction for MODFET simulation. The electron transport dynamics based on
the particle, momentum, and energy conservation equations for a heterojunction are

described[54] by:

an __
BT+V(nV)-O (6.1)

119

.31an M + V.v_nm_'_<§>_v ”.1.
q

 

 

 

 

 

at q Il(§)
V nk T
= n VW + 1) - ( B (DY) + lVlnm‘ (§) (6.2)
q no (I

8 €‘%’% TOY

M + V-n v(§ + k3T(§)'y) + V-Q = -n v-V(-q \II-x) - ii (6.3)
at T§(§)
F 3,201)
-.- —— 6.4
7m) F l/2(Tl) ( )
(Efn - Ec)
= 6.4b
kBT(§) ( )
1 e 2 3
w = 3»: (av + -2-krT(§) (6.5)
§ = w + Up . ‘ (6.6)
The three moment equations are solved together with Poisson equation
View!) = --<I(ND+ - n) ' (6.7)
where
N
M; = (2,. _ 5., . (6.8)
1 + 2e keT(§)

It should be mentioned that the conservation equations (6.1)-(6.3) presented above are
different from those we used for the simulation of MESFET’s in previous chapters in
two aspects. The first difference is the terms describing the gradient in electron affinity

and the second difference is the dependence of effective mass on position.

The model described above can be simplified as

V-evty = -q (ND -n) (6.9)

120

an

 

3;- — V-J = 0 (6.10)
3—"54 + v-s = J~E—n H” (6.11)
at 1:;
where
J = ‘4’"
= 11(an + V(nkBT)) (6.12)
s = sci-(g + kBT)-IQ (6.13)
E=-V[\V+-§-+%£lnNc]. (6.14)

Notice that the term associated with lnNc came from V7 in the momentum and energy
conservation equations (6.1) and (6.2) which can be derived by following the work of
Azoff [57]. Q is the heat flow shown in equation (2.54). The equations (6.9)-(6.14)
are the same as those used in the work of Horio for the heterojunction bipolar transis-
tor[58].

The boundary conditions for the MODFET are the same as those for the MES-
FET except for the ohmic contacts and the heterojunction. Usually there is a highly
doped region under source and drain contact to reduce the source and drain resistance.
The current across the heterojunction in the highly doped ohmic contact region is
mostly due to tunneling conduction [59]. The highly doped ohmic contact regions
behave as a small resistor r to parallel with the large heterojunction resistor R as
shown schematically in Figure 6.3(a). This will reduce the current across the hetero-
junction in non-ohmic contact region. Hence, it is expected that most current flows as
shown in Figure 6.3(b). Very little current flows across the heterojunction in the non-
ohmic region. Also, the current flow in the top layer is very small due to depletion by

the gate. This is similar to Widiger’s work[15]. Different boundary conditions have

121

Source Gate Drain

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(a)
Source Gate Drain
_
+
n Aleal_xAS :n + d3
11
Alea] £3 :u
(b) \‘ ‘ -—
GaAs : u I
Sourcel I‘GS ! LG ! LGD Drain g
AleaHAs :n d 1
AJIGaIJAs: u d2
(C) <— 1 _—
GaAs : u
A

 

Figure 6.3. Equilibrium ohmic contact boundary condition model. (a) The ohmic con-
tact behaves as a small resistance r in parallel with the large heterojunction resistance
R. (b) Most cunent flows through r, not R. (c) A simplification of (b) based on equilib-
rium ohmic contact boundary condition at (e-f) and (g-h).

122

been used in the works of Loret [60] and Shawki [54], which results in wide spread
current distribution across the heterojuction, not just confined to the source and drain
ohmic contact regions. In this study, the ohmic contacts (e-f) and (g-h) in Figure
6.3(c) are used. The potential, electron concentration and electron energy are calculated
for this vertical boundary by assuming local thermal equilibrium conditions.
Specifically, the electric field and current can terminate at these boundaries and no
current flows in the vertical direction as shown in Figure 6.3(c). This is achieved by
solving the potential \[I and the electron density n subject to the energy condition
T = T 0- The current J and the energy flux S which cross the heterojunction are
assumed to be zero. The drain current is calculated by summing the current through

the cross section A in Figure 6.3(c).
6.2 Simulation Results and Discussion

The structure used for the simulations had LGS =0.6um , LG = 0.4um ,
LCD = 0.6% , d1=50nm, d2=5nm, d3=150nm, and d4=320nm as shown in Figure
6.3. The doped #3001qu has n-type doping ND =5x1017cm'3, the undoped
Alea1.,As has n-type doping ND=1x1014cm‘3 and the undoped GaAs has n-type

’3. For GaAs, the relative dielectric constant is 12.9 and the

doping ND=lx1014cm
density of states in the conduction band is N, =4.350x1017cm‘3. For Al, Ga 1_,As, the
fraction x in Al, Ga 1.,As is 0.3 with a correspond relative dielectric constant of 12.3,
and a density of states in conduction band of NC =6.988><1017(:m"3 . The conduction
band discontinuity is 0.23 eV. The Schottky barrier potential is 0.75 V [54,56] and the

low field mobility is assumed to be 0.5 sz'lS‘l.

Figure 6.4 shows the current-voltage characteristics for VG =0.0V and VG =0.2V.
Figure 6.5 shows the conduction band edge distribution. It is seen that most of the vol-
tage drop is under the gate edge near the drain. Figure 6.6 shows the electron distribu-

tion where electrons have accumulated at the heterojunction due to the conduction

 

123

 

.8 — V65 = 0.2V

ID(A/cm)

 

 

 

 

0 I I I

0 .2 .4 .6 .8 I
V05 (Volt )

—

Figure 6.4. The current-voltage characteristics of the MODFET for VGS = 0.0V and
V05 =0.2V.

124

 

 

 

 

 

 

 

EC (eV)

750 1
.161 -
- L+88
-1 08
and V03 = 1.0 V. The source is the front left region and the drain is the front right region.

The gate extends from 0.8 pm to 1.2 pm .

Figure 6.5. The conduction band distribution (eV) over the entire device for V03 = 0.0 V

125

 

 

 

 

 

 

 

n(m‘3)
5.00 -
/'\
m 3.33e
(\J
C)
X
v 1.67"
.000 -\
.569

 

 

 

 

Figure 6.6. Electron density distribution ( m‘3 ) over the entire device for V6; = 0.0 V
and Wm = 1.0 V. The source is the front left region and the drain is the front right region.

The gate extends fiom 0.8 pm to 1.2 pm.

126

band discontinuity. The density of accumulated electrons under the gate region is less
than the density under the source and drain contacts due to the gate depletion. The
electron energy is shown in Figure 6.7 in the GaAs layer. Note that the electron
energy is relaxed to the room temperature energy at the drain ohmic contact. The drain
ohmic contact has a depth of d3=150nm and the total simulated depth is d4=320nm .
Figure 6.8 shows the longitudinal current density in the GaAs layer where most of the

current flows close to the heterojunction.

The charge-partitioning method has used to calculate the current unity gain. For
Vos=1 V, VGs =0.0V, the CP method gives gm=134mS/mrn, cg =4.7x10’1°mF/mm and
fT=45.37GHz. For VDs=1 V, VGs=0-2 V, the CP method gives gm=266mSlmm,
cg =6.44x10-1°mF/mm, and fT=65.73GI-Iz.

For the S 3A simulation of the MODFET, one change has been made to the
equilibrium boundary conditions (e-f) and (g-h) in Figure 6.3(c). The temperature at
equilibrium is fixed at room temperature, T=To, so the small signal temperature is
T=O. The MODFET is biased at VDs=1 V and VGs=0-0 V. Figure 6.9 shows the y-
parameters at the frequencies 10, 20, 30, 40, 50, 100 and 150 GHz. Using the y-
parameters, the unilateral power gain GU and the current gain A I are calculated as
shown in Figure 6.10 using equations (5.19) and (5.21). Note that the unity current
gain fT using y-parameters gives 51.48 GHz and recall that the CP method gives
45.37 Ghz. This means that the CP method underestimates by 10 percent. The simu-
lator had slow convergence for higher frequencies making the determination of f max

difficult. However, the device performance can still be studied through the unity gain
frequency f T .

127

E.(eV)

  

I

/

 

//Il

 

.388 r

.2711

.155“

Figure 6.7. Electron energy distribution (eV) over the GaAs layer for V65 = 0.0 V and

1.0 V. The source is the front left region and the drain is the fiont right region.

The gate extends fi'om 0.8 pm to 1.2 urn .

Vos =

128

 

 

 

 

 

 

 

 

 

Figure 6.8. Longitudinal electron current density distribution (A/m 2) over the entire

 

0.0 V and V03 = 1.0 V. The source is the front left region and the drain

is the front right region. The gate extends from 0.8 pm to 1.2 pm .

device for VGs

129

 

4_ 1’11

 

1'22

 

1’12

-25 1’21

 

 

 

 

Realm

Figure 6.9. The y-parameters versus frequency for the MODFET using HTM. The seven
frequencies are 10, 20, 30, 40, 50, 100 and 150 GHz.

130

 

30—3
20—

- G
Gaul10— U

((113)
0 \
-\B

-10-

 

AI

 

 

 

I I MFTTITF I I I IIFII
101 102 103

Frequency (GHz)

Figure 6.10. Gy and A, versus frequency for the MODFET using HTM. The seven
frequencies are 10, 20, 30, 40, 50, 100 and 150 GHz.

CHAPTER 7
CONCLUSIONS AND RECOMMENDATIONS

The hydrodynamic transport model applied to semiconductor device simulation is
an important area of investigation. The HTM makes a large improvement over the
DDM with only moderate increases in computation time. Additionally, it takes much
less computation time than Monte Carlo simulation and provides results that agree well

with those from Monte Carlo simulation.

Chapter 2 of this dissertation derived the moment equations from the Boltzmann
transport equation based on a parabolic band structure. Instead of solving the
Boltzmann transport equation directly, this moment equations model solves for three
characteristic quantities to describe the carrier distribution function. This model was
then applied to modeling the short channel MESFET after reviewing modeling work
by other authors based on somewhat different approximations.

Up to this date, the simulators using the HTM were primarily developed using the
finite-difference method due to its simplicity and regularity. Also, the finite-difference
method allows the application of the Scharfetter-Gummel discretization in order to
handle the exponential dependence of electron concentration on potential and tempera-
ture. The finite-element method is difficult to apply to HTM because the electron den-
sity is exponentially dependent on voltage and temperature. The box-integration
method, which is essentially a finite-difference method based on a triangular mesh or
rectangular mesh, links both methods and is a good approach to discretize the sem-
iconductor equations since the Scharffeter-Gummel technique can also be used in the

discretization. It was the box-integration method which was used in this study.

The HTM simulator was applied to the DC solution of the MESFET in Chapter 4.

The material dependent parameters including electron velocity, electron energy and

131

132

energy relaxation time were calculated using a one-particle Monte Carlo simulator.
The accuracy and the validity of the HTM simulator was checked by simulating a
GaAs MESFET in two-dimensions using both the HTM simulator and a Monte Carlo
FET simulator. Since both simulations originated from a common set of material
parameters, any differences were due to the difl’erent models. The two simulator
showed good agreement. The HTM simulator was then used as part of a study to
investigate the bias dependence of the parasitic source and drain resistances in the
MESFET. It was shown that the high energy of the electrons entering the drain region
significantly changes the drain resistance at high drain voltage values.

The AC solution of FETs using the sinusoidal steady-state. analysis (S 3A) tech-
nique applied to the HTM was developed in Chapter 5. The previous AC solution
methods for the hydrodynamic transport equation were conducted by using either the
charge-partitioning method or Fourier decomposition method. The S 3A method has an
advantage in that the accuracy of the AC solution is easier to obtain and verify. This
application of the S 3A method to the hydrodynamic transport equations is one of the
key contributions of this work. The S 3A method was then used to simulate and study

the millimeter wave performance of submicrometer GaAs MESFETS.

The HTM model was then modified and applied to simulating the DC and AC
behavior of MODFETs. The MODFET model included the effects of the heterojunc-
tion in the solution of the Poisson equation and the transport equations. The MOD-
FET study showed that the current implementation of the HTM to heterojunctions has
some limitations for the S 3A AC solutions at high frequencies. Some areas for future
work include the flow of currents across the heterojunction in MODFETs and the solu-

tion of MODFETs at very high frequencies.

APPENDICES

APPENDIX A

APPENDIX A
FINITE-ELEMENT DISCRETIZATION FOR DRIFT-DIFFUSION MODEL

The finite-element method [24,63] produces approximations ‘I’j h and n 1- h to the

exact solution for potential Wj and concentration n j at point j . The approximations

may be conveniently formulated by defining a ’shape function’ 0,- so that the approxi-

mations to w and n becomes for m degrees of freedom

‘1'” = itheiO‘J) (A1)

i=1

and

.4 = inihegOCJ) (A2)

i=1
where rm" denotes \p"(x,-y,-) and ni" denotes n”(x,-y,-). The shape functions 0,-(x,y)

are defined in two-dimension such that

6,-(xjyj) =0 fori ¢j fori = 1,2,...m (A3)
and

0,-(xj'yj) =1 fori =j fori =1,2,..,m . (A4)

The approximations \th and hi” are found using the Galerkin method. The Galerkin
method applied in two-dimensions may be defined for each element in terms of the

residual R,- , (error in" the solution), in the i"' element as
n .
EIRieildA = 0 (A5)
i=1 '

where n is the number of elements. 0,1 is the j ‘h shape function in the i ‘h element.

133

134
Consider the generalized set of non-linear partial differential equations
V-Fi (u ,Vu )-c,- (u .12) = 0 i = 1,2,3,...,n (A6)

where F represents any physical flow quantity like electric flux density D or current
density J. u denotes the unknown variables such as potential \v or electron density n.

Applying the two-dimension approximation as shown in (A1) and (A2) gives

u(x.y)= Eon-9100’) (A7)

i=1

where

a,- = norm) . (A8)
The Galerkin condition requires that

Ri(or) = [00,(v-F-c)do = 0 . (A9)
Using the identity

V-GiF = 0,-V-F + F-VO, (A10)
we have

[new-Fan = [rem-rd: - jar-vows) (All)

where the divergence theorem is used to transform the area integral over (2 into a line

integral over I‘. Using (A11), (A9) can be rewritten as
R,(ot)= jr0,F-rd1-jn(ve,-F + 0,.ch i =1,2,3,...,m . (A12)

Hence the sum over the n finite elements is given by

R (or) = imam-fiat - Ina/err + email] = 0. (A13)
i=1

The integrals over F vanish because of the boundarys and the definition of the shape
functions [24] ( for example Vw-fi’ = ary/871’ is the derivative normal to the boundary

135

which is set to zero for fiee surfaces). Hence, (A13) reduces to

R (a) = i[In(V9,-'F + 0,-c)dQ] = 0. (A14)
i=1

This matrix expression is used to solve for a,- = u (151,)?” = 1,2,...,m .
The Poisson equation is expressed as
V-eVut + q(ND -n) = 0 . (A15)

Comparing (A15) with (A6), F = ve and c = -q(ND -n) are obtained. Using
(A1),(A2) and (A14), the Poisson equation (A15) can be discretized as

[K‘Vlw‘ - qlMlle" - n") = 0 (A16)
where

x,y-V = [nave-voids (A17)
and

MU = Ina-aids . (A18)

Similarly, the continuity equation for electrons is expressed as

V-J — egg-+40 = 0 (A19)

and
J =4(nllnE +Dnvn) (A20)

where G is the net generation rate. Comparing (A19) with (A6), F = J and
c = q(3n lat )-qG are obtained. Using (Al), (A2) and (A14), the continuity equation
for electrons (A19) can be discretized as

ah

It n h_ h=.
at +[K ]n B 0 (A21)

 

41M]

where

136

Ki," = [nerve(ejttnhra:h +DnhV9j)ds (A22)
and

B,- =Inq0,-G"ds . ’ (A23)

(A16) and (A21) are the discretized equations using finite element method for Poisson

and continuity equations, respectively.

APPENDIX B

APPENDIX B
DERIVATION OF THE CURRENT DENSITY AND ENERGY
FLUX EXPRESSIONS USING THE MODFIED
SCHARFE'I'I‘ER-GUMMEL TECHNIQUE

This appendix derives in detail the modified Scharfetter-Gummel expression for
current density and energy flux expressions used in this study.

(i) General Scharfetter-Gummel discretization :
Suppose we have the following general expression between nodes i and i+1 in
the x -direction,

A 1 dn
_ = — +_ ’ B1
wJ “1w" dx ( )

where A and al are assumed to be constant between nodes i and i+1, J is any physi-
cal flux quantity which is a function of x, w is the independent variable which is also
a function of x , n is the electron concentration which exponentially depends on vol-

tage and temperature.

(Bl) can be solved by using the solution for linear first-order differential equation

as follows.

Integrating factor = exp [0 1I%dx] = exp alIT:--21W—(‘f—‘:-)dx
(—)
dx

—1
= exp [a1 [%] lnw] = exp [azlnw] (B2)

where

-1
(12 :01 [ail] . (B3)

137

138

Multiplying the integrating factor (132) to (Bl), and integrating (B1) gives

XIV-1

x'+1
nexp(a21nw)| J; = A1 I -%exp(a21nw)dx (B4)
x,

where J has been assumed constant in the interval and hence is moved out of the

integration. The right hand side becomes

 

 

- i _ .1. .,__1_(d_w_
RHS -JAjWexp(a21nw)ax _JAjww (511) dz )dx
dx
= JA w"
dW 02
( dx)
ex ((1 lnw
= A p 2 ) . (135)
dw
“2‘71?
Substituting (B5) into (B4) gives
xi+1 1 Jim
nexp(a21nw)| x. =JATexp(a21nw)l x, (B6)
1 t
02(3)

"i+1CXP [0211mm] " nicxP [azmwi]

= JA-—1—{exp [a zlnwifi] - exp azlnwi ]} (B7)

Adz”)
- 1 } (BS)
(B9)

 

 

 

14" w.
ni+1exp [azln ‘31 ] - ni = JA—fi— {exp [azln :4
1 a2(_;) 1

Wi+1

 

 

az(-dlv-) "i+lcxP[021n
dx i

J =
A Wi-I-l

exp azln - l

 

i

139

This expression is simplified by defining

 

 

 

x -a In W’ -—a lnwm (1310)
- - 2
b 2 Wi+l Wt
which gives
a n- e -x -n-
J _____1_ t+1XP( b) t . (Bll)
A exp(—xb) - 1

(B11) can be simplified further by making use of the Bernoulli function according to

exp(-xb) _ 1 _ _1_ xb
c"IX-x1») - 1 1 " ¢XP(xb) xb 1 " exp(xb)

 

 

 

 

 

 

 

= 1 I), = --1—B (x ) (B12)
xb exp(xb) - 1 xb b
1 1 ‘16 1
= = B — . BI
exp(-xb) - l -x,, exp( -x,,) - l xb ( xb) ( 3)

Therefore,

a
J = 71' {—1—} [@413 (1b) - ":3 (-xb)]

= EA—l [‘1‘] ["IB(-xb) - "i+1B(xb)]

 

xb
= .1_(fl)_l_. [n-B(-x ). n- B(x )] (B14)
A dx 1n wt 1 b i+1 b °
Wi-l-l

The following conclusion can be made

it] =al-‘13-n '1' fl (B1)

then

140

 

 

 

 

1 dw 1
J = -A—(E)—w—,-- [n;B(-x,,) - nt (1%)] (B14)
1n
Wi+1
where
In W; a1 W‘- (B3)
x = a = _-
b 2 Wi+1 _d_w_ Wi-l-l
dx
8( ) 7’ (BIS)
x = .
b c"P000 - 1
(ii) Current density J discretion:
Recall (3.20) which is
i- : VT(ct + Din + Vn . (1316)
111' T

Compared with (B1), if A is replaced by 11:, w by T, and a1 by VT(0t -I- 1), then

(B 14) becomes

 

 

 

J," = 154%) 1T- [niB (~xb) - nHIB (xb)] (B17)
Ti+1
and (B3) becomes
x, =(ot+1)1n ‘ (BIS)
Ti-I-l

The mesh point m is midway between points i and i+1. Equations (B16),(B17) and
(B18) are the same as (3.20), (3.22) and (3.23), respectively.

(iii) Energy flux S discretion:
Recall (3.28) which is

—§— = VT(ot + A5‘1)%(Tn) + V(nT) . (1319)

-u8T

141

Compared with (31),rr1 is replaced by s, A by '—1,w by T, n by nT, and a1 by

 

 

 

118
VT (or + Ali-1), then (314) becomes
41,8) cg) 1,. [(nr).B(—x.) - (nT)r+tB (xt)] (320)
In t
Ti+1
and (B3) becomes
x, = (or + A6-1)ln ‘ . (321)
Ti-i-l

Equations (B19),(B20) and (B21) are the same as (3.28), (3.29) and (3.31), respec-

tively.

APPENDIX C

APPENDIX C
RIGRID ALGORITHM

The correct allocation of the mid is a crucial issue in device simulation. The grid
mesh has a direct influence on the simulation time and solution accuracy. In order to
maintain the simulation time within reasonable bounds and to have satisfied solution
accuracy, it is desirable to design a remid algorithm to allocate fine mid in some
regions and coarse mids in others. The mesh can be remided based on a triangular
element basis as shown in PISCES[28], or based on a rectangular basis as discussed in
this appendix.

For a given rectangular mesh with nxxny mid points, there are nx-l columns
and ny-l rows as shown in Figure C.1(a). There are two phases in this remid algo-
rithm. The first one is remid by columns. The second is remid by rows.

In the first phase, the columns are remided from column 1 to column nx -1.
Each column is checked as shown in Figure C.1(b). The checking may result in two
cases: (1) the mid remains unchanged, or (2) the mid is changed as shown in Figure
C.1(b). The mid is changed when a column has any two horizontal points across
which the chosen variable changes by more than a specified tolerance. The chosen
variable can be potential 0;, electron density n ,or electron energy i. The refinement
is done by adding a vertical line inside this column as shown in Figure C.1(b) and by
evaluating the variables on the new added line using interpolation. If a line is added,
then the total number of mid points is increased by ny. Note that the mid points are
renumbered during the process of remiding due to these added new lines. After the
first phase is done, the new nx will be larger than or equal to the initial nx depending

on the remid tolerance.

142

143

column column
0 o o “X’ 1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

row 1 +1
row 2 —>
(a) :
row ny-l —>
fly 2111’ nx 0 fly
(b) (1) or (2)
(c) 0
(1)
or
(2)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure C. l. Remid algorithm for rectangular mesh. (a) The mesh contains nx-l columns
and ny-l rows, (b) remid for each column, and (c) remid for each row.

144

In the second phase, the rows are remided from row 1 to row ny-l as done in
the first phase for the columns. Similarly to the first phase, the row remid may result
in two cases as shown in Figure C.1(c). If a horizontal line is added, the total number
of mid points is increased by nx.

For example, given an initial mid with nx xny , if the remid algorithm generates 2
more lines during the column remid and generates l more line during the row remid,
then the total number of mid points is increased from nxxny to (nx +2)x(ny +1). For
MESFET simulation, since nx is larger than ny the increase in ny is more expensive
than the increase in nx because of an increase in the bandwidth of the matrix which
simrificantly increases computation time. Different tolerances can be set in column

remid and row remid to keep ny as small as possible.

LIST OF REFERENCES

pd
O

10.

LIST OF REFERENCES

C. Moglestue, "Monte Carlo particle modeling of small semiconductor devices,"
Computer Methods in Applied Mechanics and Engineering, vol. 30, pp. 173-208,
1982. ‘

R. W. Hochney and J. W. Eastwood, Computer Simulation Using Particles,
Chapter 10, New York, McGraw-Hill, 1981.

Y. Awano, K. Tomizawa, N. Hashizume and M. Kawashima, "Monte Carlo parti-
cle simulation of a GaAs short—channel MESFET," Electronics Letters, vol. 19,
pp. 20—21, 1983.

A. Yoshii, M. Tomizawa and K. Yokoyama, "Accurate modeling for submicron-
gate Si and GaAs MESFET’s using two-dimensional particle simulation," IEEE
Trans. on Electron Devices, vol. ED-30, pp. 1376-1380, 1983.

B. Camez, A. Cappy, A. Kaszynski, E. Constant, and G. Salmer,"Modeling of a
submicrometer gate field-efi‘ect transistor including effects of nonstationary elec-
tron dynamics," J. Appl. Phys, vol. 51, pp. 784-790, 1980.

R. K. Cook and J. Frey,"Two—dimensional numerical simulation of energy tran-
sport effects in Si and GaAs MESFET’s," IEEE Trans. on Electron Devices, vol.
ED-29. pp. 970-977, 1982.

C. M. Snowden and D. Loret, "Two-dimensional hot-electron models for short-
gate-length GaAs MESFE ’s," IEEE Trans. Electron Device, vol. ED-34, pp.
212-223, 1987.

Y. K. Feng and A. Hintz, "Simulation of submicrometer GaAs MESFET’s using a
full dynamic transport model," IEEE Trans. Electron Device, vol. ED-35, pp.
1419-1431, 1988.

W. R. Curtice and Y. H. Yun, "A temperature model for the GaAs MESFE ,"
IEEE Trans. Electron Device, vol. ED-28, pp. 954-962, 1981.

O. L. El-Sayed, S. El-Ghazaly, G. Salmer and M. Lefebvre, "Performance
analysis of sub-micron gate GaAs MESFETS," Solid-State Electronics, vol. 30,
pp. 643-654, 1987.

145

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

146

F. A. Buot and J. Frey, "Effects of velocity overshoot on performance of GaAs
devices, with design information," Solid-State Electronics, vol. 26, pp. 617-632,
1983.

R. Stratton, "Diffusion of hot and cold elecu'on in semiconductor barriers," Phys.
Rev., vol. 126, pp. 2002-2014, 1962.

H. H. On and T. W. Tang, "Numerical modeling of hot carriers in submicrometer
silicon BJT’s," IEEE Trans. on Electron Devices, vol. 34, pp. 1533-1539, 1987.

R. K. Cook, "Numerical simulation of hot-carrier transport in silicon bipolar
transistors," IEEE Trans. on Electron Devices, vol. 30, pp. 1103-1110, 1983.

D. Widiger, I. C. Kizilyalli, K. Hess, and J. J. Coleman, "Two-dimensional tran-
sient simulation of an idealized high electron mobility transistor," IEEE Trans.
Electron Devices, vol. ED-32, pp. 1092-1102, 1985.

K. Blotekjaer, "Transport equations for electrons in two-valley semiconductors,"
IEEE Trans. Electron Device, vol. ED-l7, pp. 38-47, 1970.

T. J. M. Boyd and J. J. Sandbom, Plasma Dynamics, Barnes & Noble, New
York, 1969.

T. J. Maloney, "Non-equilibrium electron transport in compound semiconductors,"
PhD dissertation, Cornell Univ., Ithaca, NY, 1977.

S. Kratzer, "Computer simulations of electron transport in GaAs," M.S. thesis,
Cornell Univ., Ithaca, NY, 1978.

P. A. Sandbom, A. Rao and P. A. Blakey, "An Assessment of approximate nons-
tationary charge transport models used for GaAs device modeling," IEEE Trans.
Electron Devices, vol. ED-36, pp. 1244-1253, 1989.

M. R. Friscourt, P. A. Rolland, A. Cappy, E. Constant, and G. Salmer, "Theoreti-
cal contribution to the design of millimeter-wave TEO’s," IEEE Trans. on Elec-
tron Devices, vol. 30, pp. 223-229, 1983.

T. W. Tang, "Extension of the Scharfetter-Gummel Algorithm to the energy bal-
ance equation," IEEE Trans. Electron Devices, vol. ED—31, pp. 1912-1914, 1984.

23.

24.

26.

27.

28.

29.

30.

31.

32.

33.

34.

147

D. L. Scharfetter and H. K. Gummel, 'Targe-simlal analysis of a silicon read
diode oscillator," IEEE Trans. Electron Devices, vol. ED-16, pp. 64—77, 1969.

J. J. Barnes and R. J. Lomax, "Finite-element methods in semiconductor device
simulation," IEEE Trans. Electron Devices, vol. ED-24, pp. 1082-1089, 1977 .

R. S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1962.

S. Selberherr, Analysis and Simulation of Semiconductor Devices, Springer-Verlag
Wien New York, 1984.

E. M. Buturla, P. E. Cottrell, B. M. Grossman and K. A. Salsburg, "Finite-
element analysis of semiconductor devices: the FIELDAY program," IBM J. Res.
Develop. vol. 25, pp. 218-231, 1981.

M. R. Pinto, C. S. Rafferty and R. W. Dutton, PISCES-II User’ s Manual, Stan-
ford Electronics Laboratories, 1984.

A. Gnudi, F. Odeh, P. Ciampolini, M. Rudan, and G. Baccarani, " On the discret-
ization of the hydrodynamic model of the semiconductor equations," Workshop on
Numerical Modeling of Processes and Devices for Integrated Circuits--NUPAD
11," May 9-10, 1988.

H. K. Gummel, "Self-consistent iterative scheme for one-dimensional steady state
transistor calculations," IEEE Trans. Electron Devices, vol. ED-l 1, pp. 455-465,
1964.

A. Yoshii, M. Tomizawa and K. Yokoyama, "Investigation of numerical algo-
rithms in semiconductor device simulation," Solid-State Electronics, vol. 30, pp.
813-820, 1987.

C. M. Snowden, M. J. Howes, and D. V. Morgan, "Large-signal modeling of
GaAs MESFET operation," IEEE Trans. Electron Devices, vol. ED-30, pp. 1817-
1824. 1983.

S. M. Sze, Physics 0f Semiconductor Devices, ch. 6, 2nd ed. New York: Wiley,
1981.

P. L. Hower and N. G. Bechtel, "Current saturation and small-signal characteris-
tics of GaAs FET’s." IEEE Trans. Electron Devices, vol. ED-20, pp. 213-220,
1973.

35.

36.

37.

38.

39.

40.

41.

42.

43.

45.

148

H. Fukui, "Determination of the basic device parameters of a GaAs MESFET,"
Bell Syst. Tech. 1., pp. 711-797, 1979.

P. Urien and D. Delagebeaudeuf, "New method for determining the series resis-
tance in a MESFET or TEGFET," Electron. Lett., vol. 19, pp. 702-703, 1983.

K. Lee, M. Shur, K. W. Lee, T. Vu, P. Roberts, and M. Helix, "New interpreta-
tion of ’end’ resistance measurements," IEEE Electron Device Lett., vol. EDL-5,
pp. 5-7, 1984.

S. Chaudhuri and M. Das, "On the determination of source and drain series resis-
tances of MESFET’s" IEEE Electron Device Lett., vol. EDL-5, pp. 244—246,
1984.

K. W. Lee, K. Lee, M. S. Shur, T. T. Vu, P. C. T. Roberts, M. J. Helix, "Source,
drain, and gate series resistances and electron saturation velocity in ion-implanted
GaAs FET’s," IEEE Trans. Electron Devices, vol. ED-32, pp. 987-992, 1985.

L. Yang and S. 1. Long, "New method to measure the source and drain resistance
of the GaAs MESFET," IEEE Electron Device Lett., vol. EDL-7, pp. 75-77, 1986.

Y. H. Byun, M. S. Shur, A. Peczalski and F. L. Schuermeyer, "Gate-voltage
dependence of source and drain series resistances and effective gate length in
GaAs MESFET’s," IEEE Trans. on Electron Devices, vol. 35, pp. 1241-1245,
1988.

D. H. Navon, "A theorem on heat generation in semiconductor devices," Proc. of
the IEEE, vol. 66, pp. 520-522, 1978.

M. S. Adler, "Accurate calculations of the forward drop and power dissipation in
thyristors," IEEE Trans. on Electron Devices, vol. ED-25, pp. 16-22, 1978.

Y. Awano, K. Tomizawa and H. Hashizume, "Principles of operation of short-
channel gallium arsenide field-effect transistor determined by Monte Carlo
method," IEEE Trans. on Electron Devices, vol. ED-31, pp. 448-452, 1984.

W. Fawcett, A. D. Boardman and S. Swain, "Monte Carlo determination of elec-

tron transport properties in gallium arsenide," J. Phys. and Chem. of Solids, vol.
31. PP. 1963-1990, 1970.

47.

48.

49.

50.

51.

52.

53.

54.

55.

56.

57.

 

149

K. Brennan and K. Hess, "High field tranSport in GaAs, InP and InAs," Solid-
State Electronics, vol. 27, pp. 347-357, 1984.

T. A. Grotjohn, "Compensation effects on the electron velocity in submicron
GaAs MESFET’s," IEEE Trans. on Electron Devices, vol. 35, pp. 1144—1145,
1988.

S. E. Laux, "Techniques for small-signal analysis of semiconductor devices,"
IEEE Trans. Electron Devices, vol. ED-32, pp. 2028-2037, 1985.

G. D. Vendelin, Design of amplifiers and oscillators by the S-parameter method,
pp. 32-33, Wiley, 1982.

M. B. Das and P. Schmidt, "I-Iigh-fi'equency limitations of abrupt-junction
FET’s," IEEE Trans. Electron Devices, vol. ED-20, pp. 779-792, 1973.

M. B. Steer and R. J. Trew, "High-frequency limits of millimeter-wave transis-
tor," IEEE Electron Dev. Lett., vol. EDL-7, pp. 640-642, 1986.

R. J. Trew and . B. Steer, "Millimeter-wave performance of state-of-the-art MES-
FET, MODFET and PBT transistors," Electronics Letters, vol. 23, pp. 149-151,
1987.

H. O. Vickes, "Note on unilateral power gain as applied to submicrometre transis-
tors," Electronics Letters, vol. 24, pp. 1503-1505, 1988.

T. Shawki, G. Salmer, and O. El-sayed, "MODFET 2-D hydrodynamic energy
modeling: optimization of subquarter-micro-gate structures," IEEE Trans. Electron
Devices, ED-37, pp. 21-30, 1990.

U. Ravaioli and D. K. Ferry, "MODFET ensemble Monte-Carlo model including
the quasi-two—dimensional electron gas," IEEE Trans. on Electron Devices, vol.
33, PP. 154-156, 1986.

J. Yoshida, "Classical versus quantum mechanical calculation of the electron dis-
tribution at the n-AlGaAs/GaAs heterointerface," IEEE Trans. Electron Devices,
ED-33. PP. 154-156, 1986.

E. M. Azoff, "Generalized energy-momentum conservation equations in the relax-
ation time approximation," Solid-State Electronics, vol. 30, pp. 913-917, 1987.

58.

59.

61.

62.

150

K. Horio, Y. Iwatsu, and H. Yanai, "Numerical Simulation of AlGaAs/GaAs
Heterojunction Bipolar Transistors with Various Collector Parameters," IEEE

Trans. on Electron Devices, vol. 36, pp. 617-624, 1989.

S. J. Lee and C. R. Crowell, "Parasitic source and drain resistance in high-

electron-mobility transistors," Solid-State Electronics, vol. 28, pp. 659-668, 1985.

D. Loret, "Two-dimensional numerical model for the high electron mobility
transistor," Solid-State Electronics, vol. 30, pp. 1197-1203, 1987.

R. W. Hockney, R. A. Warriner, and M. Reiser, "Two-dimensional particle
models in semiconductor-device analysis," Electronics Lett., vol. 10, pp. 483-486,
1974.

C. Moglestue and S. J. Beard, "A particle model simulation of field effect transis-
tors," ed. B. T. Browne and J. J. H. Miller, Numerical Analysis of semiconductor

Devices, Boole Press, Dublin, Ireland, pp. 232-236, 1979.