1018 IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. ED-30,NO.

9, SEPTEMBER 1983

Semiconductor Device Simulation

WOLFGANG FICHTNER,MEMBER, IEEE, DONALD J. ROSE, AND RANDOLPH E. BANK

Absfracf-The most effective way to design VLSI device structures is. to in [l]are applicable. See other papers in this issue [2] and
use sophisticated,complextwo-dimensional (2D)and three-dimensional [3] for a summary of the current state of process simula-
(3D) models. This paperand its companion 111 discusses thenumerical
simulation of such device models. Here we describe the basic semiconcluc- tion.
tor equations including severalchoices of variables. Our examples illustrate The coupled nonlinear partial differential equations de-
results obtained from finite-difference and finite-element implementations. scribing the intrinsic behavior of semiconductor devices
We stress thenecessary 3D calculationsfor small-size MOSFET's. provide a significant challenge to the scientific computing
Numerical results on inter-electrode capacitive coupling are included.
community. While most of the work in this field has been
I. INTRODUCTION performed by researchers with electrical engineering and
physics background, there has recently been increased in-
T HE development ofnew semiconductor technologies
and novel semiconductor device structures has [Link]-
tionally been guided by anexperimental approach. Starting
terest by numerical analysts in investigating this problem.
As a result advanced numerical algorithms such as sparse-
matrix techniques and adaptive multigrid methods have
from an established process sequence, fabrication steps are been applied to solve the semiconductor equations.
changed together withgeometrical (feature size) dimen- The guiding principles in the design of our software
sions. The modified process is then realized by fabricating package were robustness, speed, and easyuseraccess.
several lots. Finished devices are tested to insure their Robustness is an absolute necessity for any software crea-
performance conforms to the initial design. This approtach tor, otherwise hiscodeswill notbe used at all. In our
usually includes several iterations of the processing to computing environmentwith a large number of active users
testing loop. from different areas, robustness has priority. Nevertheless,
With the advent of increasingly complex integrated cir- execution times were of major concern to us. The search
cuits, the traditional empirical approach has becomeex- for a best device under worst case operating conditions can
pensive and time consuming. An alternative approach using result in tens of different runs, where one run can be a
sophisticated numerical simulations in process and device complete device I-V curve. The interface to a program is
development has proved to beboth costeffective and crucial for obtaining easy user access. We have attempted
reliable. to relieve the user completely from any aspects of informa-
For example, the development of a new CMOS process tion transfer between the process and the device simulator,
might involve nine lithography steps, six ion implantations grid generation and numerical fine-tuning.
and several diffusion, annealing, and oxidation steps. In a In Section 11, we present the semiconductor equations, a
medium-size computer, one cansimulate all critical process set of coupled nonlinear partial differential equations
steps and the corresponding device performance in a matter (PDE's). After we formulate the equations, we discuss
of minutes to hours. A real experiment, on the other head, boundary conditions. The description of high-field phe-
would usually take from several weeks to months. nomena and heavy doping effects is also included in this
For devices and circuits ofVLSI complexity,prozess section.
conditions are tightly coupled to the behavior of the finished Section I11 deals mainly with the simulation of intrinsic
device. Therefore, device simulation cannot be a stand alone silicon
device structures. We state our approach to
field, but has to be closely coupled to process simulation. hierarchical simulation as a costeffective albeit physical
In this paper and its companion [l] to follow we pre,lent way to analyze devices. Several examples are included in
our approach to semiconductor device simulation. We as- this section, illustrating results obtainedfrom finite-dif-
sume that the necessary input from process simulation is ference and finite-element implementations. Furthermore,
available. We do not treat process simulation extensively we presentmethodology and results from three-dimen-
here, although many of the numerical techniques presented sional (3D) simulations. We stress the necessity of 3D
calculations for small-size MOSFET's. Part of this section
Manuscript received June 2, 1983. isalsodevoted to user-oriented device simulation. One
W. Fichtner and D.J. Rose are with Bell Laboratories, Murray Hil, NJ example illustrates how we have coupled a general purpose
07974. two-dimensionalprocess simulator to our device simulation
R. E. Bank is with the University of California at San Diego, La Jolla,
CA 92403, and is partially supported underONR Contract N00014-132-K- software.
0197. Parasitic elements begin to have an increased impact on
0018-9383/83/0900-1018$01.00 01983 IEEE
Copyright held jointly: 01983 IEEE; 015183 Society for Industrial and Applied Mathematics
FICHTNER et al.: SEMICONDUCTOR DEVICE SIMULATION 1019

device and circuit performance. In Section IV, we shall and momentum coordinates to one in space coordinates
present numerical results on inter-electrode capacitive cou- alone. In this approximation, it is assumed that the re-
pling. sponse of carriers to a change in the electric field is much
faster than the effective rate of change in the field; i.e., the
11. SEMICONDUCTOR EQUATIONS device response and the dielectric relaxation processes in
the device interior are slow compared to the fundamental
In this section we formulate the partial differential equa- response of the carriers. As a consequence, the product of
tions (PDE’s) which describe the static and dynamic behav- the maximum particle velocities and relaxation times is
ior of carriers in semiconductors under the influence of negligible compared to active device dimensions, especially
external fields. We state the equations in a general form for the case of silicon.
valid for all semiconductors of practical importance, al- Assuming the above conditions hold, we can write the
though our primary interest is in silicon devices. However, basic equations of semiconductor transport in the form
the following discussion is not restricted to silicon; by using most commonly used in numerical device simulations
theappropriate material constants the equations canbe [18]-[20].
easily transformed and generalized. We state various The static and dynamic behavior of carriers under the
boundary conditions which are encountered in real device influence of external fields can be partitioned into three
simulations. Several changes of variable which simplify the groups: Maxwell’s equation, the current-density relations
PDE system are discussed. Finally, we briefly discuss how and the carrier continuity equations.
hot-carrier effects and high-doping phenomena can be a) The Maxwell equations form the basis of all elec-
incorporated. tromagnetic phenomena in semiconductors. They are given
A . EquationFormation by
dB
The equations which describe the transport of electrons VXE=--
at
and holes in semiconductor devices can be derived in a
straightforward way from the Boltzmann Transport Equa-
tion (BTE), if the motion of carriers is treated tobe
semiclassical [4]. In this simplified, but nevertheless physi-
cal description, carrier motion in a semiconductor crystal
with applied external field can be regarded as a series of where the electric displacement vector D and the electric
acceleration (treated by classical mechanics) and scattering field vector E are related by the constitutive relation
(treated by quantum mechanics) events.,Forthe electric
field strengths and temperatures encountered in small sili- D =€,E
con devices, this semiclassical picture accounts verywell with the semiconductor permittivity E , . A similar relation
for all carrier transport effects of interest (e.g., [4]-[6]). holds between the magnetic induction vector B and the
The applicability of the BTE for the description of small magnetic field vector H
semiconductor structures hasbeen recently studied by vari- B =p H
ous authors [7]-[9]. Based on the results of these papers, we
can conclude that for silicon devices with active dimensions with the permeability p . In (l),Jcond and Jtot are the total
of 0.1 pm and larger, carrier transport can be described by current and the conduction current, respectively; p is the
the semiclassical BTE approach. The method of solution of total charge density.
the BTE is guided by the active dopant concentration level, When combined with the current equations, (1) provide
which is normally well [Link] 10l6 cmP3 for small devices a complete description of carrier dynamics in semiconduc-
to [Link] current-voltage behavior. At these doping tor structures.
levels, however, it has been shown that carrier-carrier in- For our purpose, the Poisson equation is of major inter-
teractions dominate the shape of the carrier distribution est. It relates the total space-charge to the divergence of the
function [lo]-[12]. This fact allows the assumption of an electric field, namely
effective carrier temperature. Recently, Hess [13] has criti-
cally reviewed the applicability of the carrier temperature 6 , v . E = - esv2qb
=p (2)
concept for the simulation of silicon device structures. His where the electrostatic potential qb is defined via
results indicatethat for MOS devices and CCD’s, the
carrier temperature concept is well suited to describe high- E=-VJ,
field transport effects. and the total space-charge p, assuming complete ionization,
Based on these assumptions in the semiclassical picture, is given by
the operation of devices can be evaluated by solving the p=-q(n-p+N) (3)
BTE [14]-[16], although for real applications the direct
solution of the BTE is impractical. However, we can go one where N = N; - N i is the electrically active net impurity
step further and make the quasi-static local field approxi- concentration, q is the electric charge, and n and p are the
mation [17], which reduces the problem from one in space electron and hole carrier densities.
 1020 NO. ED-30, VOL.
