Transcript Tutorial: From Semi-Classical to Quantum Transport

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
Direct Solution of the Boltzmann
Transport Equation
 Particle-Based Approaches
 Spherical Harmonics
 Numerical Solution of the Boltzmann-Poisson
Problem
In here we will focus on Particle-Based
(Monte Carlo) approaches to solving the
Boltzmann Transport Equation
Ways of Solving the BTE Using MCT
Single particle Monte Carlo Technique
Follow single particle for long enough time to
collect sufficient statistics
Practical for characterization of bulk materials
or inversion layers
Ensemble Monte Carlo Technique
MUST BE USED when modeling
SEMICONDUCTOR DEVICES to have the
complete self-consistency built in
Carlo Jacoboni and Lino Reggiani, The Monte Carlo method for the solution of charge
transport in semiconductors with applications to covalent materials, Rev. Mod. Phys. 55, 645
- 705 (1983).
Path-Integral Solution to the BTE
 The path integral solution of the Boltzmann
Transport Equation (BTE), where L=Nt and
tn=nt, is of the form:
N 1
f N (t )  t  f m ( p ') Seff ( p ', p  eE ( N  m)t )e  ( N  m ) t
m 0
g m ( p  eE ( N  m)t )
K. K. Thornber and Richard P. Feynman,
Phys. Rev. B 1, 4099 (1970).
The two-step procedure is then found by
using N=1, which means that t=t, i.e.:
f1 (t )  t  f0 ( p ')Seff ( p ', p  eEt )et
p'
g0 ( p  eEt )
Intermediate function that describes
the occupancy of the state (p+eEt)
at time t=0, which can be changed
due to scattering events (SCATTER)
+
Integration over a trajectory,
i.e.probability that no scattering
occurred within time integral t
(FREE FLIGHT)
Monte Carlo Approach to Solving the
Boltzmann Transport Equation
 Using path integral formulation to the
BTE one can show that one can
decompose the solution procedure into
two components:
1. Carrier free-flights that are interrupted by
scattering events
2. Memory-less scattering events that change
the momentum and the energy of the particle
instantaneously
Particle Trajectories in Phase Space
ky
x
x
-e
x
x
x
y
Ex
x
kx
x x
x
x
Particle trajectories in k-space and real space
x
Carrier Free-Flights
 The probability of an electron scattering in a small time interval dt is
(k)dt, where (k) is the total transition rate per unit time. Time
dependence originates from the change in k(t) during acceleration by
external forces
k t   k 0  eE  v  Bt / 
where v is the velocity of the particle.
 The probability that an electron has not scattered after scattering at t = 0
is:
t
Pn (t ) 
  dt  k t  
e 0
 It is this (unnormalized) probability that we utilize as a non-uniform
distribution of free flight times over a semi-infinite interval 0 to . We
want to sample random flight times from this non-uniform distribution
using uniformly distributed random numbers over the interval 0 to 1.
Generation of Random Flight Times
Hence, we choose a random number
t
  dt  k t  
ri ,1  e 0
Ith particle
first random number
We have a problem with this integral!
We solve this by introducing a new, fictitious scattering
process which does not change energy or momentum:
ss (k )  S ( E ) (x  x) (k   k )
Generation of Random Flight Times
t
ri ,1 
(k )   i (k )
i
  dt  k t  
e 0
The sum runs over all the real
scattering processes. To this we
add the fictitious self-scattering
which is chosen to have a nice
property:new  0
ss (k )  0 
 i (k )
real
scatterers
Self-Scattering
• The use of the full integral form of the free-flight probability
density function is tedious (unless k is invariant during the
free flight).
• The introduction of self-scattering (Rees, J. Phys. Chem.
Solids 30, 643, 1969) simplifies the procedure considerably.
• The properties of the self-scattering mechanism are that it
does not change either the energy or the momentum of the
particle.
• The self-scattering rate adjusts itself in time so that the total
scattering rate is constant. Under these circumstances, one
has that:
t
  k t   self k t 
P t dt  e
  dt 
0
dt  e t dt
Self-Scattering
• Random flight times tr may be generated from P(t) above using
the direct method to get:
r e
 t r
1
1
tr   ln1  r    lnr 


where r is a uniform random between 0 and 1 (and therefore r
and 1-r are the same).
•  must be chosen (a priori) such that > (k(t)) during the entire
flight.
• Choosing a new tr after every collision generates a random walk
in k-space over which statistics concerning the occupancy of the
various states k are collected.
Free-Flight Scatter Sequence for
Ensemble Monte Carlo
Particle time
scale
n 1
2
3
4
5
6
  







 


 
 
 


 
 



 


N
ti ,1
ti ,1



 





 

  
However, we need a
second time scale,
which provides the
times at which the
ensemble is
“stopped” and
averages are
computed.
 = collisions

dte=dtau
no
yes
dte ≥ t?
dt2 = dte
dt2 = t
Call drift(dt2)
yes
Free-Flight
Scatter
Sequence
dte ≥ t?
dte2 = dte
Call scatter_carrier()
Generate free-flight dt3
dtp=t-dte2
no
dt2 = dtp
yes
dt3 ≤ dtp?
dt2 = dt3
Call drift(dt2)
dte2=dte2+dt3
dte=dte2
R. W. Hockney and J. W.
Eastwood, Computer Simulation
Using Particles, 1983.
yes
dte < t ?
no
dte=dte-t
dtau=dte
Choice of Scattering Event Terminating
Free Flight
o At the end of the free flight ti, the type of scattering which ends the
flight (either real or self-scattering) must be chosen according to the
relative probabilities for each mechanism.
o Assume that the total scattering rate for each scattering mechanism
is a function only of the energy, E, of the particle at the end of the
free flight
  self E   i E   ac E   pop E   
