# 1 Boltzmann equation - Quantum Transport Theory Group

#### Document technical information

Format pdf
Size 231.3 kB
First found May 22, 2018

#### Document content analysis

Category Also themed
Language
English
Type
not defined
Concepts
no text concepts found

#### Transcript

```ECE539 - Advanced Theory of Semiconductors and Semiconductor Devices
Numerical Methods and Simulation / Umberto Ravaioli
Review of Conventional Semiconductor Device Models
Based on Partial Differential Equations
1
Boltzmann equation
The standard semi–classical transport theory is based on the Boltzmann equation
X
∂f
eE
+ v · ∇r f +
· ∇k f =
∂t
h̄
k0
[S(k0 , k)f (r, k0 , t)[1 − f (r, k, t)]
−[S(k, k0 )f (r, k, t)[1 − f (r, k0 , t)]
(1)
where r is the position, k is the momentum, f (k, t) is the distribution function (for instance, a
Fermi–Dirac in equilibrium), v is the group velocity, E is the electric field, S(k, k0 ) is the transition
probability between the momentum states k and k0 , and [1 − f (k0 , t)] is the probability of non–
occupation for a momentum state k0 . The summation on the right hand side is the collison term,
which accounts for all the scattering events. The terms on the left hand side indicate, respectively,
the dependence of the distribution function on time, space (explicitly related to velocity), and
momentum (explicitly related to electric field).
The Boltzmann equation is valid under assumptions of semi–classical transport: effective mass
approximation (which incorporates the quantum effects due to periodicity of the crystal); Born approximation for the collisions, in the limit of small perturbation for the electron–phonon interaction
and instantaneous collisions; no memory effects, i.e. no dependence on initial condition terms. The
phonons are usually treated as in equilibrium, although the condition of non–equilibrium phonons
may be included through an additional equation.
Analytical solutions of the Boltzmann equation are possible only under very restrictive assumptions. Direct numerical methods for device simulation have been limited by the complexity of the
equation, which in the complete 3–D time–dependent form requires seven independent variables
for time, space and momentum. In recent times, more powerful computational platforms have
spurred a renewed interest in numerical solutions based on the spheroidal harmonics expansion of
the distribution function. To–date, most semiconductor applications have been based on stochastic
solution methods (Monte Carlo), which involve the simulation of particle trajectories rather than
the direct solution of partial differential equations.
The vast majority of device simulations are normally based on the numerical solution of approximate models which are related to the Boltzmann equation, coupled to Poisson’s equation for
self–consistency.
In the simplest approach, the collision term on the right hand side of (1) is substituted with a
phenomenological term
feq − f (r, k, t)
(2)
τ
where feq indicates the (local) equilibrium distribution function, and τ is a microscopic relaxation
time. It is very useful to express the distribution function in terms of velocity, rather than momentum, since it will be easier to calculate electrical currents. In equilibrium we may use the
Maxwell–Boltzmann distribution function
1
feq (r, v) = n(r)
2πkB To
m∗
−3/2
m∗ |v|2
exp −
2kB To
!
(3)
where n(r) is the carrier density, To is the lattice temperature, kB is the Boltzmann constant, and
m∗ is the effective mass. The use of (3) for semiconductors is justified in equilibrium as long as
degeneracy is not present. The carrier density n(r) is directly related to the distribution function
according to
Z
n(r) =
dvf (r, v)
(4)
which is of general applicability. The significance of the momentum relaxation time can be understood if the electric field is switched off instantaneously and a space–independent distribution is
considered. The resulting Boltzmann equation is then
∂f
feq − f
=
(5)
∂t
τ
which shows that the relaxation time is a characteristic decay constant for the return to the equilibrium state.
2
Drift–diffusion model
The popular drift-diffusion current equations can be easily derived directly from the Boltzmann
equation. Let’s consider a steady state situation and for simplicity a 1–D geometry. With the use
of a relaxation time approximation as in (2) the Boltzmann equation becomes
eE ∂f
∂f
feq − f (v, x)
+v
=
(6)
∗
m ∂v
∂x
τ
Here, we have used the relation m∗ v = h̄k, which is valid for a parabolic energy band. Note that
the charge e has to be taken with the proper sign of the particle (positive for holes and negative
for electrons). A general definition of current density is given by
Z
J(x) = e
vf (v, x)dv
(7)
where the integral on the right hand side represents the first moment of the distribution function.
This definition of current can be related to (6) after multiplying both sides by v and integrating
over v. From the left hand side we get
1
J(x)
vfeq dv − vf (v, x)dv = −
(8)
τ
eτ
The equilibrium distribution function is symmetric in v, and its integral is zero. Therefore, we have
from (6)
Z
eτ
J(x) = −e ∗ E
m
Integrating by parts we have
Z
Z
Z
∂f
d
v
dv − eτ
∂v
dx
∂f
v
dv = [vf (v, x)]∞
−∞ −
∂v
|
{z
}
0
2
Z
Z
v 2 f (v, x)dv
(9)
f (v, x)dv = −n(x)
(10)
and we can write
Z
D
v 2 f (v, x)dv = n(x) v 2
E
(11)
where v 2 is the average of the square of the velocity defined as
D
v
2
E
1
=
n
Z
v 2 f (v, x)dv
(12)
For a purely 1–D treatment, the − 23 exponent in (3) may be replaced with − 12 , while the appropriate
thermal kinetic energy becomes kB2T instead of 3k2B T .
2
eτ
The drift–diffusion equations are derived introducing the mobility µ = m
∗ and replacing v
BT
with its average equilibrium value km
∗ , therefore neglecting thermal effects. The diffusion coefficient
D = µkBe To (Einstein’s relation) is also introduced, and the resulting drift–diffusion current is
dn
dx
dp
Jp = qp(x)µp E(x) − qDp
dx
Jn = qn(x)µn E(x) + qDn
(13)
(14)
where q is used to indicate the absolute value of the electronic charge. Although no direct assumptions on the non–equilibrium distribution function f (v, x) was made in the derivation of (13) and
(14), in effect, the choice of equilibrium (thermal) velocity means that the drift–diffusion equations
are only valid for very small perturbations of the equilibrium state (low fields). The validity of the
drift–diffusion equations is empirically extended by introducing field–dependent mobility µ(E) and
diffusion coefficient D(E), obtained from empirical models or detailed calculations.
3
Physical Limitations on Numerical Drift–Diffusion Schemes
The complete drift–diffusion model is based on the following set of equations
1. Current equations
Jn = qnµn E + qDn ∇n
(15)
Jp = qpµp E − qDp ∇p
(16)
2. Continuity equations
∂n
1
= ∇·Jn + Un
∂t
q
∂p
1
= − ∇·Jp + Up
∂t
q
3. Poisson’s equation
∇2 V =
+
q(n − p + NA− − ND
)
(17)
(18)
(19)
where Un and Up are the net generation–recombination rates. The continuity equations (17) and
(18) are the conservation laws for the carriers. A numerical scheme which solves the continuity
equations should
3
1. Conserve the total number of particles inside the device simulated.
2. Respect local positivity of carrier density. Negative density is unphysical.
3. Respect monotonicity of the solution (i.e. it should not introduce spurious space oscillations).
Conservative schemes are usually achieved by subdivision of the computational domain into
patches (boxes) surrounding the mesh points. The currents are then defined on the boundaries of
these elements, thus enforcing conservation (the current exiting one element side is exactly equal
to the current entering the neighboring element through the side in common). For example, on
a uniform 2–D grid with mesh size ∆, the electron continuity equation could be discretized in an
explicit form as follows
n(i, j; k + 1) − n(i, j; k)
∆t
=
J x (i + 21 , j; k) − J x (i − 12 , j; k))
q∆
1
y
J (i, j + 2 ; k) − J y (i, j − 21 ; k)
+
+ U (i, j; k)
q∆
(20)
This simple approach has certainly practical limitations, but is sufficient to exemplify the idea
behind the conservative scheme. With the present convention for positive and negative components,
it is easy to see that in the absence of generation-recombination terms, the only contributions to
the overall device current arise from the contacts. Remember that, since electrons have negative
charge, the particle flux is opposite to the current flux. The actual determination of the current
densities appearing in (20) will be discussed later.
When the equations are discretized, using finite differences for instance, there are limitations
on the choice of mesh size and time step:
1. The mesh size ∆x is limited by the Debye length.
2. The time step is limited by the dielectric relaxation time.
A mesh size must be smaller than the Debye length where one has to resolve charge variations
in space. A simple example is the carrier redistribution at an interface between two regions with
different doping levels. Carriers diffuse into the lower doped region creating excess carrier distribution which at equilibrium decays in space down to the bulk concentration with approximately
exponential behavior. The space decay constant is the Debye length
s
LD =
kB T
q2N
(21)
where N is the doping. In GaAs and Si, at room temperature the Debye length is approximately
400 Å when N ≈ 1016 cm−3 and decreases to about only 50 Å when N ≈ 1018 cm−3 .
The dielectric relaxation time is the characteristic time for charge fluctuations to decay under
the influence of the field that they produce. The dielectric relaxation time may be estimated using
tdr =
qN µ
(22)
To see the importance of respecting the limitations related to the dielectric relaxation time, imagine
to have a space fluctuation of carrier concentration which is going to relax to equilibrium with a
law
4
∂∆n
∆n(t = 0)
=−
∂t
tdr
(23)
A finite difference discretization of this equation gives at the first time step
∆n(∆t) = ∆n(0) −
∆t∆n(0)
tdr
(24)
Clearly, if ∆t > tdr , the value of ∆n(∆t) is negative, which means that the actual concentration
is evaluated to be less than the equilibrium value, and it is easy to see that the solution at higher
time steps will decay oscillating between positive and negative values of ∆n. An excessively large
∆t may lead, therefore, to nonphysical results. In the case of high mobility the dielectric relaxation
2
time can be very small. For instance, a sample of GaAs with a mobility of 6, 000 Vcm−s and doping
1018 cm−3 has approximately tdr ≈ 10−15 s, and in a simulation the time step ∆t should be chosen
to be considerably smaller.
4
Steady State Solution of Bipolar Semiconductor Equations
The general semiconductor equations may be rewritten as
∇·(∇V ) = q(n − p + NB )
∂n
∇·Jn = qU (n, p) + q
∂t
∂p
∇·Jp = −qU (n, p) + q
∂t
kB T
Jn = qµn −n∇V +
∇n
q
kB T
Jp = qµp −p∇V −
∇p
q
(25)
(26)
(27)
(28)
(29)
with NB = NA − ND . Note that the above equations are valid in the limit of small deviations from
equilibrium, since the Einstein relations have been used for the diffusion coefficient, normally valid
for low fields or large devices. The generation–recombination term U will be in general a a function
of the local electron and hole concentrations, according to possible different physical mechanisms,
to be examined later in more detail. We will consider from now on steady state, with the time
dependent derivatives set to zero.
The semiconductor equations constitute a coupled nonlinear set. It is not possible, in general,
to obtain a solution directly in one step, but a nonlinear iteration method is required. The two
more popular methods for solving the discretized equations are the Gummel’s iteration method and
the Newton’s method. It is very difficult to determine an optimum strategy for the solution, since
this will depend on a number of details related to the particular device under study.
There are in general three possible choices of variables.
1. Natural variable formulation (V, n, p)
2. Quasi–Fermi level formulation (V, φn , φp ), where the quasi–Fermi levels derive from the following definition of carrier concentration out of equilibrium (non–degenerate case)
5
q(V − φn )
kB T
q(φp − V )
p = ni exp
kB T
n = ni exp
(30)
(31)
3. Slotboom formulation (V, Φn , Φp ) where the Slotboom variables are defined as
−qφn
= ni exp
k T
B qφp
= ni exp
kB T
Φn
Φp
(32)
(33)
The Slotboom variables are therefore related to the carrier definitions (30) and (31), and the
extension to degenerate conditions is cumbersome.
Normally, there is a preference for the quasi–Fermi level formulation in steady state simulation,
and for the natural variables n and p in transient simulation.
4.1
Normalization and Scaling
For the sake of clarity, all formulae have been presented without the use of simplifications or normalization. It is however common practice to perform the actual calculation using normalized units
to make the algorithms more efficient, and in cases to avoid numerical overflow and underflow. It is
advisable to input the data in M.K.S. or practical units (the use of centimeters is for instance very
common in semiconductor practice, instead of meters) and then provide a conversion block before
and after the computation blocks to normalize and denormalize the variables. It is advisable to use
consistently scaling, rather than set certain constants to arbitrary values. In computational physics,
it is common practice to use h̄ = 1, for instance. This simplifies a lot the derivation of formulae,
but if this is not done carefully, when denormalization is applied to recover quantitative results in
M.K.S. or other appropriate units, it may be very difficult to verify which terms are multiplied by
h̄, h̄2 , h̄3 , etc. Problems may arise when formulae are taken from the published literature, since
c.g.s. and M.K.S. units are practically interchangeable, unless the dielectric constant is involved.
In c.g.s. units, the vacuum permittivity is o = 1 or 1/4π but it is about 8.85 × 10−12 [F/m] in
M.K.S. units. A simple change of centimeters into meters and grams into kilograms, will leave
many orders of magnitude unaccounted for when the dielectric permittivity is involved. According
to our experience, the majority of mistakes in computational applications prepared by beginners
are caused by improper scaling or choice of units. Great care and systematic dimensional analysis
of the formulae is suggested, at least until a good feeling is developed for the orders of magnitudes
of the variables involved in a computation.
The most common scaling factors for normalization of semiconductor equations are listed in
Table I.
4.2
Gummel’s Iteration Method
Gummel’s method solves the equations with a decoupled procedure. If we choose the quasi–Fermi
level formulation, we solve first a nonlinear Poisson’s equation. The potential obtained is substituted into the continuity equations, which are now linear, and are solved directly to conclude the
6
TABLE I - Scaling Factors
Space
Potential
Carrier concentration
Diffusion coefficient
Intrinsic Debye length (N = ni )
Extrinsic Debye length (N = Nmax )
Thermal voltage
Intrinsic concentration
Maximum doping concentration
Practical unit
Maximum diffusion coeff.
L=
q
kB T
q2 N
V ∗ = kBq T
N = ni
N = Nmax
2
D = 1 cm
s
D = Dmax
M = VD∗
R = DN
L2
L2
T = D
Mobility
Gen–Recomb
Time
iteration. The result in terms of quasi–Fermi levels is then substituted back into Poisson’s equation
until convergence is reached. In order to check for convergence, one can calculate the residuals
obtained by positioning all the terms to the left hand side of the equations and substituting the
variables with the iteration values. For exact solution values, the residuals should be zero. Convergence is assumed when the residuals are smaller than a set tolerance. The rate of convergence
of the Gummel method is faster when there is little coupling between the different equations.
The computational cost of one Gummel iteration is one matrix solution for each carrier type
plus one iterative solution for the linearization of Poisson’s equation. Note that in conditions of
equilibrium (zero bias) only the solution of Poisson’s equation is necessary, since the equilibrium
Fermi level is constant and coincides with both quasi–Fermi levels.
We give some examples of the quasi–linearization of Poisson equation, as necessary when Gummel’s method is implemented. Let’s consider the 1–D case in equilibrium first. As said earlier,
one has to solve only Poisson’s equation, since exact expressions for the carrier concentrations are
known. In the non–degenerate case, (30) and (31) are substituted into Poisson’s equation to give
d2 V
q
qV
qV
=
ni exp(−qφn ) exp
− ni exp(qφp ) exp −
+ NA − ND
(34)
2
dx
kB T
kB T
In equilibrium it is φn = φp = 0 (Fermi level taken as reference for the potential energy). Furthermore, the equation may be scaled by using the (minimum) extrinsic Debye length for the space
coordinate x, and the thermal voltage kB T /q for the potential V . Indicating with V̄ and x̄ for the
normalized potential and space coordinates
d2 V̄
ni
NA − ND
=
exp(V̄ ) − exp(−V̄ ) +
(35)
dx̄2
N
ni
The equilibrium Poisson’s equation (35) can be solved with the following quasi–linearization
procedure
1. Set initial guess for the potential V̄ .
2. Write the potential at the next iteration as V̄new = V̄ + δV , and write Poisson’s equation for
V̄new with the above substitution
d2 V̄
d2 δV
ni
NA − ND
+
=
exp(V̄ ) exp(δV ) − exp(−V̄ ) exp(−δV ) +
2
2
dx̄
dx̄
N
ni
7
(36)
3. Use the linearization exp(±δV ) ≈ 1 ± δV and discretize (36). This equation has a tridiagonal
matrix and is readily solved for δV (i).
ni 2
δV (i − 1) − 2 +
∆ x[exp(V̄ (i)) − exp(−V̄ (i))] δV (i) + δV (i + 1) =
N
(37)
n1 2
NA − ND
∆ x exp(V̄ (i)) − exp(−V̄ (i)) +
N
ni
−V̄ (i − 1) + 2V̄ (i) − V̄ (i + 1) +
4. Check for convergence. The residual of (35) is calculated and convergence is achieved if the
norm of the residual is smaller than a preset tolerance. If convergence is not achieved, return
to step 2. In practice one might simply check the norm of the error
||δV ||2 ≤ T ol
or
||δV ||∞ ≤ T ol
(38)
Note that for the solution of Poisson’s equation the boundary conditions are referenced to the
equilibrium Fermi level. One may use the separation between the Fermi level and the intrinsic
Fermi level at the contacts for the boundary conditions.
After the solution in equilibrium is obtained, the applied voltage is increased in steps ∆V ≤
kB T /q. Now the scaled nonlinear Poisson equation becomes
d2 V
ni
NA − ND
=
exp(−φn ) exp(V ) − exp(φp ) exp(−V ) +
2
dx
N
ni
(39)
where the quasi–Fermi levels are also normalized. The continuity equations, as long as the Einstein’s
relations are valid, may be written as
∂V
kB T ∂
q(V − φn )
−qµn n
+ qµn
ni exp
∂x
q ∂x
kB T
∂V
kB T
q
∂V
∂φn
−qµn n
+ qµn
n
−
∂x
q
kB T ∂x
∂x
∂φn
−qµn n
∂x q(V − φn ) ∂φn
−qµn ni exp
kB T
∂x
qV
−kB T ∂
−qφn
−qµn ni exp
exp
kB T
q
∂x
kB T
Jn =
=
=
=
=
(40)
which can be compacted, including quasi–Fermi level normalization as
Jn = an (x)
∂
exp(−φn )
∂x
(41)
∂
exp(φp )
∂x
(42)
A similar formula is obtained for the holes
Jp = ap (x)
and the continuity equations are
8
∂
∂
an (x)
exp(−φn ) = qU (x)
∂x
∂x
∂
∂
ap (x)
exp(φp ) = qU (x)
∂x
∂x
(43)
(44)
The continuity equations may be discretized with a straightforward finite difference approach
(here for simplicity with uniform mesh)
an (i+ 12 )[Φn (i+1)−Φn (i)]
∆x
an (i− 1 )[Φn (i)−Φn (i−1)]
2
−
∆x
=U
(45)
∆x
where the Slotboom variables have been used for simplicity of notation. Note that the inner derivative has been discretized with centered differences around the points (i ± 12 ) of the interleaved mesh.
Variables on the interleaved mesh must be determined very carefully, using consistent interpolation
schemes for potential and carrier density, as discussed later. The continuity equations give the
tridiagonal system
1
1
1
1
an (i − )Φn (i − 1) − [an (i + ) + an (i − )]Φn (i) + an (i + )Φn (i + 1) =
2
2
2
2
∆2 xU (i)
1
1
1
1
ap (i − )Φp (i − 1) − [ap (i + ) + ap (i − )]Φp (i) + ap (i + )Φp (i + 1) =
2
2
2
2
−∆2 xU (i)
(46)
(47)
The decoupled iteration solves now Poisson’s equation (39), initially with a guess for the quasi–
Fermi levels. The voltage distribution obtained for the previous voltage considered is normally a
good initial guess for the potential. Since the quasi–Fermi levels are inputs for Poisson’s equation,
the quasi–linearization procedure for equilibrium can be used again. The potential is then used to
update the an (i) and ap (i), and (46) and (47) are solved to provide new quasi–Fermi level values
for Poisson’s equation, and the process is repeated until convergence is reached. The generation–
recombination term depends on the electron and hole concentrations, therefore it has to be updated
at each iteration. It is possible to update the generation–recombination term also intermediately,
using the result of (46) for the electron concentration.
The examples given here to illustrate the Gummel’s approach are limited to the nondegenerate
case. If field dependent mobility and diffusion coefficients are introduced, minimal changes should
be necessary, as long as it is still justified to use of Einstein’s relations. Extension to nonuniform
mesh is left as an exercise for the reader.
In the 2–D case, the quasi–linearized Poisson’s equation becomes
ni
− 4+h
[Φn (i, j) exp(V (i, j)) + Φp (i, j) exp(−V (i, j))] δV (i, j)
N
+[δV (i − 1, j) + δV (i + 1, j) + δV (i, j − 1) + δV (i, j + 1)] =
2
4V (i, j) − V (i − 1, j) − V (i + 1, j) − V (i, j − 1) − V (i, j + 1) +
NA + NB
2 ni
h
Φn (i, j) exp(V (i, j)) + Φp (i, j) exp(−V (i, j)) +
N
ni
9
(48)
The normalized mesh size is h = ∆x = ∆y. As before, the thermal voltage kB T /q has been used to
normalize the potential V and the quasi–Fermi levels φn and φp included in the Slotboom variables
Φn,p = exp(±φn,p ).
The continuity equations with the form ∇·(a(x, y)∇Φ) = ±U (x, y) are discretized as
1
1
1
1
1
− a(i + , j) + a(i − , j) + a(i, j + ) + a(i, j − ) Φ(i, j) + a(i + , j)Φ(i + 1, j)
2
2
2
2
2
1
1
1
+a(i − , j)Φ(i − 1, j) + a(i, j + )Φ(i, j + 1) + a(i, j − )Φ(i, j − 1) = ±h2 U (i, j)
2
2
2
4.3
(49)
Newton’s Method
Newton’s method is a coupled procedure which solves the equations simultaneously, through a
generalization of the Newton–Raphson method for determining the roots of an equation. We
rewrite (25)–(27) in the residual form
WV (V, n, p) = 0 Wn (V, n, p) = 0 WV (V, n, p) = 0
(50)
Starting from an initial guess Vo , no , and po , the corrections V , ∆n, and ∆p are calculated from
the Jacobian system
∂WV
∂V
∂Wn
∂V
∂Wp
∂WV
∂n
∂Wn
∂n
∂Wp
∂n
∂V
∂WV
∂p
∂Wn
∂n
∂Wp
∂p
W
∆V V
∆n = − Wn
Wp
∆p (51)
which is obtained by Taylor expansion. The solutions are then updated according to the scheme
V (k + 1) = V (k) + ∆V (k)
n(k + 1) = n(k) + ∆n(k)
(52)
p(k + 1) = p(k) + ∆p(k)
where k indicates the iteration number. In practice, a relaxation approach is also applied to avoid
excessive variations of the solutions at each iteration step.
The system (51) has 3 equations for each mesh point on the grid. This indicates the main
disadvantage of a full Newton iteration, related to the computational cost of matrix inversion (one
may estimate that a 3N × 3N matrix takes typically 20 times longer to invert than an analogous
N × N matrix!). On the other hand convergence is usually fast for the Newton method, provided
that the initial condition is reasonably close to the solution, and is in the neighborhood where the
solution is unique. There are several viable approaches to alleviate the computational requirements
of the Newton’s method. In the Newton–Richardson approach, the Jacobian matrix in (51) is
updated only when the norm of the error does not decrease according to a preset criterion.
In general, the Jacobian matrix is not symmetric positive definite, and fairly expensive solvers
are necessary. Iterative schemes have been proposed to solve each step of Newton’s method by
reformulating (51) as
∂W
V
∂V
∂Wn
∂V
∂Wp
∂V
0
∂Wn
∂n
∂Wp
∂n
0
0
∂Wp
∂p
|
W
V
= − Wn
Wp
{z }
∆V
∆n
∆p
k+1
10
0
− 0
0
∂WV
∂n
0
0
∂WV
∂p
∂Wn
∂p
0
|
{z }
∆V
∆n
∆p
k
(53)
Since the matrix on the left hand side is lower triangular, one may solve decoupling into three
systems of equations solved in sequence. First, one solves the block of equations (again, one for
each grid point)
∂WV
∂WV
∂WV
(∆V )k+1 = −WV −
(∆n)k −
(∆p)k
∂V
∂n
∂p
and the result is used in the next block of equations
(54)
∂Wn
∂Wn
∂Wn
(∆n)k+1 = −Wn −
(∆V )k+1 −
(∆p)k
∂n
∂V
∂p
Similarly, for the third block
(55)
∂Wp
∂Wp
∂Wp
(∆p)k+1 = −Wp −
(∆V )k+1 −
(∆n)k+1
(56)
∂p
∂V
∂n
The procedure achieves a decoupling of the equations as in a block Gauss-Seidel iteration, and can
be intended as a generalization of the Gummel method. A block–SOR method is obtained if the
left hand sides are premultiplied by a relaxation parameter. This iteration procedure has better
performance if the actual variables are (V, φn , φp ).
In general, Gummel’s method is preferred at low bias because of its faster convergence and
low cost per iteration. At medium and high bias the Newton’s method becomes more convenient,
since the convergence rate of Gummel’s method becomes worse as the coupling between equations
becomes stronger at hogher bias. But since Gummel’s method has a fast initial error reduction, it
is often convenient to couple the two procedures, using Newton’s method after several Gummel’s
iterations. Remember that it is very important for the Newton’s iteration to start as close as
possible to the true solution. Close to convergence, the residual in Newton’s iteration should
decrease quadratically from one iteration to the other.
5
Generation and Recombination
The Shockley–Reed–Hall model is very often used for the generation–recombination term due to
trap levels
USRH =
np − n2i
h
τp n + ni exp
q(Et −Ei )
kB T
i
h
+ τn p + ni exp
q(Ei −Et )
kB T
i
(57)
where Et is the trap energy level involved and τn and τp are the electron and hole lifetimes. Surface
1
rates may be included with a similar formula, in which the lifetimes are substituted by sn,p
where
sn,p is the surface recombination velocity.
The Auger recombination may be accounted for by using the formula
UAug = Cn [pn2 − nn2i ] + Cp [np2 − pn2i ]
(58)
where Cn and Cp are appropriate constants. The Auger effect is for instance very relevant in the
modeling of highly doped emitter regions in bipolar transistors.
The generation process due to impact ionization can be included using the field–dependent rate
a∞
n exp
crit βn
− EnE
|Jn |
UI =
+
q
11
a∞
p exp
E crit
− pE
βp
|Jp |
(59)
6
Time-dependent simulation
The time–dependent form of the drift–diffusion equations can be used both for steady–state and
transient calculations. Steady–state analysis is accomplished by starting from an initial guess,
and letting the numerical system evolve until a stationary solution is reached, within set tolerance
limits. This approach is seldom used in practice, since now robust steady–state simulators are
widely available. It is nonetheless an appealing technique for beginners since a relatively small
effort is necessary for simple applications and elementary discretization approaches. If an explicit
scheme is selected, no matrix solutions are necessary, but it is normally the case that stability is
possible only for extremely small timesteps.
The simulation of transients requires the knowledge of a physically meaningful initial condition, which can be obtained from a steady-state calculation. The same time–dependent numerical
approaches used for steady-state simulation are suitable, but there must be more care for the
boundary conditions, because of the presence of displacement current during transients. In a transient simulation to determine the steady–state, the displacement current can be neglected because
it goes to zero when a stationary condition is reached. Therefore, it is sufficient to impose on
the contacts the appropriate potential values provided by the bias network. In a true transient
regime, however, the presence of displacement currents manifests itself as a potential variation at
the contacts, superimposed to the bias, which depends on the external circuit in communication
with the contacts. Neglect of the displacement current in a transient is equivalent to the application
of bias voltages using ideal voltage generators, with zero internal impedance. In such a situation,
the potential variations due to displacement current drop across a short circuit, and are therefore
cancelled. In this arrangement, one will observe the shortest possible switching time attainable
with the structure considered, but in practice an external load and parasitics will be present, and
the switching times will be normally longer. A simulation neglecting displacement current effects
may be useful to assess the ultimate speed limits of a device structure.
When a realistic situation is considered, it is necessary to include a displacement term in the
current equations. It is particularly simple to deal with a 1–D situation. Consider a 1–D device
with length W and a cross–sectional area A. The total current flowing in the device is
∂E(x, t)
(60)
∂t
The displacement term makes the total current constant at each possition x. This property can be
exploited to perform an integration along the device
ID (t) = In (x, t) + A
1
ID (t) =
W
Z W
In (x, t)dx +
0
A ∂V ∗
W ∂t
(61)
where V ∗ (t) is the total voltage drop across the structure, with the ground reference voltage applied
at x = W . The term A
W is called cold capacitance. The 1–D device, therefore, can be studied as
the parallel of a current generator and of the cold capacitance which is in parallel with the (linear)
load circuit. At every time step, V ∗ has to be updated, since it depends on the charge stored by
the capacitors.
To illustrate the procedure, consider a simple Gunn diode in parallel with an RLC resonant
load containing the bias source. Calling Co the parallel of cold and load capacitance, it is
I(t) = Co
dV ∗ (t)
+ Io (t)
dt
12
(62)
where Io (t) is the particle current given by the first term on the right hand side of (61), calculated
at the given time step with drift–diffusion (or any other suitable scheme). It is also
V ∗ (t)
V ∗ (t) − Vb
−
dt
R
L
Upon time differencing this last equation, with the use of finite differences we obtain
I(t) = −
Z
∆t
Co
V ∗ (t + ∆t) − V ∗ (t)
∆t
I(t + ∆t) = I(t) −
− [V ∗ (t) − Vb ]
R
L
V ∗ (t + ∆t) = V ∗ (t) + [I(t) − Io (t)]
(63)
(64)
(65)
This set of difference equations allows one to update the boundary conditions for Poisson’s equation
at every time step to fully include displacement current.
A robust approach for transient simulation should be based on the same numerical apparatus
established for purely steady-state models. It is usually preferred to use fully implicit schemes,
which require a matrix solution at each iteration, because the choice of the timestep is more
likely to be limited by the physical time constants of the problem rather than by stability of the
numerical scheme. In order to estimate the timestep limits, let’s assume a typical electron velocity
v = 107 cm/s and a spatial mesh ∆x = 0.01µm. The C.F.L. condition necessary to resolve correctly
a purely drift process on this mesh requires ∆t ≤ ∆x/v = 10−15 s. As calculated earlier, this value is
not too far from typical values of the dielectric relaxation time in practical semiconductor structures.
When dealing with unipolar devices, as often used in many microwave applications, it is possible to formulate very simple time–dependent drift–diffusion models, which can be solved with
straightforward finite difference techniques and are suitable for small student projects. If we can
neglect the generation–recombination effects, the 1–D unipolar drift–diffusion model is reduced to
the following system of equations
∂n
∂t
2
d V
dx2
d
d
d
[nvd (E)] +
D(E) n
dx
dx
dx
q(n − ND )
= −
(66)
=
(67)
where vd (E) = −µn (E)E is the drift velocity. There are two physical processes involved: drift
(advection) expressed by the first term on the right hand side of (66), and diffusion described
by the second term. The continuity equation (66) is an admixture of competing hyperbolic and
parabolic behavior whose relative importance depends on the local electric field strength.
The system (66) and (67) can be used for both transient or steady state conditions if the
simulation is run ∂n/∂t = 0. A basic simple algorithm could consists of the following steps
1. Guess the carrier distribution n(x).
2. Solve Poisson’s equation to obtain the field distribution.
3. Compute one iteration of the discretized continuity equation with time step ∆t. v(E) and
D(E) are updated according to the local field value.
4. Check for convergence. If convergence is obtained, stop. Otherwise, go back to step (2)
updating the charge distribution.
13
This is an uncoupled procedure, since (66) and (67) are not solved simultaneously. Usually,
explicit methods are used for computational speed. The time step must respect the limitations due
to the C.F.L. condition (related to the advective component) and to the dielectric relaxation time.
A simple discretization scheme could employ an explicit finite difference approach
n(i; k + 1) = n(i; k) +
+
∆t
{[vd (i − 1; k)n(i − 1; k) − vd (i; k)n(i; k)]
∆x
1
D(i; k)[n(i − 1; k) − 2n(i; k) + n(i + 1; k)]};
∆x
n(i; k + 1) = n(i; k) +
vd < 0
(68)
∆t
{[vd (i; k)n(i; k) − vd (i + 1; k)n(i + 1; k)]
∆x
1
D(i; k)[n(i − 1; k) − 2n(i; k) + n(i + 1; k)]};
vd > 0
(69)
∆x
where we have introduced upwinding for the drift term and we have assumed that the diffusion
coefficient is slowly varying in space. There are of course many other possible explicit and implicit
discretizations. Such simple finite difference approaches are in general a compromise which cannot
provide at one time an optimal treatment of both advective and diffusive components. Because
of spatially varying drift velocity, spurious diffusion and dispersion are present. This could be
mitigated by using a nonuniform grid discretization, where the mesh size is locally adapted to
achieve vd = ∆x/∆t everywhere, which would involve interpolation to the new grid-points. The
discretization for a diffusive process is better behaved with a fully implicit scheme (if the CrankNicholson approach is used, one needs to make sure that spurious oscillations in the solution do not
develop). On the other hand, the fully implicit algorithm for advection is not conservative. From
these conflicting requirements, it emerges that it would be beneficial to split the drift and diffusion
processes, and apply an optimal solution procedure to each. There are 1–D situations where this
is known to be nearly exact. In well known experiments, a small concentration of excess carriers is
generated in a semiconductor sample with a uniform electric field, and the motion of the centroid of
the carrier envelope can be studied independently of the diffusive spread of the spatial distribution
around the centroid itself. For an initial Gaussian distribution in space, a simple analytical solution
shows that drift and diffusion can be treated as a sequential process, each using the total duration
of the observation as simulation time. In analogy with this, the 1–D continuity equation can be
solved in two steps, for instance
+
∆t
;
vd < 0
(70)
∆x
n(j, i + 1) = n∗ (j, i + 1) + D(j)[n(j − 1, i + 1) − 2n(j, i − 1)
∆t
+n(j + 1, i + 1)] 2
(71)
∆ x
where again a simple explicit upwinding scheme is used for the drift, while a fully implicit scheme
is used for the diffusion.
n∗ (j, i + 1) = n(j, i) + vd (j)[n(j − 1, i) − n(j, i))
7
Scharfetter-Gummel approximation
The discretization of the continuity equations in conservation form, requires the determination of
the currents on the mid–points of mesh lines connecting neighboring grid nodes. Since the solutions
are accessible only on the grid nodes, interpolation schemes are needed to determine the currents.
14
For consistency with Poisson’s equation, it is common to assume that the potential varies
linearly between two neighboring nodes. This is equivalent to assume a constant field along the
mesh lines, and the field at the mid–point is obtained by centered finite differences of the potential
values. In order to evaluate the current, it is also necessary to estimate the carrier density at the
mid–points. The simplest approximation which comes to mind is to also assume a linear variation
of the carrier density, by taking the arithmetic average between to neighboring nodes. This simple
approach is only acceptable for very small potential variation between the nodes, and indeed is
exact only if the field between two nodes is zero, which implies the same exact carrier density on
the two points.
In order to illustrate this, let’s consider a 1–D mesh where we want to discretize the electron
current
dψ
dn
+ qDn
(72)
dx
dx
Here, the field is explicitly expressed by the derivative of the potential. The discretization on the
mid–point of the mesh line between nodes xi and xi+1 is given by
Jn = qµn n −
ψi+1 − ψi
ni+1 − ni
+ qDn
2
2
∆x
∆x
In the simple approach indicated above, the carrier density ni+ 1 is expressed as
Ji+ 1 = − qµn ni+ 1
(73)
2
ni+1 + ni
(74)
2
In (73), the assumed linearity of the potential between meshes, is implied by the use of the
centered finite differences to express the field on the mid–point. We can now rewrite (73) including
the approximation in (74) as
ni+ 1 ≈
2
Ji+ 1
=
2
ni+1
−q
µn ψi+1 − ψi
Dn
+ q
2
∆x
∆x

−

 µn ψi+1 − ψi
Dn 

ni 
+ q
q 2
∆x}
|
{z∆x }
| {z
a
(75)
b
If we assume a condition where Jn = 0 (equilibrium) and a >> b (negligible diffusion) it is
easy to see that positivity of the carrier density is not guaranteed, since the solution oscillates as
ni+1 ≈ −ni . Also, it can be shown that for stability we need to have ψi+1 − ψi > 2kB T /q, which
requires very small mesh spacing to be verified.
The approach by Scharfetter and Gummel (1969) has provided an optimal solution to this
problem, although the mathematical properties of the proposed scheme have been fully recognized
much later. We consider again a linear potential variation between neighboring mesh points, which
is consistent with the use of finite differences to express the field. We express the current in the
interval [xi ; xi+1 ] as a truncated expansion about the value at the mid–point xi+ 1
2
∂
Jn (x)
(76)
2
2 ∂x
From (76) we obtain a first order differential equation for Jn which can be solved to provide n(x)
in the mesh interval, using as boundary conditions the values of carrier density ni and ni+1 . We
obtain
Jn (x) = Jn (xi+ 1 ) + (x − xi+ 1 )
15
n(x) = [1 − g(x, ψ)] ni + g(x, ψ)ni+1 ;
x ∈ [xi ; xi+1 ]
(77)
where g(x, ψ) is the growth function
ψi+1 − ψi x − xi
1 − exp
kB T /q
∆x
g(x, ψ) =
ψi+1 − ψi
/ 1 − exp
kB T /q
(78)
The result in (77) can be used to evaluate n(xi+ 1 ) for the discretization of the current in (73). It
2
is easy to see that only when ψ(i + i) − ψ(i) = 0 we have
1
1
ni + ni+1
n1+ 1 = (1 − )ni + ni+1 =
(79)
2
2
2
2
The continuity equation can be easily discretized on rectangular uniform and nonuniform meshes
using the above results for the currents, because the mesh lines are aligned exactly.
For a more in depth reading on advanced aspects of drift–diffusion simulation, we recommend
to consult the book by S. Selberherr, cited below.
Some classic references:
H.K. Gummel, “A self–consistent iterative scheme for one–dimensional steady state transistor calculation,” IEEE Transactions on Electron Devices, vol. ED-11, pp.455–465, 1964.
A. DeMari, “An accurate numerical steady state one–dimensional solution of the p–n junction,”
Solid–state Electronics, vol. 11, pp. 33–59, 1968.
A. DeMari, “An accurate numerical one–dimensional solution of the p–n junction under arbitrary
transient conditions,” Solid–state Electronics, vol. 11, pp. 1021–2053, 1968.
S. Selberherr, Analysis and Simulation of Semiconductor Devices, Springer–Verlag, 1984.
D. L. Scharfetter and D. L. Gummel, “Large signal analysis of a Silicon Read diode oscillator,”
IEEE Transaction on Electron Devices, vol. ED-16, pp.64–77, 1969.
J. W. Slotboom, “Iterative scheme for 1 and 2–dimensional d.c. transistor simulation,” Electronics
Letters, vol. 5, pp. 677–678, 1969.
J.W. Slotboom, “Computer–aided two–dimensional analysis of bipolar transistors,” IEEE Transactions on Electron Devices, vol. ED-20, pp. 669–679, 1973.
16
```