DEVICES,
ELECTRON
IEEE TRANSACTIONS
ON 9, SEPTEMBER 1983

b) The current equations for electrons and holes are given c) The continuity equations for electrons and holes are
by given by
J , = qp,nE + qDnvn (4)
J p = W p P E- (5)
where p n and p p are the electron and hole mobility, and D,, 1
- V . J p - G + R + - =aP
O
and Dp are the corresponding diffusion coefficients. 130th 4 at
mobilities and diffusion coefficients depend on the electric where G incorporates generation phenomena, such as im-
field. pact ionization or carrier generation by external radiation
Undernondegenerate conditions, the diffusion coeffi- and R describes recombination processes. We shall com-
cients and the mobilities are related by the Einstein [Link]
ment on specific expressions for these terms in Subsection
PI, ~211,c221 11-D.
Equations (l),(15), (16) with either (4) and ( 5 ) or (13)
and (14) describe the current flow in the semiconductor
and determine the electrical performance of the devices.
where k is Boltzmann’s constant and T is the [Link] In this paper we are concerned with the numerical
temperature. solution of the steady-state semiconductor equations. Thus
For nondegenerate materials, if the electrostatic poten- we summarize (1)-(6) as
tial is taken with reference to the intrinsic Fermi level E;,
the carrier densities are given by the Boltzmann app:roxi-
~ ~ v 7=. qE( p - n N ) + (17)
mation

Equations (17)-(20) are usuallynormalized intodimen-

where n i is the intrinsic density and E , = q+F is the Fermi sionless form. We have followed the work of de Mari [23];
level position under equilibrium conditions. From (7) and several other ways exist and have been reported e.g., [24]. If
(8), it follows immediately, that we express all spatial dimensions in terms of the intrinsic
Debye length,
n p = n 2i . (9)
Undernonequilibrium conditions, the electron and hole
concentrations will departfrom their equilibrium values
and they can no longer be represented via the single
all densities in units of n; and all voltage terms in units of
quantity E,. However, we may introduce two new energy
k T / q , (17)-(20) read
parameters E,,, and E,, , which allow us to write n and ,D in
a form similar to that in (7) and (8), namely -Au+n-p-k=O (21)
-0.j =v.j =g-r
P (22)
jn = - p n [n v u - o n ] = - p,nVu (23)
jp=-pp[p~u+~p]=-ppp~w
(24)
where E,, = q+,, and E,, = q+p. The quantities E,,, and with = eU-0 p = ew-Ua
E,, are the quasi-Fermi levels (also called imrefs), and c&, (25)
and are the corresponding quasi-Fermi potentials. Here u, u , and w are the normalized potential and the
The simple relationship in (9) is replaced by electron and hole quasi-Fermi level, respectively, and k is
the normalized impurity distribution. In the transition from
(7)-(20) to (21)-(25), we have kept the symbols n , p , and
p n and p p . The reader should be aware that in (21)-(25),
these are actually normalized variables.
and the difference in the quasi-Fermi levels indicates the
The PDE system (21)-(25) is posed on a region (in 2D)
deviation of the np-product from its equilibrium value.
such as shown in Fig. 1. For MOStechnology the gate
Rewriting (4) and ( 9 , using (6), (lo), and (ll),yields for
contact is insulated from the channel between the source
the current densities in terms of quasi-Fermi levels
and drain contacts in the siliconregionbySiO,. In this
4 = - qPnnV+n (13) oxide region, (21) is augmented by
Jp = - 9PpPv+p. (14) - <,AU= 0 . (21a)
FICHTNER et al.: SEMICONDUCTOR DEVICE SIMULATION 1021

to the early work of Gummel [25] on one-dimensional

bipolar transistor modeling. This algorithm has since been
extensively used both in its original one-dimensional form
and also in two dimensions,by Slotboom [26], [27], Mock
[28], [29], Heimeier [30], Manck [31], and others. The major
difference between the various attempts to solve the semi-
conductor equations lies in the choice of variables and the
treatment of the carrier continuity equations [32]. In our
companion paper, we present a discussion of the various
possibilities that exist.
The basis for Gummel's algorithm [25]is as follows.
Fig. 1. Idealized cross section of a planar n-channel MOSFET. Estimates are made for the quasi-Fermi level distributions
u and w, and a solution for u is obtained from Poisson's
Equation (24)is normalized form, and c 0 is the ratio of equation. The solution for u so found is then used in .the
permittivities in silicondioxide and silicon. transport equations to improve the initial estimates of w
and v. The cycle is repeated until some convergence crite-
B. BoundaryConditions rion is satisfied. The important point is that the equations
Equations (21)-(25) summarize the coupled system of are decoupled. We shall call this solution process Block
PDE's describing the semiconductor device. It remains to Nonlinear Iteration (BNI). Although Mock [33]has claimed
specify boundary conditions for a particular geometry. that under certain circumstances the algorithm does not
Fig. 1 shows a representative example a two-dimensional converge, the scheme itself has proved to be of consider-
cross section of an n-channel MOS device. For the ohmic able value.
contacts of the source, drain, and bulk we derive Dirichlet Some workers have chosen the carrier concentration n
boundaryconditions as follows. Assuming space-charge and p as variables in which to apply the Gummel algo-
neutrality on these contacts ( p = 0 in (23)) we obtain rithm, for example de Mari [23], Petersen [34], Manck [31].
This choice of variables is attractive as it can be seen that if
n--p-k=0. (26)
an estimate for u is known, (22) and (23) are linear in n and
Carrier equilibrium defines the condition p , and that given an estimate for n and p , then (21)is
pn = 1. (27) linear in u. Consequently, only successive solutions of
linear problems are required. However, because of the large
Recalling (25) we obt,ain u = w = u F , where u F is the range of values to be accommodated, care must be taken in
normalized equilibrium ]Fermilevel(see (7) and (8)). Fi- the discretization scheme,as small errors may produce
nally u F is taken to be equal to the applied voltage at the inadmissible (and unphysical) negative carrier concentra-
contact; e g , for the drain tions. In addition, large spatial variations make an ade-
u=w=uDs (28) quate numerical representation of derivatives difficult.
where uDs is the normalized drain-to-source voltage. Others [26], [27], [35], [36], [37] used exponentials of the
From (26) and (27) we find quasi-Fermi levels as dependent variables