where the rates due to the real scattering mechanisms are typically
stored in a lookup table.
o A histogram is formed of the scattering rates, and a random number
r is used as a pointer to select the right mechanism. This is
schematically shown on the next slide.
Choice of Scattering Event Terminating
Free-Flight
We can make a table of the scattering processes at
the energy of the particle at the scattering time:
Self
0
5
1  2  3  4  5
4
1  2  3  4
3
1  2  3
2
1
1  2
1 E ti 
r0
Selection process
for scattering
Look-up table of scattering rates:
Store the total
scattering rates
in a table for a
grid in energy


4E
3E
2E
E
0
1
2
3
………
Choice of the Final State After Scattering
 Using a random number and probability
distribution function
4
6
x 10
Arbitrary Units
5
4
3
2
1
0
0
0.5
1
1.5
2
2.5
3
3.5
Polar angle
 Using analytical expressions (slides that follow)
1. Isotropic scattering processes
cos   1  2r ,   2 r
2. Anisotropic scattering processes (Coulomb, POP)
kz
k
0
Step 1:
Determine 0 and 0
ky
0
kx
kz’ k’≠k for
inelastic
kz’
Step 2:
Assume
rotated
coordinate
system
k
k

ky’
ky’
kx’
kx’
Step 3:
k’

ky’
ky’
kx’

kx’
Step 3:
perform scattering
1     1  2 
cos  
r

cos   1 
, 
2r
1  4k L (1  r )
2 2
D
=2r for both
Step 4:
kxp = k’sin cos, kyp = k’sin*sin, kzp = k’cos
Return back to the original coordinate system:
kx = kxpcos0cos0-kypsin0+kzpcos0sin0
ky = kxpsin0cos0+kypcos0+kzpsin0sin0
kz = -kxpsin0+kzpcos0

2 Ek  Ek  0 
Ek  Ek  0
Coulomb

2
POP
Representative Simulation Results From
Bulk Simulations - GaAs
Conduction bands
L-valley [111]
X-valley [100]
-valley
k-vector
Define scattering mechanisms for each valley
-valley table
-Mechanism1
-Mechanism2
-…
-MechanismN
L-valley table
-Mechanism1
-Mechanism2
-…
-MechanismNL
X-valley table
-Mechanism1
-Mechanism2
-…
-MechanismNx
Valence bands
Call specified
scattering mechanisms subroutines
Simulation Results Obtained by
D. Vasileska’s Monte Carlo Code.
Renormalize scattering tables
parameters initialization
readin()
15
10
intervalley gamma to X
14
10
scattering table construction
sc_table()
Scattering Rate [1/s]
1
0.9
0.8
Cumulative rate [1/s]
carriers initialization
init()
histograms calculation
histograms()
0.7
0.6
0.5
polar optical phonons
13
10
12
10
11
10
intervalley gamma to L
acoustic
0.4
0.3
10
10
0
0.1
0.2
0.3
0.4
0.2
0.5
0.6
0.7
0.8
0.9
1
Energy [eV]
0.1
Free-Flight-Scatter
free_flight_scatter()
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Energy [eV]
histograms calculation
histograms()
Optional
40
35
90
Initial Distribution of the
wavevector along the y-axis
that is created with the
subroutine init()
80
70
Arbitrary Units
Arbitrary Units
30
write data
write()
Initial Energy Distribution created
with the subroutine init()
25
20
15
60
50
40
30
10
20
t  t  t
no
5
Time t exceeds maximum
simulation time tmax
?
yes
0
-6
10
-4
-2
0
2
Wavevector ky [1/m]
4
6
0
0
0.05
0.1
0.15
0.2
8
x 10
Energy [eV]
0.25
0.3
0.35
parameters initialization
readin()
15
10
intervalley gamma to X
14
10
scattering table construction
sc_table()
Scattering Rate [1/s]
1
0.9
0.8
Cumulative rate [1/s]
carriers initialization
init()
histograms calculation
histograms()
0.7
0.6
0.5
polar optical phonons
13
10
12
10
11
10
intervalley gamma to L
acoustic
0.4
0.3
10
10
0
0.1
0.2
0.3
0.4
0.2
0.5
0.6
0.7
0.8
0.9
1
Energy [eV]
0.1
Free-Flight-Scatter
free_flight_scatter()
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Energy [eV]
histograms calculation
histograms()
Optional
40
35
90
Initial Distribution of the
wavevector along the y-axis
that is created with the
subroutine init()
80
70
Arbitrary Units
Arbitrary Units
30
write data
write()
Initial Energy Distribution created
with the subroutine init()
25
20
15
60
50
40
30
10
20
t  t  t
no
5
Time t exceeds maximum
simulation time tmax
?
yes
0
-6
10
-4
-2
0
2
Wavevector ky [1/m]
4
6
0
0
0.05
0.1
0.15
0.2
8
x 10
Energy [eV]
0.25
0.3
0.35
Transient Data
5
4
x 10
3.5
velocity [m/s]
3
2.5
2
1.5
1
0.5
0
0
1
time [s]
2
-11
x 10
Conduction bands
Steady-State
Results
L-valley [111]
X-valley [100]
-valley
k-vector
Gunn Effect
Valence bands
5
2
x 10
10000
Conduction Band Valley Occupancy
1.8
Drift Velocity [m/s]
1.6
1.4
1.2
1
0.8
0.6
0.4
gamma valley occupancy
8000
7000
6000
5000
4000
3000
L valley occupancy
2000
X valley occupancy
1000
0.2
0
9000
0
1
2
3
4
Electric Field [V/m]
5
6
7
5
x 10
0
0
1
2
3
4
Electric Field [V/m]
5
6
7
5
x 10
Particle-Based Device Simulations
In a particle-based device simulation
approach the Poisson equation is
decoupled from the BTE over a short time
period dt smaller than the dielectric
relaxation time
Poisson and BTE are solved in a time-marching
manner
During each time step dt the electric field is
assumed to be constant (kept frozen)
Particle-Mesh Coupling
The particle-mesh coupling scheme consists of the
following steps:
- Assign charge to the Poisson solver mesh
- Solve Poisson’s equation for V(r)
- Calculate the force and interpolate it to the
particle locations
- Solve the equations of motion:
dr 1
 kE kt  ;
