Transcript Semi-Classical Transport Theory Outline:  What is Computational Electronics?  Semi-Classical Transport Theory  Drift-Diffusion Simulations  Hydrodynamic Simulations  Particle-Based Device Simulations  Inclusion of Tunneling.

Semi-Classical Transport
Theory
Outline:
 What is Computational Electronics?
 Semi-Classical Transport Theory
 Drift-Diffusion Simulations
 Hydrodynamic Simulations
 Particle-Based Device Simulations
 Inclusion of Tunneling and Size-Quantization Effects in Semi-Classical Simulators
 Tunneling Effect: WKB Approximation and Transfer Matrix Approach
 Quantum-Mechanical Size Quantization Effect
 Drift-Diffusion and Hydrodynamics: Quantum Correction and Quantum
Moment Methods
 Particle-Based Device Simulations: Effective Potential Approach
 Quantum Transport
 Direct Solution of the Schrodinger Equation (Usuki Method) and Theoretical
Basis of the Green’s Functions Approach (NEGF)
 NEGF: Recursive Green’s Function Technique and CBR Approach
 Atomistic Simulations – The Future
 Prologue
Semi-Classical Transport Theory
 It is based on direct or approximate solution of the
Boltzmann Transport Equation for the semi-classical
distribution function f(r,k,t)
f k 1
F
V
  k E  r f k   k f k  3  dk  f k  1  f k  k k  f k 1  f k    kk  
t
8
which gives one the probability of finding a particle in
region (r,r+dr) and (k,k+dk) at time t
 Moments of the distribution function give us information
about:
 Particle Density
 Current Density
 Energy Density
D. K. Ferry, Semiconductors, MacMillian, 1990.
Semi-Classical Transport Approaches
1.Drift-Diffusion Method
2.Hydrodynamic Method
3.Direct Solution of the Boltzmann Transport
Equation via:
Particle-Based Approaches – Monte Carlo
method
Spherical Harmonics
Numerical Solution of the Boltzmann-Poisson
Problem
C. Jacoboni, P. Lugli: "The Monte Carlo Method for Semiconductor Device Simulation“,
in series "Computational Microelectronics", series editor: S. Selberherr; Springer, 1989, ISBN: 3211-82110-4.
1. Drift-Diffusion Approach
Constitutive Equations
Poisson
  V    p  n  N D  N A 
Continuity Equations
n
1
   Jn  Un
t
q
p
1
    J p U p
t
q
Current Density Equations
S. Selberherr: "Analysis and Simulation of
Semiconductor Devices“, Springer, 1984.
dn
J n  qn( x) n E ( x)  qDn
dx
dn
J p  qp ( x)  p E ( x)  qD p
dx
Numerical Solution Details
 Linearization of the Poisson equation
 Scharfetter-Gummel Discretization of the
Continuity equation
 De Mari scaling of variables
 Discretization of the equations
 Finite Difference – easier to implement but requires
more node points, difficult to deal with curved interfaces
 Finite Elements – standard, smaller number of node
points, resolves curved surfaces
 Finite Volume
Linearized Poisson Equation
φ→φ + δ where δ= φnew - φold
Finite difference discretization:
 Potential varies linearly between mesh points
 Electric field is constant between mesh points
Linearization → Diagonally-dominant
coefficient matrix A is obtained

d 2V new eni V old /VT

e
2
dx
VT





eni V old /VT V old /VT
eni
d 2V new
V old / VT
V old / VT

e
e
 C / ni   e
e
2
dx


eni V old /VT V old /VT
V old / VT
new
e
V 
e
e
 C / ni 





eni V old /VT V old /VT old
e
e
V
VT
  V new  V old

Scharfetter-Gummel Discretization of the
Continuity Equation
 Electron and hole densities n and p vary exponentially
between mesh points → relaxes the requirement of
using very small mesh sizes
 The exponential dependence of n and p upon the
potential is buried in the Bernoulli functions
Din1/ 2
2
 Din1/ 2
 Vi 1  Vi 
B
 ni 1   2
 VT 
 
 Vi  Vi 1  Din1/ 2
B
 2
 VT  
 Vi  Vi 1  
Din1/ 2
B
  ni  2

 VT  
 V V 
B  i 1 i  ni 1  U i
 VT 
Din1/ 2
2
 Din1/ 2
 Vi  Vi 1 
B
 pi 1   2
 VT 
 
 Vi 1  Vi  Din1/ 2
B
 2
 VT  
 Vi 1  Vi 
Din1/ 2
B
 pi  2

 VT 
 V V 
B  i i 1  pi 1  U i
 VT 