n=k/2+/(k2/4)+1 jp=-p P P~ w = - pe-' V e " = - pPe-'vw (30)

which combined with (25) and (28) produces the Dirichdet j,, = - p,,nvv = p,e've-' = p,,e"vv. (31)
condition at the drain
The BNI algorithm produces linear transport equations in
u=uDs+ln(k/2+\i'o+1). (29)
the new variables v = e-" and o = ew,although the Poisson
The other contacts are treated similarly. equation is nonlinear. An advantage of this scheme is that
Along the side boundaries, homogeneous Neumann (re- the quantities v and w are nearly constant over large
flecting) boundary conditions insure that no current can regions of the device, thus allowing a reduction of mesh-
flow in orout of the device along these edges. This points. The range of numerical values, however, is greater
statement is equivalent to equality of drift and diffusion than that in n and p . In our companion paper [l],we
along the sides. suggest several strategies to scale the original equations
The gate contact is treated as a Dirichlet boundary (30) and (31).
condition with an appropriate work function for the Mock[28], Sutherland [38], and others [39],[40]have
material. Along the side of the gate oxide homogeneous adopted the so-called stream function as a variable in the
Newmann boundary conditions are used. continuity equations. The main assumption is thatno
generation and recombination process exist, i.e.,
C. Choice of Variables and Earlier Work
The origin of semiconductor device modeling using
v .j , = 0 (32)
numerical instead of analytical methods can be traced back v .j p = 0. (33)
1022 IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. ED-30,NO. 9, SEPTEMBER 1983

This divergence-free flow of the electron and hole cur- simultaneous solution of (21) to (23)[45]. By ordering
rent can be enforced by introducing a stream-function 8 as the physical variables u, u , and w in blocks, we can utilize
a vector potential for j , e.g., for the electrons the advantages of botha global Newton strategy and
efficient linear equation solution (for details, see [l]).
j, = j,v x e, 134)
where & is the total current. For a two-dimensional prob- D. High -Field Phenomena and Heavy Doping Effects
lem the components of j , are, respectively In (22) and (23), the terms g and r describe the variable
generation-recombination phenomena which occur in sili-
con devices. The relative importance of recombination
phenomena depends on the particular device; e.g., for
bipolar devices, recombination strongly influences the cur-
rent gain, but for MOS devices, it is usually not of impor-
Identifying the x and y components of (24)with the tance. On the other hand, it is imperative to include the
corresponding expressions of (35) and (36) we have avalanche generation term under highfield conditions.
Furthermore, under strong breakdown conditions, recom-
A