dt 
dk qEr 

dt

Laux, S.E., On particle-mesh coupling in Monte Carlo semiconductor device simulation,
Computer-Aided Design of Integrated Circuits and Systems, IEEE Transactions on,
Volume 15, Issue 10, Oct 1996 Page(s):1266 - 1277
Assign Charge to the Poisson Mesh
1. Nearest grid point scheme
2. Nearest element cell scheme
3. Cloud in cell scheme
Force interpolation
The SAME METHOD that is used for the
charge assignment has to be used for the
FORCE INTERPOLATION:
1 Vp 1  Vp Vp  Vp 1 

F  ri   q W  ri  rp E  rp  , Ep   

e
w
2  x p
x p 
p
x wp
xp-1
x ep
xp
xp+1
Treatment of the Contacts
From the aspect of device physics, one can distinguish
between the following types of contacts:
(1) Contacts, which allow a current flow in and out of the
device
- Ohmic contacts: purely voltage or purely current
controlled
- Schottky contacts
(2) Contacts where only voltages can be applied
Calculation of the Current
The current in steady-state conditions is
calculated in two ways:
 By counting the total number of particles that
enter/exit particular contact
 By using the Ramo-Shockley theorem
according to which, in the channel, the current
is calculated using
e N
I
v x (i ),

dL i 1
Current Calculated by Counting the Net Number of
Particles Exiting/Entering a Contact
Source
Gate
Drain
Mesh node
Electron
Dopant
Electrons that naturally came out in time interval dt (N1)
Electrons that were deleted (N2)
Electrons that were injected (N3)
dq = q(N1+N2-N3), q(t+dt)=q(t) + dq, current equals the slope of q(t) vs. t
Device Simulation Results for MOSFETs:
Current Conservation
VG=1.4 V, VD=1 V
0.6
Drain
contact
source
contact
drain contact
Source
contact
[mA/mm]
(a)
6000
5000
D
4000
3000
Current I
net # of electrons
exiting/entering contact
7000
2000
1000
W = 0.5 mm
G
0
1
1.5
2
2.5
3
time [ps]
Cumulative net number of particles
Entering/exiting a contact for a 50 nm
Channel length device
X. He, MS, ASU, 2000.
(b)
0.5
0.4
0.3
0.2
0.1
0
0
50
100
150
Distance [nm]
Current calculated using Ramo-Schockley
formula
e N
I
v x (i ),

dL i 1
Simulation Results for MOSFETs: Velocity
and Enery Along the Channel
2.5x10
7
2x10
7
1.5x10
7
1x10
7
5x10
6
Average Kinetic Energy
Along the Channel
Average energy [eV]
Drift velocity [cm/s]
Mean Drift Velocity
Along the Channel
(a)
0
0
50
100
0.6
0.4
0.3
0.2
0.1
0
150
Distance [nm]
(b)
0.5
0
50
100
150
Distance [nm]
VD = 1 V, VG = 1.2 V
Velocity overshoot effect observed
throughout the whole channel length of
the device – non-stationary transport.
For the bias conditions used average
electron energy is smaller that 0.6 eV
which justifies the use of non-parabolic
band model.
Drain current [mA/mm]
Simulation Results for MOSFETs:
IV Characteristics
0.5
The differences between
the Monte Carlo and the
Silvaco simulations are
due to the following
reasons:
VG = 1.4 V
2D MCPS
0.4
VG = 1.0 V
0.3
• Different transport
models used (nonstationary transport is
taking place in this device
structure).
VG = 1.4 V
1.2 V
0.2
1.0 V
0.8 V
0.1
Silvaco simulations
0
0
0.2
0.4
0.6
0.8
Drain voltage V [V]
D
X. He, MS, ASU, 2000.
1
• Surface-roughness and
Coulomb scattering are
not included in the
theoretical model used in
the 2D-MCPS.
Simulation Results For SOI MESFET
Devices – Where are the Carriers?
Lg =60, 100nm
SOI
MOSFET
Nd = 1019 Na =3-10x 1015 Nd = 1019
Lg =60, 100nm
Nd = 1019 Nd = 3-10x1015
Oxide Layer
Oxide Layer
Si Substrate
Si Substrate
SOI
MESFET
Nd =1019
Fig. 3.1 Schematic cross-sections of a) the SOI MOSFET and b) the SOI
MESFET devices that have been simulated.
Applications:
Micropower circuits based
on weakly inverted
MOSFETs
Digital Watch
Pacemaker
Implantable
cochlea and retina
Low-power RF electronics.
3.2. The electron
distributions
in c)Dev.
the 60Lett.,
nm SOI8171
MOSFET
and d) the 60 nm
T.J.Fig.
Thornton,
IEEE
Electron
(1985).
Proper Modeling of SOI MESFET Device
Gate current calculation:
 WKB Approximation
 Transfer Matrix Approach