x
Bernoulli function: B( x)  x
e 1
Scaling due to de Mari
Variable
Scaling Variable
Space
Intrinsic Debye length (N=ni)
Formula
 k BT
L
q2 N
Extrinsic Debye length (N=Nmax)
Potential
Thermal voltage
V* 
Carrier concentration
Intrinsic concentration
N=ni
Maximum doping concentration
N=Nmax
Practical unit
cm 2
D 1
s
Diffusion coefficient
k BT
q
Maximum diffusion coefficient
D = Dmax
Mobility
M 
D
V*
Generation-Recombination
R
DN
L2
T 
L2
D
Time
Numerical Solution Details
Governing
Equations
ICS/BCS
Continuous
Solutions
Discretization
Finite-Difference
Finite-Volume
System of
Algebraic
Equations
Equation
(Matrix)
Solver
Discrete
Nodal
Values
Tridiagonal
φi (x,y,z,t)
SOR
p (x,y,z,t)
Finite-Element
Gauss-Seidel
Spectral
Krylov
Boundary Element
Multigrid
Approximate
Solution
n (x,y,z,t)
Hybrid
D. Vasileska, EEE533 Semiconductor Device and Process Simulation Lecture Notes,
Arizona State University, Tempe, AZ.
Numerical Solution Details
Poisson solvers:
 Direct
Gaussian Eliminatioln
LU decomposition
 Iterative
 Mesh Relaxation Methods
• Jacobi, Gauss-Seidel, Successive over-Relaxation
Advanced Iterative Solvers
• ILU, Stone’s strongly implicit method, Conjugate gradient
methods and Multigrid methods
G. Speyer, D. Vasileska and S. M. Goodnick, "Efficient Poisson solver for semiconductor
device modeling using the multi-grid preconditioned BICGSTAB method", Journal of
Computational Electronics, Vol. 1, pp. 359-363 (2002).
Complexity of linear solvers
Time to solve
model problem
(Poisson’s
equation) on
regular mesh
n1/2
n1/3
2D
3D
Sparse Cholesky:
O(n1.5 )
O(n2 )
CG, exact
arithmetic:
O(n2 )
O(n2 )
CG, no precond:
O(n1.5 )
O(n1.33 )
CG, modified IC:
O(n1.25 )
O(n1.17 )
CG, support trees:
Multigrid:
O(n1.20 ) -> O(n1+ ) O(n1.75 ) -> O(n1.31 )
O(n)
O(n)
LU Decomposition
• If LU decomposition exists, then for a tri-diagonal matrix A,
resulting from the finite-difference discretization of the 1D
Poisson equation, one can write
 a1
b
 2




c1
a2
...
c2
...
bn 1
...
a n 1
bn
where 1  a1 ,
  1
 
  2

cn 1  
a n  
k 
1
 3 ...
...
...
n
 1




1 
bk
,  k  a k   k ck 1 ,
 k 1
c1
2
c2
...



...

... cn 1 
 n 
k  2,3,..., n
Then, the solution is found by forward and back substitution:
g1  f1 ,
gn
xn 
,
n
gi  f i  i gi 1 , i  2,3,...,n,
gi  ci xi 1
xi 
, i  n  1,...,2,1
i
Numerical Solution Details of the Coupled
Equation Set
Solution Procedures:
 Gummel’s Approach – when the constitutive
equations are weekly coupled
 Newton’s Method – when the constitutive
equations are strongly coupled
 Gummel/Newton – more efficient approach
D. Vasileska and S. M. Goodnick, Computational Electronics, Morgan & Claypool,
2006.
D. Vasileska, S. M. Goodnick and G. Klimeck, Computational Electronics: From
Semiclassical to Quantum Transport Modeling, Taylor & Francis, in press.
Schematic Description of the Gummel’s
Approach
Original Gummel’s scheme
initial guess
of the solution
solve
Poisson’s eq.
n
converged?
Modified Gummel’s scheme
initial guess
of the solution
Solve Poisson’s eq.
Electron eq.
Hole eq.
n
y
y
Solve electron eq.
Solve hole eq.
Update
generation rate
y
n
converged?
y
converged?
n
converged?
Constraints on the MESH Size and TIME
Step
The time step and the mesh size may correlate
to each other in connection with the numerical
stability:
 The time step t must be related to the plasma
frequency
e2n
p 
sm *
 The Mesh size is related to the Debye length
LD 
 sVT
qN max
 Large devices where the velocity of the carriers is