J,-
ae,a (ne-")
= p,e'- bination phenomena become important due to the large
aY ax (37)
increase of electron and hole densities. For complex struc-
ae,
4-=- a tures, it is not always possible to choose the "right" physics
ax p,e'- aY ( n e - ' ) a priori for the simulation, since the intrinsic behavior of
the structure might change for different operating points.
using a / a x ( n e - " ) = - e - u ( a u / a x ) and (24). In Section III-B, we present as an illustrative example the
By dividing both sides of these intermediate equations case of parasitic structure whichis typical for a CMOS
by p,,e', and differentiating the first equation partly with bulk technology.
respect to y , the second partially with respect to x, and We can partition the various phenomena according to
adding the two, we obtain
g - r = (g-r)thermal+(g)avf(g-Y)Auger t40)
where ( g-r)thermal describes thermal processes in the
volume and at the surface (Shockley-Read-Hall processes),
A similar equation holds for the holes. The iterative proce- the second term describes generation of carriers dueto
dure starts with initial guesses for u , e,, and [Link] the impact ionization and ( g-r)Augerstands for Auger phe-
nonlinear Poisson equation is solved, followed by a solu- nomena [46].
tion of (39) and the equivalent hole equation. The stre:am- For silicon, the dominant recombination mechanism is
function approach has recently been generalized to include an indirect process involving a trap center in the energy
generation-recombination processes [41]. gap. This process can be well modeled by the expression
One alternative to the above mentioned schemesis to
work in the quasi-Fermi levels themselves, as defined in (9)
and (10). This is a natural choice because, in addition to a
range compression, the current equations simplify to mme where r,, and rp are the electron-hole lifetimes and n ,and p,
extent and also the boundary conditions are directly ex- are associated with defect levels in the energy gap. The
pressed in u and w and are therefore simple. The arising most effective trap center occurs in the middle of the gap,
difficulty is that the continuity equations are nonlinear in u and in this case n , = p , = ni.
and w when an estimate for u is available, so that with this Surface recombination effects are modeled similar to
choice of variables all equations are nonlinear. However, (40),where the lifetimes are approximated by reciprocal
with a proper nonlinear equation solver this is a small price surface velocities.
to pay. If a simultaneous solution scheme is used, this The term gav describes the number of generated elec-
choice of variables seems preferable. tron-hole pairs per unit volume and time
A problem usually encountered with the decoupled solu-
tion scheme is slow convergence of the overall solution
under certain physical conditions such as boundary condi-
tions in high-current [Link], this implies a where a , and ap are the ionization coefficients and j,, and j p
strong coupling between the Poisson and the continuity are the current densities in (24) and (25).
equations. The coefficients a , and ap depend on the electric field.
Historically, Manck [42] recognizedthis problem first For multidimensional simulations, we treat the field depen-
and tried to avoid it bysolving the equations simulta- dence of a by the relations
neously using a LSOR-scheme to solve the linearized
system. Adler [43], and Buturla et al. [44] presented s8Au-
tions of the coupled system using sparse matrix techniques.
Recently, we have presented an alternative approach to the
FICHTNER et a / . : SEMICONDUCTOR DEVICE SIMULATION 1023