for piece-wise linear
potentials
Interface-Roughness:
 K-space treatment
 Real-space treatment
Goodnick et al., Phys. Rev. B 32, 8171
(1985)
Output Characteristics and Cut-off
Frequency of a Si MESFET Device
400
12
110
10
-200
10
10
-400
Vgs = 0.3V
d
I [µA/fµTm]
Cutoff Frequency
[Hz]
10
200
Vgs = 0.4V
-600
9
10
10
-1000
-0.2
0.1
V
= 0.6V
Lg =25nm simulated results
gs
Lg =50nm simulated results
Lg = 90nm simulated results
0
0.2Experimental
0.4 results
0.6
0.8
Lg = 0.6um
Lg = 50nm Projected Experimental results
V
7
10
=0.5V
gs
-800
8
V
1
ds
[Volts]
100
10
Drain Current I [ µA/µm]
d
Tarik Khan, PhD, ASU, 2004.
1
1000
Output Characteristics and Cut-off
Frequency of a Si MESFET Device
400
12
110
10
-200
10
10
-400
Vgs = 0.3V
d
I [µA/fµTm]
Cutoff Frequency
[Hz]
10
200
Vgs = 0.4V
-600
9
10
10
-1000
-0.2
0.1
V
= 0.6V
Lg =25nm simulated results
gs
Lg =50nm simulated results
Lg = 90nm simulated results
0
0.2Experimental
0.4 results
0.6
0.8
Lg = 0.6um
Lg = 50nm Projected Experimental results
V
7
10
=0.5V
gs
-800
8
V
1
ds
[Volts]
100
10
Drain Current I [ µA/µm]
d
Tarik Khan, PhD, ASU, 2004.
1
1000
Modeling of SOI Devices
When modeling SOI devices lattice
heating effects has to be accounted for
In what follows we discuss the following:
 Comparison of the Monte Carlo, Hydrodynamic
and Drift-Diffusion results of Fully-Depleted SOI
Device Structures*
 Impact of self-heating effects on the operation