smaller than the saturation velocity
x 10
25 nm FD SOI nMOSFET
90 nm FD SOI nMOSFET
5
5
SOI nMOSFETx 10
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
x 10
100 nm
140 nm
45 nm FD SOI nMOSFET
60 nm FD SOI nMOSFET
120 nm FD SOI nMOSFET
140 nm FD SOI nMOSFET
100
5 nm FD SOI nMOSFET
5
5
x 10
2.5 2.5 x 10 2.5 x 10
channel
Tgate=300K
source
drain
Tgate=400K
2
Tgate=600K
2
2
1.5 1.5
1
1
0.5 0.5
Velocity (m/s)
25 nm
nm5FD
Velocity (m/s)
0
When is the Drift-Diffusion Model Valid?
5
2.5
2
180 nm FD SO
5
x 10 x 10
2.5 Tgate=300K
Tgate=400K
Tgate=300K
Tgate=600K
Tgate=400K
2
Tgate=600K
5
2.5
x 10
2
1.5
1.5
1.5
1.5
1
1
1
1
0.5
0.5
0.5
0.5
0
0
0
0
0
0
0
0
0
1.35
0.6
2.4
1.41.2 2.8 1.8
4.2
0
0.8
1.6
2.4
0 2.5 0.9 5 1.8 7.5 2.70
0 0.451 0 0.91.2
2
3 0 3.6 0
-8
-7 -7
-7 -7
-7
-7
-7
x (m)
x (m)
x(m)x (m)
(m)
x (m)
x (m)
x (m) x x
x 10
10
x 10
x 10
x 10
x 10
x 10
x 10
 The validity of the method can be extended for velocity