thus taking only the field component parallel to the current conditions. In the case of MOS devices, these operating
flow into account. conditions range from the subthreshold region to high
Auger recombination is the reverse process toimpact current conditions in the saturation region, possibly includ-
ionization, involving three particles. The recombination of ing impact ionization effects. If the device operates in the
an electron-hole pair [Link] energy for a third particle. subthreshold region, it is essentially close to turn-off and a
The probability for this process is proportionalto the very small current is flowing between source and drain (see
density of involved particles, i.e., Fig. 1). This current flowis caused by diffusion only,
resulting in a small perturbation in the quasi-Fermi level
(g-r)*uger=(n:-pn)(Cnn+~p~) (44) distributions, i.e.,
with the Auger coefficients C, and Cp. It is obvious from v .j n = 0 (48)
(44) that Auger recombination will only become important
if the electron and the hole density are high in the same v . j p = 0. (49)
area. Therefore, constant quasi-Fermi levels can be assumed for
In various silicon devices, the active device regions con-low-current operating conditions. This fact has been utilized
tain doping levels above 10" ~ m - This ~ . is particularlysuccessfully by various authors [21], [49].
true for bipolar devices, where heavily doped emitter-base The assumption of aconstant minority-carrier quasi-
regions have a crucial influence on the transistor behavior. Fermi level breaks down under operating conditions result-
For essentially all deuice,s, in which high doping levels areing in medium to high current flow. For our hypothetic
necessary, the transport of minority carriers is essential. n-channel example, this means that (48) no longer holds.
At these high doping levels, several phenomena occur However, if generation phenomena such as impact ioniza-
which are unique. These phenomena can be classified in 1) tion effects can be neglected, the majority-carrier current in
effects concerning the bandstructure of silicon, 2) an the substrate is zero, i.e.
increase in recombination with a corresponding decrease in
0 . j so
minority-carrier lifetime, and 3) questions concerning the P (50)
validity of Maxwell-Boltzmann statistics. for a n-channel device and the hole quasi.-Fermi level w is
In this paper, we shall not dwell on the physics of high constant. This reduces the problem to one of solving only
doping phenomena a complicated and still active area of two equations, namely (21) and (22) with (23). This case of
research (for a review, see [47]). constant quasi-Fermi level for majority carriers is the
To incorporate these high doping effects into our system standard case for most MOS simulations.
(21) to (25), we have to take into account Auger recombi- Under bias conditions resulting in majority-current flow,
nation, (44), and correct the intrinsic carrier concentration,
all three equations have to be solved. For MOS-like struc-
n,, which is replaced by an effective intrinsic concentration
tures, these situations usually arise under high bias condi-
tions with impact ionization.
The concept of modeling hierarchies can be generalized
for arbitrary semiconductor structures. Dirks and Engl[50]
where AE, describes bandgap narrowing. The effective have successfully used their stepsolving strategy to analyze
carrier concentration has to be used in (10) and (11) bipolar structures.
The utilization of the hierarchical solution scheme to a
complex simulation problem can result in dramatic reduc-
tion of CPU-time needed. The bottom-up philosophy allows
to start with the lowestlevel solution ( e g , a Poisson
(47) solution), which is a good starting approximation to the
next level, and so on. Fig. 2 presents a schematic summary
Experimental values for n,, are available in the literature of the hierarchy concept. The motivation for the use of the
[481. hierarchical solution schemeiseasy to understand when
the alternative for obtaining the same information-a full
111. INTRINSIC DEVICESIMULATION
solution of the coupled set of equations-is considered.
In this section, we present some practical aspects of In the next subsection, we present two simulation exam-
numerical device simulation. Several examples should il- ples where we utilized the hierarchy scheme.
lustrate our points. These examples have been constructed
from typical simulations;of MOS (metal-oxide-semiconduc- B. Intrinsic Device Simulation
tor) devices. In Subsection III-A, we discuss hierarchical We have developed a software package for the solution
device simulation, and Subsections III-Band III-C give of (21)-(25) which is capable of dealing with most two-
two-dimensional and three-dimensional examples. dimensional device structures arising in modern IC tech-
nology. The package consists of two parts, a finite-dif-
A. Hierarchical Device Simulation ference (FD) code and a finite-element (FE) code. In our
In simulations one is normally interested in the behavior environment, with many users having widely varying prob-
of semiconductor devices over a wide range of operating lems, this has proved to be a good solution to the tradeoff
1024 IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. ED-30,NO. 9, SEPTEMBER 1983

DISTANCEALONG GATE (pm)

0
0

Fig. 2. Hierarchical device simulation,

0.4

DRAIN (IV)
-E 0.2
I

X
L
0.3

04

0.5

(b)
Fig. 4. (a) Surface and (b) contour plot of the electron density.

equations. This allows us to take some advantage of vector

processing capabilities available in our CRAY-1.
The finite-element code is used for device structures with
nonrectangular geometry. Our implementation is based on
Bank's PLTMG [51] package. We use elements of triangu-
lar shape with linear basis functions. This package utilizes
a multigrid strategy with adaptive mesh refinement. For
details on the discretization, we refer the reader again to
our companion paper (see [l, sec. 3.51). The implementa-
tion of our FE package is inherently more costly in the
assembly process than the FD code.
In the following, we compare both implementations on
an example of a planar MOSFET structure with a 1-pm

,L 0 5 (pm)
DISTANCE

(b)
gate length, 0.25-pm junctions and a gate oxide thickness
of 250 A. This device has a constant substrate doping of
lo1' cm-3 and has not been implanted in the channel. We
chose bias voltages of V g s = V,, = 0 V and v d , = 1V. Under
these conditions, the device operates close to punchthrough
Fig. 3. (a) Surface plot of the electrostatic potential, (b) contour plol as mode with part of the total current flowing through the
seen from the bulk for Vgs= Vhs= 0 and Vds= 1V. bulk. The total current is I d = 0.7 mA ( W = 30 pm). The
device operation is illustrated in Figs. 3-6. Fig. 3 shows a
surface plot (a) and the corresponding contour plot (b) of
between general software and computer resource [Link]- the electrostatic potential. We note that the depletion edge
ments. ( - 0.3-V line) has been pushed deep into the substrate. Fig.
The finite difference code is used for device structures 4 presents the surface plot (a) and contour plot (b) of the
comprised of unions of rectangles. We use a tensor product electron concentration (note the different depth scales).
mesh, which is the same for all three equations. The b83x Electric field plots and a gate view of the lateral current
method (see [l],Section 3.1)isused to discretize the density are shown in Figs. 5 and 6.
FICHTNER et ul.: SEMICONDUCTOR DEVICE SIMULATION 102s

TABLE I
function
Plug -
lIXll
!valuations
IlUV
-
0 0.11253~+04 [Link]@OC+00 0
I 0.51084eCO3 0.38984e+00 2
2 0.60836~+02 0 .I 8 ~ + 0 0 I
3 0.20109e+02 0.28712c-01 I
0 4 0.65792~+01 0.8618oS-01 I
5 0.16669~+01 0.49291c-02 I
6 0.19886c+W 0.22268c-02 I
7 0.37518s-02 0.34863e-03 1
0 0.65112e+06 [Link]+w 1
I 0.6405oS+06 O.I5065S+OI 4
2 0.50558~+06 0.56559~-01 3
(a) 3 0.15807~+06 0.13187c-01 1
4 0.47666~+05 0.77781~-02 1
LATERAL DISTANCE ( p m ) I 5 0.14562~+05 0.58847~-02 I
I 00 6 0.50547c+04 0.51 113e-02 I
7 0.17475~+04 0.43824~-02 I
8 0.53914e+03 0.34663e-02 1
9 0.1199Oe+03 0.21785e-02 1
10 O.l0935c+O2 0.77098~-03 I

11
0 0.33154e+03
1 0.14364e+03
2 2 0.26625~+02
3 0.31638c+Ol
4 O.l339Oc+W
0 0.89318~+01 0.00WOefW I
~ 0.32373s+Ol 0.61739~-02 I
0.21459e+w
2 0.1003oS-02 I
3 0.10522e-02 0.63353~-04
-
1
0 O.l9269e+Ol [Link]+w 1
O.l8654e+W 0.34266~42 I
SUBSTRATE
2 0.2WJoS-02 0.29986~-03 I
0 0.68586e+W o.W000c+w
0.50 0.36458~-01 0.20029c-02
0.10573c-03
2 0.95605~44
0 0.2935&+(11 0.00000e+00
Fig. 5. (a) Surfaceplot of thelateralelectric field, and (b) thecorre- O.llM4c-02
sponding contour plot of the lateral electric field (kV/cm) near the 2 0 . 6 1 4 2 5 ~ 4 5 0.211ooe-04
drain contact.
Fermi levels (iteration = 0). This is equivalent to the lowest
complexity in our solution hierarchy, and it gives a good
initial condition for the electrostatic potential. Then, the
electron continuity equation .is solved, alternating with the
Poisson equation, corresponding to next :level of complex-
ity in the solution hierarchy. The need €or a stable algo-
rithm to treat the nonlinear equations becomes obvious
from the result in Table I. In order for the equations to
converge rapidly, it is imperative to select the right step-
size. The algorithm used allows a proper choice by
reevaluating the right hand sides,whichis a relatively
cheap process.
Fig. 6 . Surface plot of the electron current density, as seen through the The use of the hierarchical solution methodology can
gate.
result in a large reduction of CPU time. In this example,
eleven BNI iterations (“plugs”) would have been necessary
Although the current is relatively high under these bias without the initial Poisson solution. For operating points
conditions, no impact ionization occurs due to the modest resulting in high current or in cases requiring a solution of
drain bias (V,, = 1 V). Therefore, we have only solved the all three equations (such as avalanche breakdown)the BNI
Poisson equation and the electron continuity equation. We iteration strategy would be replaced by aL (coupled) New-
used a block nonlinear iteration scheme (BNI, see earlier) ton-like iteration (see [l]and [45]).
since our FE code presently does not allow a fully coupled, The result obtained from the FE code is similar to the
simultaneous scheme. For both cases, our choice of vari- FD case. An initial triangulation with NV= 847 vertices
ables results in a nonlinear Poisson equation, and a linear and N T - 1 5 3 2 triangles was used, which should be com-
Continuity equation (see [l eq. (3.38)]).The finite difference pared with a total number of 4070 gridpoints of the F D
run was performedonone gridlevelwhereas the finite run.
element run used to multigrid levels with uniform triangle The power of our FE package lies in its’ability to model
refinement. complex nonrectangular structures. To illustrate this point,
In Table I, we illustrate several performance aspects of we have chosen a structure which arises in twin-tub CMOS
the FD, calculation. The inonlinear iteration scheme for the technology [54]. Fig. 7 schematically shows a cross section.
Poisson equation is the one we proposed in [53]. Initially, Depending on the kind of dopant, this structure can con-
wesolve the Poisson eq,uation assumingconstant quasi- tain several parasitic devices, and it can operate in various
1026 IEEE TRANSACTIONS ON
DEVICES,
ELECTKON VOL. ED-30,NO. 9, SEPTEMBER 1983

In the following, we present the basic methodology and

several results of three-dimensional numerical solutions of
the nonlinear system

r
n+
- EPI
- SUBSTRATE
Fig. 7. Schematiccrosssection of theparasiticdevices
CMOS technology.
which arise in
[ I
d u ,u,w)
h(z)= g2(u,u,w) =O
d u ,U,W)

corresponding to the coupled PDE system in (21)-(25).

Fig. 11 shows the basic geometry of a planar MOSFET
(49)

in three dimensions. We have employed a“standard”

seven-point finite difference scheme on tensor-product non-
uniform grids to discretize the linearized equations. To
account for the coupling along directions of current flow,
the unknowns are arranged into blocks for different planes
in the z-direction. The number p of planes varies between
15 and 30, depending on the operating conditions. To
illustrate the interdependence of the variables, Fig. 12
shows the coupling graph G(h’), h’ being the Jacobian of
(49) for p = 4. Each node represents approximately 2K
unknowns, 10K nonzeros and has the structure of a con-
ventional 2D grid graph. Since the CRAY-1 has no virtual
memory capabilities, we were forced to utilize out-of-core
storage in the calculations and to recalculate matrix ele-
ments. In analogy to the 2D case, the linear systems of the
,
,
r FIELD O X I D E form
Hx=-h (50)

have been solved by block iteration. The sub-blocks corre-

(4 sponding to the variables u , u, and w have been solved by
Fig. 8. (a)Parasiticdevices:p-channel MOSFET, (b)lateralbipolar an iterative method (vectorized SOR with RB-ordering).
p-n-p transistor,(c) vertical n-p-n transistor, and(d) SCR-structurc:. The results of 3D simulations have been successfully
used to design small, submicrometer-size devices fabricated
modes for different bias conditions. In Fig. 8(a) and 8(b), using a 1-pm NMOS process. As an example, we consider
the parasitic element is either a p-channel field transistor the enhancement-driver and the depletion-load device in
(a) oralateralp-n-p device (b). Situations c and d are this technology. Depletion devices are more susceptible to
equivalent to a SCR-structure and a vertical n-p-n device, geometry effects than enhancement devices, because of the
respectively. internal doping distribution. Fig. 13 shows a surface plot of
The automatically generated initial grid is shown in Fig. the two-dimensional net impurity concentration profile in
9, together with the doping profile for the case in Fig. 8(a) this device. In order to shift the threshold voltage to a
and (b). Fig. 10 presents surface plots of the electrostatic proper negative value, the channel of this deviceis im-
potential for the cases in Fig. 8(a) and 8(b). planted with a shallow, low-dose arsenic implant. The
maximum arsenic concentration in the channel is 3.1017
C. Numerical Simulation of Three - Dimensional Structwes ~ a junction depth of 800 A. Depending on the
~ r n -with
Silicondevices are inherently three-dimensional (3D) operating conditions, the device works either in bulk mode
structures. For most problems, however, the behavior, of ( Vgs< V,,) or in surface mode (V,, > V,,). Fig. 14 com-
devices can be modeled in either one or two dimensions. pares 2D and 3D results with experimental data for en-
This is normally the case for bipolar devices with a vertical hancement and depletion devices with a gatelength of 1p m
emitter-base-collector structure or for MOS devices with a and different widths. In the case of the driver device we
“wide” short-channel or a “small” long-channel geometry. can observe a moderate short-channel effect (illustrated by
For MOS deviceswith aspect ratios close to one, 3D the increase in Id at higher Vd$)and a proper dependence
simulation becomes necessary, especially in the sub- of I , on the channel width. In the case of the load device,
threshold regime and under high-field conditions. however, we can observe a large short-channel effect (the
The simulation of three-dimensional device structures current increases by two orders of magnitude for high vd,).
presents a formidable problem requiring large compul:ing The width dependence of Id is again linear, if the width is 2
resources. Nevertheless, severalsuccessful attempts have pm or larger. For smaller devices, however, the current is
been reported in the engineering literature [55]-[57]. reduced drastically. As we can see from agreement between
FICHTNER et d.: SEMICONDUCTORDEVICE SIMULATION 1027

10’9,
1 0’8 - - loi9
I 0’ 7/ - 10’6

10‘6

1015

I 0’4

(b)
Fig. 9. (a)Initialtriangulationand (b) doping concentration

the simulated and experimental results, 2D simulation is and an aluminum level at the top (conductors 3 and 4).
certainly adequate for devices with W > 2 pm. Smaller Since the situation is symmetric, we shown only half of the
devices, however, can only be modeled by a full 3D calcu- geometry.
lation. Similar conclusions have been obtained by the The capacitances associated with this system of conduc-
authors in the simulation of avalanche effects in depletion tors are given by the relation
devices [56].
Qi=Cj,(v,-v,)+Ci,(v,--~)+
Iv. SIMULATIOJV OF CIRCUIT INTERCONNECTS
+ CiiV;
+ * * + Ci,(V, - V , ) (51)
The ability to calculate accurately the capacitances of
interconnection wires is of major importance as the circuit where Q j , y , and Cii are the charge, thle potential and the
density increases and wire length and width decrease. For self-capacitance of the ith conductor, respectively; Cij is
very large circuits, the total delay on an integrated circuit the coupling capacitance between the ith and the j t h
chip isgiven to a large extent by the capacitance of the conductor.
interconnection wires. Generally, for a N-conductor problem, we obtainan
In the past, capacitances of conductors and wires have N X N capacitance matrix C which is ;simply defined as a
been calculated using the parallel plateapproximation, collection of all capacitances, analogous to (51). Histori-
thus neglecting effects such as fringing and coupling. As cally, the diagonal elements of the capacitance are called
transistors are scaled down, the width and the length of the coefficients of capacitance while the off-diagonal terms are
conductor are also reduced according to the scaling factor. called coefficients of electrostatic induction.
However, the thickness of the conductor and that of the If Maxwell’s equations hold, then (l),
dielectric are not scaled down, but reduced by different
aB
factors. This fact has two serious consequences. It not only vXEf-=O
increases the two-dimensional fringe effect but results in a at (52)
considerable increase of coupling capacitance where lines can be rewritten in terms of a vector potential A using
are neighboring. B=vXA (53)
A typical situation for today’s integrated circuits is shown
in Fig. 15. In this example, two levels of interconnects are to yield
used: a polysilicon level, made up of conductors 1 and 2, vx(vxA)=O (54)
1028 IEEE TRANSACTIONS O N ELECTRONDEVICES, VOL. E D - 3 0 , NO. 9, SEPTEMBER 1983

VIEW A I 2 3 4

DRAIN SOURCE n n n n

P+
1 p - T U B ,;Q-TuB
-______---*

n-SUBSTRATE
1 Pi

Fig. 12. Couplinggraphbetweenvariableplanes

Fig. 13. Surfaceplot of netimpuritydistribution of adepletionload

MOSFET.

B A G and since v X ( - v # )= 0, we obtain

(b)
Fig. 10. Surface plot of electrostatic potential distribution for parasitic
devices in Fig. 8(a) and(b).
In terms of the charge and current density, we can calcu-
late the electric field at time t and at r

where p' and J' are the charge and current densities in the
volume element dv' at time [ t - ( ( r- r'l/v)] with the prop-
agation velocity v. The exact calculation of E ( r , t ) would
involve the knowledge of the retarded source distributions.
In integrated circuit problems, however, we can ignore
retardation effects because the sources of interest are only
microns away from each other, leading to propagation
"ES
times much shorter than the time scales of circuit opera-
Fig. 11. Basic geometry of a planar MOSFET in three dimensions, tion.
 FICHTNER et ui.: SEMICONDUCTOR DEVICE SIMULATION 1029

-a
-0
-5
-4
-6

-7
I L =r p m
vqs=o'6v

01
- -8

-9

-IO

dep1e;idn device, 2D and 3D are the same for W=15,5 , and 2pm.

I e.,

I . - {+ , J+ CIZ
fC f I
Fig. 15. Schematiccrosssectionillustratingadoublelevelinterconnect
geometry.
Fig. 17. Capacitancebetweenconductorsas ar function of conductor
width/insulation thickness ( W / H ) for two thickness-to-height ( T / H )
Estimating the magnitude of the two terms for E leads to ratios.
the result that the second term may be neglected. There-
fore, we can calculate E from
Different techniques have been applied to analyze the
E=-v$. (57)
interconnect problem. For simple conductor arrangements,
It is therefore possible to use a solution of Laplace's analytical expressions [58] and numerical techniques [59]
equation have been successfully utilized. Systems with conductors of
nonrectangularshapeare difficult to simulate with the
-€VIc,=O ( 5 8 ) above techniques. We have used a multigrid finite element
to calculate the electric field. The total charge on a conduc-method to solve (58) for arbitrary geometries. The multi-
tor is then obtained by integrating the normal field compo- grid method lends itself naturally to problems of this kind.
nent along the conductor boundary. Fig. 16 shows the finite element triangulation at increased
1030 IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. E D - ~ O ,NO. 9, SEPTEMBER 1983

levels of refinementaround center conductor 1. Fig. 16(a) NJ: vanNostrand, 1950.

is the initial triangulation is automatica~~y generated. w. vanRoosbroeck, Be0 $>st. Tech. J . , vel. 29,
P.
560,
1950.
[20] R. C. Prim, 111, Bell Syst. Tech. J . , vol.30, p. 1174,1951.
Fig. 16(b) and (c) show the adaptively refined meshes. [21]A. H. Marshak, Solid-State Electron., vol. 21, p. 429,1978.
A typical result of a capacitance calculation isshow11 in [221 R. A. Smith, Semiconductors, 2nded. New York:Cambridge
University Press, 1978, ch. 7.2.
Fig* l7 for a three-conductor are [23]A. deMari, Solid-State Electron., vol.11, p. 1021, 1968;see also
defined in the inset of the figure. The self-capacitance C,, Solid-State Electron., vol. 11, p. 33, 1968.
the capacitance to ground C,, and the coupling capaci- [24] F. de la Moneda, IEEE Trans. Circuit Theory, [Link]-20, p. 666,
1973.
‘12 have been by biasing the center ‘On- [25] H. K. Gummel, IEEE Trans. Electron Devices, vol. ED-11, p. 455,
ductor at 1 V and by grounding all the others. The figure 1964.
shows the basic problem that arises if we reduce conductor [26] J. W. Slotboom, Electron. Lett., vol. 5 , p. 677,1969.
W. Slotboom, IEEE Trans. Electron. Devices, vol. ED-20, p. 669,
sizes. ~~l~~ a certain W / H ratio, the coupling capacitance [27] J.1973.
c,, increases above the capacitance-to-ground CII. ( h e [28]M. S. Mock, Solid-State Electron., vol. 16, p. 601, 1973.
possible way to make C, smaller is to reduce the conduc- [291 M. S. Mock, J . Eng. Math., vol. 7, p. 193, 1973.
tors thickness (see T / H = 0.5 in the figure), at the expense [30] H. H. Heimeier, IEEE Trans. Electron Devices, vol. ED-20, p. 708,
1973.
of an increased conductor resistance. [31] 0. [Link], IEEE Trans. Electron Devices, vol.
ED-22, p. 339, 1975.
v. CONCLUSION [32] D. ScharfetterH. and
V O ~ .ED-16, p. 64, 1969.
K. Gummel, IEEE Trans. Electron Devices,
The development ofnew device technologies has tradi- [33]M. s. Mock, Solid-State Electron., vol. 15,p.1, 1972.
tionally been guided by an experimental approach, The llSe [341 0. G. Peterson, R. A. Rikoski, and w . w. C o u h Solid-state
Electron., vol. 16, p. 239, 1973.
of software tools to aid process and device engineers [Link]- [351 p, Dubo&, Electron, Lett,, vel. 6, p, 53, 1970,
sents a challenging alternative to the original approach. We [36] E. J. Zaluska, Electron. Lett., vol. 9, p. 599,1973.
feel that the simulation of semiconductor structures by [371 c. R. Jesshope, E.D. Zaluska, and H. [Link], Electron.
accurate numerical models is an efficient way to [Link] [381 Lett., vol. 11, p. 285, 1975.
A, D, Sutherland, DARPA Rep, ECOM-75-1344-F-supp,
process conditions and the corresponding device perfor- [39] T. Toyabe, K. Yamaguchi, S. Asai, and M. S. Mock, in IEDM Tech.
mance. Dlg., paper 2006, Dec. 1977.
1401 E. Takeda. H. Kume, T. Tovabe, and S. Asai. IEEE Trans. Electron
REFERENCES
Devices, vol. ED-24, p. 611,-1982.
D. N. Pattanayak, R. A. Williams, and J. G. Poksheva, Appl. Phys.
R. E. Bank, D. J. Rose, and W. Fichtner, “Numerical methods for Lett., vol. 41, p. 459, 1982.
semiconductor device simulator,” this issue, pp. 1031-1041, 1983. 0. Manck, H. Heimeier, and W. End, IEEE Trans. Electron De-
See papers under Process Simulation, same issue. vices, vol. ED-21, p. 403, 1974.
W. Fichtner, Process Simulation in VLSI Technology, S. M. Sze, I:%. M. S . Adler, in Numerical Analysis of Semiconductor Devices, B. T.
New York: McGraw-Hill, 1983. Browne and J. J. H. Miller,[Link],Ireland: Boole Press,
[Link], High Field Transport in Semiconductors. NewYork: 1979.
Academic Press, 1967; see also Handhook on Semiconductors, vol I, E. M. Buturla, P. E. Cottrell, B. M. Grossman, and K. A. Salsburg,
[Link], E d . NewYork: North Holland, 1982. IBM J . Res. Develop., vol. 25, p. 218, 1981.
K. Seeger, Semiconductor Physics, [Link]:Springer, [Link] D. J. Rose,in Elliptic Problem Solves, M. H.
1982. Schultz,[Link]:AcademicPress, 1981.
G. Bauer, in Springer Tracts in Modern Physics, vol. 74, G. Hebcr, J, S. Blakemore, Semiconductor Statistics. NewYork:Pergamon
Ed. NewYork:Springer, 1974. Press, 1962.
D. K. Ferry, in VLSI Electronics- Microstructure Science, vol. I, R. [Link],R. J. vonOverstraeten,andH. J. deMan, in
N. G. Einsprech, Ed. NewYork:AcademicPress, [Link] also Advances in Electronics and Electron Physics, vol. 55, p. 77, 1981.
D. K. Ferry, I$. Hess, and P. Vogl, in VLSZ Electronics- Mic~,o- S . [Link], G. R. Srinivasan, and I. Antipov, in IEDM Tech. Dig.,
structure Science, vol. 11. paper 11.2, 1980.
K. Hess, in Advances in Electronics and Electron Physics, vol. (9, J. A. Greenfield and R. W. Dutton, IEEE Truns. Electron Devices,
Peter W. Hawkes, Ed. New York: Academic Press, 1982. V O ~ ED-27,
. p. 1520, 1980.
K. K. Thornber, J . Appl. Phys., vol. 51, p. 2127, 1980. H. Dirks and W. Engl,in Numerical AnalysisforSemiconducfor
C. J. Hearn, Proc. Phys. Soc., vol. 86, 881, 1965; see also Proc. Ph)ls. Devices, B. [Link] andJ. J. [Link], Eds. Dublin,Ireland:
SOC.,vol. 88, p. 407, 1966; Phys. Status Solidi ( h ) , vol. 64, p. 527, Boole Press, 1981.
1974 and Physics of Nonlinear Transport in Semiconductors, D. IC. R. E. Bank and A . H.
Sherman, in Sparse
Matrix
Ferry,J. R. Barkerand C. Jacoboni,[Link]:Plenum Proceedings-Z978, I. S. Duff andG. W. Stewart,[Link]-
Press, 1980. delphia, PA: SIAM, 1979.
R. Stratton, Proc. Roy. Soc., vol. A242, p. 355, 1957; see also Prcc. L. A. Hageman and D. M. Young, Applied Iterative Methods. New
Roy. Soc., vol. A246, p. 406, 1958. York: Academic Press, 1981.
H. Frohlich, Proc. [Link]., vol. A188,p. 521, 1947; see also If. R. E. Bank and D. J. Rose, Numer. Muth., vol. 37, p. 279, 1981.
Frohlich and B. V. Paranjape, Proc. Phys. Soc. (London), vol. B69, L. C. Parillo,L. K. Wang,R. D. Swenumson,R. L. Field, R. C.
p. 21, 1956. Melin, and R. A. Levy, in IEDM Tech. Dig., paper 29.3, 1982.
K. Hess, in Physics of Nonlinear Transport in Semiconductors, D. EL. E. M. Buturla, P. E. Cottrell, B. M. Grossman, B. M. Lawlor, C. T.
Ferry, J.R. Barkerand C. Jacoboni,Eds. New York:P1enu.n McMullen, and K. A. Salsburg, in IEEE Solid-State Circuits Conf.
Press, 1980. Dig. Tech. Papers, 1980, pp. 76.
P. A. Lebwohl and P. J. Price, Solid-State Communications, vol. 9,p. W. Fichtner, R. L. Johnston, and D. J. Rose. in Proc. 1981 Device
1221, 1971. Research Conf., Santa Barbara, CA., 1981.
Lebwohl. A. Yoshii, H. Kitazawa, M. Tomizawa, S. Hoiguchi, and [Link],
R. [Link] J. [Link], ComRuter Simulution Using IEEE Truns. Electron Devices, vol. ED-29, p. 184, 1982.
Particles. NewYork:McGraw-Hill, 1981, ch. 10. W. H. Chang, IEEE Trans. Microwave Theory Tech., vol. MTT-24,
H. D. Rees, J. [Link]. Chem. Solids, vol. 30, p. 643, 1969. 1976.
W. Shockley, Electrons and Holes in Semiconductors. Princeton, A. Ruehli, IBM J . Res. Develop., vol. 23, p. 626, 1979.
-