of the same generations of Fully-Depleted SOI
Devices
*D. Vasileska. K. Raleva and S. M. Goodnick, IEEE Trans. Electron Dev., in press.
FD-SOI Devices:
Monte Carlo vs. Hydrodynamic vs. Drift-Diffusion
1.6
1.4
Source
tsi
Drain
LS
Lgate
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%
Source/drain doping = 1020 cm-3 and 1019 cm-3 (series resistance (SR) case)
Channel doping = 1E18 cm-3
Overshoot= (IDHD-IDDD)/IDDD (%) at on-state
Drain Current
Drain [mA/um]
Drain Current
Current
[mA/um][mA/um]
tox
Gate oxide
2.5
1.2
1
2
3
0.8
2.5
1.5
0.6
90 nm
2
0.4
1
DD SR
HD SR
Monte Carlo
0.2
1.5
0.5
0
1
0
0.4
0.6
0.8
Drain Voltage [V]
0
0.50
0
0.2
0
0.2
0.2
DD SR
0.8 HD SR1
Drain Voltage [V] Monte Carlo
0.4
0.6
0.4
0.6
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
DD SR
1HD SR1.2
Monte Carlo
0.8
1.4
1.2
1
FD-SOI Devices:
Monte Carlo vs. Hydrodynamic vs. Drift-Diffusion
1.6
1.4
Source
tsi
Drain
LS
Lgate
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%
Source/drain doping = 1020 cm-3 and 1019 cm-3 (series resistance (SR) case)
Channel doping = 1E18 cm-3
Overshoot= (IDHD-IDDD)/IDDD (%) at on-state
Drain Current
Drain [mA/um]
Drain Current
Current
[mA/um][mA/um]
tox
Gate oxide
2.5
1.2
1
2
3
0.8
2.5
1.5
0.6
2
0.4
1
DD SR
HD SR
Monte Carlo
0.2
1.5
0.5
0
1
25 nm
0
0.4
0.6
0.8
Drain Voltage [V]
0
0.50
0
0.2
0
0.2
0.2
DD SR
0.8 HD SR1
Drain Voltage [V] Monte Carlo
0.4
0.6
0.4
0.6
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
DD SR
1HD SR1.2
Monte Carlo
0.8
1.4
1.2
1
FD-SOI Devices:
Monte Carlo vs. Hydrodynamic vs. Drift-Diffusion
1.6
1.4
Source
tsi
Drain
LS
Lgate
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%
Source/drain doping = 1020 cm-3 and 1019 cm-3 (series resistance (SR) case)
Channel doping = 1E18 cm-3
Overshoot= (IDHD-IDDD)/IDDD (%) at on-state
Drain Current
Drain [mA/um]
Drain Current
Current
[mA/um][mA/um]
tox
Gate oxide
2.5
1.2
1
2
3
0.8
2.5
1.5
0.6
2
0.4
1
DD SR
HD SR
Monte Carlo
0.2
1.5
0.5
0
1
0
0
0.50
0
0
0.2
0.4
0.2
0.2
0.6
0.8
Drain Voltage [V] 14
nm
DD SR
0.8 HD SR1
Drain Voltage [V] Monte Carlo
0.4
0.6
0.4
0.6
Drain Voltage [V]
Silvaco ATLAS simulations performed by Prof. Vasileska.
DD SR
1HD SR1.2
Monte Carlo
0.8
1.4
1.2
1
FD-SOI Devices:
Why Self-Heating Effect is Important?
1. Alternative materials (SiGe)
2. Alternative device designs (FD SOI, DG,
TG, MG, Fin-FET transistors
FD-SOI Devices:
Why Self-Heating Effect is Important?
L~
300nm
dS
A. Majumdar, “Microscale Heat Conduction
in Dielectric Thin Films,” Journal of Heat
Transfer, Vol. 115, pp. 7-16, 1993.
experimental data
full lines: BTE predictions
dashed lines: empirical model
thin lines: Sondheimer
60
100nm
40
50nm
30nm
20
20nm
300
400
500
Temperature (K)
 /2
600



 a  2 z 
a
 ( z)   0 (T )  sin  1  exp  
cosh


d
2

(
T
)cos

2

(
T
)cos





0

3
 (T )  0 (300 / T )
135
 0 (T ) 
a  bT  cT 2
W/m/K
20
17.5
15
Si
BOX
10nm
300K
400K
600K
12.5
10
7.5
5
0
2
4
6
8
10
Distance from Si/gate oxide interface (nm)
Thermal conductivity (W/m-K)
80
Thermal conductivity (W/m-K)
Thermal conductivity (W/m/K)
Conductivity of Thin Silicon Films –
Vasileska Empirical Formula
80
70
60
50
40
30
20
10
20nm
30nm
50nm
100nm
0
0
20
40
60
80
100
Distance from Si/gate oxide interface (nm)
Particle-Based Device Simulator That
Includes Heating
Define device structure
Generate phonon temperature
dependent scattering tables
Initial potential, fields, positions and
velocities of carriers
Average and smooth: electron
density, drift velocity and electron
energy at each mesh point
end of MCPS
phase?
t=0
t = t + t
Acoustic and Optical Phonon
Energy Balance Equations Solver
Transport Kernel (MC phase)
no
t = n t?
yes
Field Kernel (Poisson Solver)
end of
simulation?
yes
end
Heating vs. Different Technology
Generation
Acoustic Phonon Temperature Profiles
T=300K on gate
T=400K on gate
25 nm FD SOI nMOSFET (Vgs=Vds=1.2V)
25 nm FD SOI nMOSFET (Vgs=Vds=1.2V)
3
source contact
drain contact
3
500
source region
20
12
35
21
300
35
21
300
80 nm FD SOI nMOSFET (Vgs=Vds=1.5V)
10
20
30
40
50
60
70
45 nm FD SOI nMOSFET (Vgs=Vds=1.2V)
20 50
100
60
40
80150
100
90
60 nm
nm FD
FD SOI
SOI nMOSFET
nMOSFET (Vgs=Vds=1.5V)
(Vgs=Vds=1.2V)
3
3
200120
4
50
100
150
200
300
70
600
500
400
20 50
3
500
500
21
400
15
400
300
39
250
100 nm FD SOI nMOSFET (Vgs=Vds=1.5V)
120 nm 60
FD SOI nMOSFET
40
80
100(Vgs=Vds=1.8V)
120
140
20
8020
nm FD SOI
30nMOSFET
40 (Vgs=Vds=1.5V)
50
60
45 nm FD SOI nMOSFET (Vgs=Vds=1.2V)
10
100
60
40
80150
100
300
200120
90
60 nm FD SOI nMOSFET (Vgs=Vds=1.5V)
(Vgs=Vds=1.2V)
21
15
39
400
drain region
13
300
6003
3
500
20
12
400
27
3
500
drain contact
8
400
8
13
3
3
600
source contact
160
180
27
300
3
4
23
28
600
500
23
500
28
400
43
52
52
300
43
300
600
500
400
400
50
100
150
200
100 nm FD SOI nMOSFET (Vgs=Vds=1.5V)
120 nm 60
FD SOI nMOSFET
40
80
100(Vgs=Vds=1.8V)
120
140
20
250
160
180
400
50
0
100
120
150
200
240
250
300
360
140 nm FD SOI nMOSFET (Vgs=Vds=1.8V)
4
4
600
33
33
400
60
60
50
100
150
200
250
300
350
400
4
600
41
41
400
76
100
200
300
x (nm)
400
500
100100 150 150200
250
200
140 nm FD SOI nMOSFET (Vgs=Vds=1.8V)
300
250
350
300
600
700
600
500
500
400
400
300
600
400
50
100
150
200
250
300
350
400
180 nm FD SOI nMOSFET (Vgs=Vds=1.8V)
180 nm FD SOI nMOSFET (Vgs=Vds=1.8V)
4
76
5050
300
600
400
100
200
300
x (nm)
400
500
Higher Order Effects Inclusion in ParticleBased Simulators
 Degeneracy – Pauli Exclusion Principle
 Short-Range Coulomb Interactions
 Fast Multipole Method (FMM)
V. Rokhlin and L. Greengard, J. Comp. Phys., 73, pp. 325-348
(1987).
 Corrected Coulomb Approach
W. J. Gross, D. Vasileska and D. K. Ferry, IEEE Electron Device
Lett. 20, No. 9, pp. 463-465 (1999).
 P3M Method
Hockney and Eastwood, Computer Simulation Using Particles.
MOTIVATION
Potential, Courtesy of Dragica Vasileska, 3D-DD Simulation, 1994.
MOTIVATION
150
Width [nm]
140
130
120
110
100
60
80
100
120
140
Length [nm]
Current Stream Lines, Courtesy of Dragica Vasileska, 3D-DD Simulation, 1994.
The ASU Particle-Based Device
Simulator
Short-Range Interactions (1) Corrected
Coulomb Approach
and Discrete/Unintentional
Discrete Impurities
(2) P3M Algorithm
Dopants
(Short-Range Interaction) (3) Fast Multipole
Method (FMM)
Effective
Potential
Quantum
Mechanical
(Space Quantization)
Size-quantization
Effects
(1)
(2)
3DMCDS
Long-range Interactions
Efficient 3D Poisson
(3D Poisson
Equation Solvers
Equation Solver)
Ferry’s Effective
Potential Method
Quantum Field
Approach
3D Monte
CarloEquations
Boltzmann
Transport
Transport
KernelTransport Kernel)
(Particle-Based
Monte Carlo
Statistical Enhancement:
Event Biasing Scheme
Significant Data Obtained
Between 1998 and 2002
MOSFETs - Standard Characteristics
2x10
300
D
200
D
G
V =0.5 [V], V =1.0 [V]
D
250
7
(a)
V =0.5 [V], V =1.0 [V]
V =1.0 [V], V =1.0 [V]
Drift velocity [cm/s]
Electron energy [meV]
350
G
L = 80 nm
G
150
100
1.5x10
G
V =1.0 [V], V =1.0 [V]
7
D
G
L =80 nm
1x10
7
5x10
6
G
50
0
40
60
80
100
Distance [nm]
120
140
0
40
60
80
100
120
140
Distance [nm]
 The average energy of the carriers increases when going from the
source to the drain end of the channel. Most of the phonon scattering
events occur at the first half of the channel.
 Velocity overshoot occurs near the drain end of the channel. The
sharp velocity drop is due to e-e and e-i interactions coming from the
drain.
W. J. Gross, D. Vasileska and D. K. Ferry, "3D Simulations of Ultra-Small MOSFETs with Real-Space
Treatment of the Electron-Electron and Electron-Ion Interactions," VLSI Design, Vol. 10, pp. 437-452 (2000).
2.5x10
7
2x10
7
1.5x10
7
with e-e and e-i
mesh force only
V =V =1.0 V
D
1x10
7
5x10
6
Electron energy [meV]
Drift velocity [cm/s]
MOSFETs - Role of the E-E and E-I
G
0
-5x10
6
-1x10
7
0
source
channel
40
80
drain
120
400
with e-e and e-i
mesh force only
300
V =1 V, V =1 V
D
200
100
channel drain
0
100 110 120 130 140 150 160 170 180
160
Length [nm]
Length [nm]
800
800
mesh force
only
600
500
400
300
200
100
0
0.12
with e-e and e-i
700
Energy [meV]
Energy [meV]
700
Individual electron
trajectories over
time
G
600
500
400
300
200
100
0.13
0.14
0.15
0.16
Length [nm]
0.17
0.18
0
0.12
0.13
0.14
0.15
0.16
Length [nm]
0.17
0.18
MOSFETs - Role of the E-E and E-I
Mesh force only
With e-e and e-i
V =0.5 V, V =0.8 V
10
-2
10
-3
G
D
Source
Channel
Drain
0
50
100 150 200 250 300 350 400
Energy [meV]
G
Electron distribution
(arb. units)
Electron distribution
(arb. units)
V =0.5 V, V =0.8 V
10
-2
10
-3
D
source
channel
drain
0
50
100 150 200 250 300 350 400
Energy [meV]
Short-range e-e and e-i interactions push some
of the electrons towards higher energies
D. Vasileska, W. J. Gross, and D. K. Ferry, "Monte-Carlo particle-based simulations of deep-submicron nMOSFETs with real-space treatment of electron-electron and electron-impurity interactions," Superlattices
and Microstructures, Vol. 27, No. 2/3, pp. 147-157 (2000).
Degradation of Output Characteristics
 The short range e -e and e -i interactions have significant influence on
the device output characteristics.
 There is almost a factor of two decrease in current when these two interaction terms are considered.
60
[mA]
with corrected Coulomb
mesh force only
increasing V
D
70
80
G
50
Drain current I
Drain current I
D
[mA]
80
LG = 50 nm, WG = 35 nm, NA = 5x1018 cm-3
Tox = 2 nm, VG = 11.6 V (0.2 V)
LG = 35 nm, WG = 35 nm, NA = 5x1018 cm-3,
Tox = 2 nm, VG = 11.6 V (0.2 V)
40
30
20
10
0
0.0
0.2
0.4
0.6
0.8
Drain voltage V
D
[V]
1.0
1.2
with corrected Coulomb
70
mesh force only
60
increasing V
50
G
40
30
20
10
0
0
0.2
0.4
0.6
Drain voltage V
0.8
D
1
1.2
[V]
W. J. Gross, D. Vasileska and D. K. Ferry, "Ultra-small MOSFETs: The importance of the full Coulomb
interaction on device characteristics," IEEE Trans. Electron Devices, Vol. 47, No. 10, pp. 1831-1837 (2000).
Mizuno result:
(60% of the fluctuations)
Stolk et al. result:
(100% of the fluctuations)
DepthDistribution
of the charges
Fluctuations in the
surface potential
Fluctuations in the
electric field
MOSFETs - Discrete Impurity Effects
Approach 1 [1]:
4 q 3  T
Si B ox
sVth 
ox
2
Approach 2 [2]:
4 4q 3  
Si B
sVth 

3
k T N 
;  B  B ln A 
q
Leff Weff
 ni 
4N
A
T  4 NA
k BT / q
 ox 
 4q Si  B N A ox  Leff Weff
[1] T. Mizuno, J. Okamura, and A. Toriumi, IEEE Trans. Electron Dev. 41, 2216 (1994).
[2] P. A. Stolk, F. P. Widdershoven, and D. B. Klaassen, IEEE Trans. Electron Dev. 45, 1960
(1998).
40
60
100
s => approach 1
Vth
s => approach 2
Vth
s => our simulation results
5
s => approach 1
Vth
s => approach 2
Vth
s => our simulation results
0
18
1x10
3x10
10
40
Vth
40
30
s
15
Vth
60
Vth
[mV]
20
Vth
25
s => approach 1
Vth
s => approach 2
Vth
s => our simulation results
50
[mV]
80
30
s
s
Vth
[mV]
35
20
20
Vth
18
5x10
18
Doping density N
7x10
-3
A
[cm ]
18
0
0
1
2
3
4
Oxide thickness T
ox
[nm]
5
10
20
40
60
80
100
120
Device width [nm]
140
Depth Correlation of sV
 To understand the role that the position
of the impurity atoms plays on the
threshold voltage fluctuations, statistical
ensembles of 5 devices from the low-end,
center and the high-end of the
distribution were considered.
200
180
160
140
120
100
80
60
40
20
0
1.4
(a)
L =50 nm, W =35 nm
G
G
1.3
N =5x10
18
-3
cm , T =3 nm
ox
A
1.2
high-end
1.1
low-end
center
1
0.9
160
220
200
180
300
280
260
240
atoms
dopant
channel in
of atoms
Numberof
Number
the channel
1
5 samples of average
(b)
5 samples at
maximum
270
260
250
240
230
220
210
200
190
180
170
VVTT correlation
correlation
5 samples at
minimum
160
Number
Numberof
of devices
Devices
 Significant correlation was observed
between the threshold voltage and the
number of atoms that fall within the first 15
nm depth of the channel.
Threshold voltage
voltage [V][V]
Threshold
T
L =50 nm, W =35 nm, T =3 nm
G
G
0.8
N =5x10
18
A
cm
-3
0.6
Moving slab
0.4
ND
range
ND
0.2
depth
Number of Atoms in Channel
Number of atoms in the channel
ox
0
0
5
10
15
Depth [nm]
Depth
[nm]
20
25
Fluctuations in High-Field
Characteristics
1.5 10
center
10
low-end
5
high-end
0
160
180
200
220
240
260
Number of channel dopant atoms
Number of atoms in the channel
(c)
7
6
280
1
0.8
center
5 10
15
(a)
low-end
1 10
VG = 1.5 V, VD = 1 V
7
Correlation
Correlation
[cm/s]
Drift
Drift velocity
velocity [cm/s]
 Significant correlation was observed
between the drift velocity (saturation
current) and the number of atoms that fall
within the first 10 nm depth of the
channel.
(b)
current
Drain
Drain
current[mA]
[mA]
 Impurity distribution in the channel also
affects the carrier mobility and saturation
current of the device.
20
VG=1.5 V, VD=1 V
high-end
LG=50 nm, WG=35 nm
18
-3
NA=5x10 cm
0
160
180
200
220
240
260
280
Number of channel dopant atoms
Number
of atoms in the channel
average velocity
correlation
drain current
correlation
0.6
0.4
LG=50 nm
0.2
Tox=3 nm
WG=35 nm
18
-3
NA=5x10 cm
0
0
5
10
15
Depth [nm]
Depth
[nm]
20
25
Current Issues in Novel
Devices – Unintentional Dopants
THE EXPERIMENT …
Results for SOI Device
Size Quantization Effect (Effective Potential):
0.6
Threshold Voltage [V]
Experimental
0.5
Simulation
0.4
0.3
0.2
0.1
0
-0.1
-0.2
2
4
6
8
10
12
14
16
Channel Width [nm]
S. S. Ahmed and D. Vasileska, “Threshold voltage shifts in narrow-width SOI devices due to
quantum mechanical size-quantization effect”, Physica E, Vol. 19, pp. 48-52 (2003).
Results for SOI Device
0.18
Velocity
Energy
Average Velocity [m/s]
110000
0.16
0.14
90000
0.12
70000
0.1
50000
0.08
0.06
Dip due to the presence of
the impurity. This affects the
transport of the carriers.
30000
0.04
10000
0.02
-10000
0
0
20
40
60
Distance Along the Channel [nm]
Due to the unintentional dopant both the electrostatics
and the transport are affected.
80
Average Kinetic Energy [eV]
130000
Results for SOI Device
Unintentional Dopant:
D. Vasileska and S. S. Ahmed, “Narrow-Width SOI Devices: The Role of Quantum Mechanical
Size Quantization Effect and the Unintentional Doping on the Device Operation”, IEEE
Transactions on Electron Devices, Volume 52, Issue 2, Feb. 2005 Page(s):227 – 236.
33.01%
27.18%
11.11%
6
5
4
47.62%
42.85%
26.19%
Drain
7
3
2
32.54%
26.98%
11.90%
1
0
0
10
20
30
40
Distance Along the Channel [nm]
50
1
42.06%
26.19%
19.46%
47.62%
42.85%
26.19%
34.13%
16.66%
9.93%
2
3
4
Drain
8
0
Source
9
Distance Along the Depth [nm] .
10
Source
Distance Along the Width [nm] .
Results for SOI Device
5
6
7
0
10
20
30
40
Distance Along the Channel [nm]
Channel Width = 10 nm
VG = 1.0 V
VD = 0.1 V
50
2
86.30%
1
96.09%
87.39%
88.26%
69.57%
0
0
10
20
30
40
Distance Along the Channel [nm]
50
1
88.48%
88.26%
76.09%
96.76%
96.09%
88.26%
81.09%
79.78%
59.78%
2
3
4
Drain
96.76%
67.39%
0
Source
3
86.96%
Drain
86.52%
4
Distance Along the Depth [nm]
5
Source
Distance Along the Width [nm]
Results for SOI Device
5
6
7
0
10
20
30
40
Distance Along the Channel [nm]
Channel Width = 5 nm
VG = 1.0 V
VD = 0.1 V
50
Results for SOI Device
60%
Impurity position varying along the
center of the channel
Current Reduction
50%
V G = 1.0 V
40%
V D = 0.2 V
30%
20%
Source end
Drain end
10%
0%
0
10
20
30
40
50
Distance Along the Channel [nm]
Impurity located at the very source-end, due to the availability of Increasing
number of electrons screening the impurity ion, has reduced impact on the
overall drain current.
Results for SOI Device
Threshold Voltage [V]
0.6
Experimental
Simulation (QM)
Discrete single dopants
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
2
4
6
8
10
12
14
16
Channel Width [nm]
D. Vasileska and S. S. Ahmed, IEEE Transactions on Electron Devices, Volume 52, Issue 2, Feb. 2005
Page(s):227 – 236.
S. Ahmed, C. Ringhofer and D. Vasileska, Nanotechnology, IEEE Transactions on, Volume 4, Issue 4, July 2005
Page(s):465 – 471.
D. Vasileska, H. R. Khan and S. S. Ahmed, International Journal of Nanoscience, Invited Review Paper, 2005.
Results for SOI Device
Electron-Electron Interactions:
1.E+02
PM
FMM
2.5E+05
2.0E+05
V G = 1.0 V
V D = 0.3 V
1.5E+05
1.0E+05
5.0E+04
Source
Drain
Distribution Function [a.u.]
Electron Velocity [m/s]
3.0E+05
PM
FMM
1.E+01
V G = 1.0 V
V D = 0.3 V
1.E+00
1.E-01
1.E-02
1.E-03
1.E-04
0.0E+00
0
20
40
60
80
100
Distance Along the Channel [nm]
0
0.2
0.4
0.6
0.8
1
Electron Kinetic Energy [eV]
D. Vasileska and S. S. Ahmed, IEEE Transactions on Electron Devices, Volume 52, Issue 2, Feb. 2005
Page(s):227 – 236.
S. Ahmed, C. Ringhofer and D. Vasileska, Nanotechnology, IEEE Transactions on, Volume 4, Issue 4, July 2005
Page(s):465 – 471.
D. Vasileska, H. R. Khan and S. S. Ahmed, International Journal of Nanoscience, Invited Review Paper, 2005.
Summary
 Particle-based device simulations are the most desired
tool when modeling transport in devices in which velocity
overshoot (non-stationary transport) exists
 Particle-based device simulators are rather suitable for
modeling ballistic transport in nano-devices
 It is rather easy to include short-range electron-electron
and electron-ion interactions in particle-based device
simulators via a real-space molecular dynamics routine
 Quantum-mechanical effects (size quantization and
density of states modifications) can be incorporated in
the model quite easily with the assumption of adiabatic
approximation and solution of the 1D or 2D Schrodinger
equation in slices along the channel section of the device