saturated devices as well by introduction of electric field
dependent mobility and diffusion coefficient:
Dn(E) and μn(E)
1.8
x (m
What Contributes to The Mobility?
Scattering Mechanisms
Defect Scattering
Crystal
Defects
Neutral
Impurity
Carrier-Carrier Scattering
Alloy
Ionized
Lattice Scattering
Intervalley
Intravalley
Acoustic
Deformation
potential
Piezoelectric
Optical
Nonpolar
Acoustic
Polar
D. Vasileska and S. M. Goodnick, "Computational Electronics", Materials Science
and Engineering, Reports: A Review Journal, Vol. R38, No. 5, pp. 181-236 (2002)
Optical
Mobility Modeling
Mobility modeling can be separated in
three parts:
 Low-field mobility characterization for bulk
or inversion layers
 High-field mobility characterization to
account for velocity saturation effect
 Smooth interpolation between the low-field
and high-field regions
Silvaco ATLAS Manual.
(A) Low-Field Models for Bulk Materials
Phonon scattering:
- Simple power-law dependence of the
temperature
- Sah et al. model:
acoustic + optical and intervalley phonons
combined via Mathiessen’s rule
Ionized impurity scattering:
- Conwell-Weiskopf model
- Brooks-Herring model
(A) Low-Field Models for Bulk
Materials (cont’d)
Combined phonon and ionized impurity scattering:
- Dorkel and Leturg model:
temperature-dependent phonon scattering +
ionized impurity scattering + carrier-carrier
interactions
- Caughey and Thomas model:
temperature independent phonon scattering +
ionized impurity scattering
(A) Low-Field Models for Bulk
Materials (cont’d)
- Sharfetter-Gummel model:
phonon scattering + ionized impurity scattering
(parameterized expression – does not use the
Mathiessen’s rule)
- Arora model:
similar to Caughey and Thomas, but with
temperature dependent phonon scattering
(A) Low-Field Models for Bulk
Materials (cont’d)
Carrier-carrier scattering
- modified Dorkel and Leturg model
Neutral impurity scattering:
- Li and Thurber model:
mobility component due to neutral impurity
scattering is combined with the mobility due to
lattice, ionized impurity and carrier-carrier
scattering via the Mathiessen’s rule
(B) Field-Dependent Mobility
The field-dependent mobility model provides smooth transition
between low-field and high-field behavior
 (E) 
0
 1 / 
  E
1   0  
  vsat  
 = 1 for electrons
 = 2 for holes
vsat is modeled as a temperature-dependent quantity:
2.4  10 7
vsat (T ) 
cm/s
 TL 
1  0.8 exp

 600 
(C) Inversion Layer Mobility Models
 CVT model:
 combines acoustic phonon, non-polar optical
phonon and surface-roughness scattering (as
an inverse square dependence of the
perpendicular electric field) via Mathiessen’s
rule
 Yamaguchi model:
 low-field part combines lattice, ionized impurity
and surface-roughness scattering
 there is also a parametric dependence on the
in-plane field (high-field component)
(C) Inversion Layer Mobility Models
(cont’d)
Shirahata model:
uses Klaassen’s low-field mobility model
takes into account screening effects into the
inversion layer
has improved perpendicular field dependence
for thin gate oxides
Tasch model:
the best model for modeling the mobility in
MOS inversion layers; uses universal mobility
behavior
Generation-Recombination Mechanisms
Classification
One step
(Direct)
Two Energy-level
particle consideration
•
•
•
•
• Shockley-Read-Hall
(SRH) generationrecombination
• Surface generationrecombination
Two-step
(indirect)
Impact
ionization
Three
particle
Auger
Photogeneration
Radiative recombination
Direct thermal generation
Direct thermal recomb.
Pure generation process
•
•
•
•
Electron emission
Hole emission
Electron capture
Hole capture
Hydrodynamic Modeling
In small devices there exists nonstationary transport and carriers are
moving through the device with velocity
larger than the saturation velocity
 In Si devices non-stationary transport occurs
because of the different order of magnitude of
the carrier momentum and energy relaxation
times
 In GaAs devices velocity overshoot occurs due
to intervalley transfer
T. Grasser (ed.): "Advanced Device Modeling and Simulation“, World Scientific
Publishing Co., 2003, ISBN: 9-812-38607-6 M.
M. Lundstrom, Fundamentals of Carrier Transport, 1990.
Velocity Overshoot in Silicon
kz
7
ky
g
kx
Drift velocity [cm/s]
2.5x10
1 kV/cm
5 kV/cm
10 kV/cm
50 kV/cm
7
2x10
7
1.5x10
1x107
6
5x10
0
6
-5x10
0
0.5
1
f
1.5
2
2.5
3
3.5
4
time [ps]
Scattering mechanisms:
• Acoustic deformation potential scattering
• Zero-order intervalley scattering (f and gphonons)
• First-order intervalley scattering (f and gphonons)
X. He, MS Thesis, ASU, 2000.
Energy [eV]
0.25
0.2
1 kV/cm
5 kV/cm
10 kV/cm
50 kV/cm
0.15
0.1
0.05
0
0
0.5
1
1.5
2
2.5
time [ps]
3
3.5
4
How is the Velocity Overshoot Accounted
For?
In Hydrodynamic/Energy balance
modeling the velocity overshoot effect is
accounted for through the addition of an
energy conservation equation in addition
to:
Particle Conservation (Continuity Equation)
Momentum (mass) Conservation Equation
Hydrodynamic Model due to Blotakjer
• Constitutive
More convenient set
of balance equations
is in terms
Equations:
Poisson
+ of n, vd
and w:
 n 
 n 
   (nvd )   
t
 t coll
 vd 
v
2
1

2
  d  (m * vd ) 
 nw  nm * v d 
t
m*
3nm * 
2


eE  ( vd ) 


m *  t coll
2

m * v d 
 w 
2
 
  vd  w 
   nvd 
  w 



t
3n 
kB 
2



 (w ) 
 eE  vd  

 t coll
Closure
• To have a closed set of equations, one either:
(a) ignores the heat flux altogether
(b) uses a simple recipe for the calculation of the heat flux:
nq   T ,
5k B2 nT

2m * v (w )
• Substituting T with the density of the carrier energy, the
momentum and energy balance equations become:
 npd 
2 
1
 (npd ) 
2
   (nvd pd )   nw  nm * v d   enE  

t
3 
2

 t co

 nw 
2 
 
1
2 
   nvd w   nvd 
  w  m * v d 
t
3
k
2
Momentum Relaxation Rate
• The momentum rate is determined by a steady-state MC
calculation in a bulk semiconductor under a uniform bias
electric field, for which:
vd
eE  vd


t
m *  t
eE


  p (w )vd  0

m*
coll
eE
 p (w ) 
m * vd
K. Tomizawa, Numerical Simulation Of
Submicron Semiconductor Devices.
Energy Relaxation Rate
• The emsemble energy relaxation rate is also determined by a
steady-state MC calculation in a bulk semiconductor under a
uniform bias electric field, for which:
w
 w 
 eE  vd  
 eE  vd  w (w  w 0 )  0

t
 t coll
eE  vd
w (w )  
w  w0
K. Tomizawa, Numerical Simulation Of
Submicron Semiconductor Devices.
Validity of the Hydrodynamic Model
tox
Gate oxide
Source
90 nm device
tsi
Drain
Lgate
2.5
LD
tBOX
BOX
feature
14 nm
25 nm
90 nm
Tox
1 nm
1.2 nm
1.5 nm
VDD
1V
1.2 V
1.4 V
Overshoot
EB/HD
233% / 224%
139% / 126%
31% /21%
Overshoot EB/DD
with series resistance
153%/96%
108%/67%
39%/26%
Drain Current [mA/um]
LS
2
1.5
1
DD
EB
HD
DD SR
EB SR
HD SR
0.5
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
SR = series resistance
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1.4
Validity of the Hydrodynamic Model
25 nm device
tox
Gate oxide
7
tsi
Drain
LS
Lgate
DD
HD
EB
DD SR
EB SR
HD SR
6
2.5
LD
tBOX
BOX
Drain Current [mA/um]
Drain Current [mA/um]
Source
5
2
4
3
1.5
feature
14 nm
25 nm
90 nm
Tox
1 nm
1.2 nm
1.5 nm
VDD
1V
1.2 V
1.4 V
Overshoot
EB/HD
233% / 224%
139% / 126%
31% /21%
Overshoot EB/DD
with series resistance
153%/96%
2
1
DD
EB
HD
DD SR
EB
1 SR
HD SR
1
0.5
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
108%/67%
39%/26%
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
SR = series resistance
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1.2
1.4
Validity of the Hydrodynamic Model
14 nm device
12
Source
tsi
Drain
LS
Lgate
Drain Current [mA/um]
Drain[mA/um]
Current [mA/um]
Drain Current
10
7
tox
Gate oxide
DD
EB
HD
DD SR
DD
EB SR
HD
HD SR
EB
DD SR
EB SR
HD SR
8
6
2.5
LD
tBOX
BOX
6
5
2
4
4
3
1.5
2
feature
14 nm
25 nm
90 nm
Tox
1 nm
1.2 nm
1.5 nm
VDD
1V
1.2 V
1.4 V
Overshoot
EB/HD
233% / 224%
139% / 126%
31% /21%
Overshoot EB/DD
with series resistance
153%/96%
2
0
1
0
0.2
0.4
1
0.5
0
DD
0.8
EB
HD
DD SR
EB
1 SR
HD SR
0.6
Drain Voltage [V]
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
108%/67%
39%/26%
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
SR = series resistance
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1
1.2
1.4
Failure of the Hydrodynamic Model
14 nm
14
25 nm
0.3 ps
8
10
7
6
1020 cm-3
2
0
0
0.2
0.2 ps
6
0.3 ps
0.1 ps
0.4
0.2 ps
2
5
0.1 ps
0.1 ps
1019 cm4-3
4
90 nm
0.3 ps
0.2 ps
Drain Current [mA/um]
8
Drain Current [mA/um]
Drain Current [mA/um]
12
1020 cm-3
3
2
0.6
0.8
1
Drain Voltage [V]1
0
0
0.2
0.4
1.5
10
20
-3
cm
10
19
-3
cm
1
1019 cm-3
0.5
0.6
0.8
1
1.2
Drain Voltage [V]
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1.4
Failure of the Hydrodynamic Model
14 nm
14
25 nm
0.3 ps
8
10
7
6
1020 cm-3
2
0
0
0.2
0.2 ps
6
0.3 ps
0.1 ps
0.4
0.2 ps
2
5
0.1 ps
0.1 ps
1019 cm4-3
4
90 nm
0.3 ps
0.2 ps
Drain Current [mA/um]
8
Drain Current [mA/um]
Drain Current [mA/um]
12
1020 cm-3
3
2
0.6
0.8
1
Drain Voltage [V]1
0
0
0.2
0.4
1.5
10
20
-3
cm
10
19
-3
cm
1
1019 cm-3
0.5
0.6
0.8
1
1.2
Drain Voltage [V]
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1.4
Failure of the Hydrodynamic Model
14 nm
14
25 nm
0.3 ps
8
10
7
6
1020 cm-3
2
0
0
0.2
0.2 ps
6
0.3 ps
0.1 ps
0.4
0.2 ps
2
5
0.1 ps
0.1 ps
1019 cm4-3
4
90 nm
0.3 ps
0.2 ps
Drain Current [mA/um]
8
Drain Current [mA/um]
Drain Current [mA/um]
12
1020 cm-3
3
2
0.6
0.8
1
Drain Voltage [V]1
0
0
0.2
0.4
1.5
10
20
-3
cm
10
19
-3
cm
1
1019 cm-3
0.5
0.6
0.8
1
1.2
Drain Voltage [V]
0
0
0.2
0.4
0.6
0.8
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
1
1.2
1.4
Summary
 Drift-Diffusion model is good for large MOSFET devices,
BJTs, Solar Cells and/or high frequency/high power
devices that operate in the velocity saturation regime
 Hydrodynamic model must be used with caution when
modeling devices in which velocity overshoot, which is a
signature of non-stationary transport, exists in the device
 Proper choice of the energy relaxation times is a
problem in hydrodynamic modeling