Download Fortran Code for Calculating Electric Field and Charge Density in a Semiconductor Device and more Study notes Theatre in PDF only on Docsity!
6. Particle-Based Device Simulation Methods
Figure 3 illustrated various levels of approximation in describing charge transport
within a hierarchical structure ranging from the exact quantum mechanical solution of the
n -particle problem at the bottom, to analytical 1D phenomenological modeling used in
circuit simulation at the top. The exact quantum mechanical solution of even a few
particle system is a challenging computational task, and clearly impossible for a
semiconductor device with typical free-carrier electron densities that are on the order of
17
cm
or more. Hence, simplifying approximations are necessary.
For conventional semiconductor devices, such as bipolar junction transistors
(BJTs) and field effect transistors (FETs), the device behavior has been adequately
described within the semi-classical model of charge transport, since the characteristic
dimensions are typically at length scales much larger than those over which quantum
mechanical phase coherence is maintained. Hence a particle-based description is
adequate as described within the Boltzmann equation framework, and approximations
thereof. As device dimensions continue to shrink, the channel lengths are now
approaching the characteristic wavelength of particles (the de Broglie wavelength at the
Fermi energy, for example), and quantum effects are expected to be increasingly
important. It has in fact been well known for 30 years that quantum confinement effects
occur for electrons in the inversion layers of Si metal-oxide-semiconductor field effect
transistor (MOSFET) devices. However, at room temperature and under strong driving
fields, such quantum effects have usually been found to be second order at best in terms
of the overall device behavior. However, as discussed in Section ??, it is not clear that
this situation will persist as all spatial dimensions are reduced, and consideration of
quantum effects, such as tunneling and interference, may in fact dominate.
As mentioned above, the classical description of charge transport is given by the
BTE in the hierarchy of Figure 3. The BTE is an integral-differential kinetic equation of
motion for the probability distribution function for particles in the 6-dimensional phase
space of position and (crystal) momentum
Coll
t
f t
E f t f t
t
f t
(, ,)
( ) ( , , ) ( , , )
( , , ) 1 rk
r k
F
k r k
r k
k r k
where
f ( r , k , t )
is the one-particle distribution function. The right hand side is the rate
of change of the distribution function due to randomizing collisions, and is an integral
over the in-scattering and the out-scattering terms in momentum (wavevector) space.
Once
f ( r , k , t )
is known, physical observables, such as average velocity or current, are
found from averages of f. Equation is semi-classical in the sense that particles are
treated as having distinct position and momentum in violation of the quantum uncertainty
relations, yet their dynamics and scattering processes are treated quantum mechanically
through the electronic band structure (discussed in Section ??) and the use of time
dependent perturbation theory. Through moment expansion of the BTE, a set of
approximate partial differential equations in position space, similar to those arising in the
field of fluid dynamics, are obtained leading to the so called hydrodynamic model for
charge transport, discussed in Section ??. The simplification of the hydrodynamic model
to include just the continuity equation and the current density written in terms of the local
electric field and concentration gradients leads to the so-called drift-diffusion model, also
discussed in Section ??. Finally, the reduction of the drift-diffusion model to one
dimensional nonlinear analytical expressions allows for the development of lumped
parameter behavioral models suitable for circuit level simulation of many individual
devices as well as passive elements.
The BTE itself is an approximation to the underlying many body classical
Liouville equation, and quantum mechanically by the Liouville-von Neumann equation
of motion for the density matrix. The main approximations inherent in the BTE are the
assumption of instantaneous scattering processes in space and time, the Markov nature of
scattering processes (i.e. that they are uncorrelated with the prior scattering events), and
the neglect of multi-particle correlations (i.e. that the system may be characterized by a
single particle distribution function). In semi-classical simulation, some of these
assumptions are relaxed through the use of molecular dynamics techniques discussed in
Section ?? (in the context of device simulations). However, the inclusion of quantum
effects such as particle interference, tunneling, etc. which take one further down the
hierarchy of Figure 3 is more problematic in the semi-classical Ansatz , and is an active
area of research today as device dimensions approach the quantum regime.
without scattering after the previous collision at t = 0, and then suffer a collision in a
time interval dt around time t. The probability of scattering in the time interval dt around
t may be written as [ k ( t )] dt , where [ k ( t )] is the scattering rate of an electron or hole of
wavevector k. The scattering rate, [ k ( t )], represents the sum of the contributions from
each individual scattering mechanism, which are usually calculated using perturbation
theory, as described later. The implicit dependence of [ k ( t )] on time reflects the change
in k due to acceleration by internal and external fields. For electrons subject to time
independent electric and magnetic fields, Eq. ?? may be integrated to give the time
evolution of k between collisions as
e t
t
E v B
k k
0 ,
where E is the electric field, v is the electron velocity (given by Eq. ??, and B is the
magnetic flux density. In terms of the scattering rate, [ k ( t )], the probability that a
particle has not suffered a collision after a time t is given by
t
t d t
0
exp k
.
Thus, the probability of scattering in the time interval dt after a free flight of time t may
be written as the joint probability
P t dt t t dt dt
t
0
k exp k .
Random flight times may be generated according to the probability density P ( t ) above
using, for example, the pseudo-random number generator implicit on most modern
computers, which generate uniformly distributed random numbers in the range [0,1].
Using a direct method (see, for example [Error: Reference source not found]), random
flight times sampled from P ( t ) may be generated according to
r P t dt
r
t
0
where r is a uniformly distributed random number and t r
is the desired free flight time.
Integrating with P ( t ) given by above yields
r
t
r t d t
0
1 exp k
Since 1- r is statistically the same as r , may be simplified to
r
t
r t d t
0
ln k .
Equation is the fundamental equation used to generate the random free flight time after
each scattering event, resulting in a random walk process related to the underlying
particle distribution function. If there is no external driving field leading to a change of k
between scattering events (for example in ultrafast photoexcitation experiments with no
applied bias), the time dependence vanishes, and the integral is trivially evaluated. In the
general case where this simplification is not possible, it is expedient to introduce the so
called self-scattering method [
vi
], in which we introduce a fictitious scattering mechanism
whose scattering rate always adjusts itself in such a way that the total (self-scattering plus
real scattering) rate is a constant in time
t t
self
k k
,
where self
[ k ( t ´)] is the self-scattering rate. The self-scattering mechanism itself is
defined such that the final state before and after scattering is identical. Hence, it has no
effect on the free flight trajectory of a particle when selected as the terminating scattering
mechanism, yet results in the simplification of such that the free flight is given by
t r
r
ln
1
.
The constant total rate (including self-scattering) is chosen a priori so that it is larger
than the maximum scattering encountered during the simulation interval. In the simplest
case, a single value is chosen at the beginning of the entire simulation (constant gamma
method), checking to ensure that the real rate never exceeds this value during the
simulation. Other schemes may be chosen that are more computationally efficient, and
which modify the choice of at fixed time increments [
vii
].
6.1.3 Ensemble Monte Carlo Simulation
The algorithm above may be used to track a single particle over many scattering
events in order to simulate the steady-state behavior of a system. Transient simulation
requires the use of a synchronous ensemble of particles in which the algorithm above is
repeated for each particle in the ensemble representing the system of interest until the
simulation is completed. Figure 7 illustrates an ensemble Monte Carlo simulation in
which a fixed time step, t , is introduced to which the motion of all the carriers in the
system is synchronized. The crosses ( ) illustrate random, instantaneous, scattering
events, which may or may not occur during one time-step. Basically, each carrier is
simulated only up to the end of the time-step, and then the next particle in the ensemble is
treated. Over each time step, the motion of each particle in the ensemble is simulated
independent of the other particles. Nonlinear effects such as carrier-carrier interactions
or the Pauli exclusion principle are then updated at each times step, as discussed in more
detail below.
The non-stationary one-particle distribution function and related quantities such
as drift velocity, valley or subband population, etc., are then taken as averages over the
ensemble at fixed time steps throughout the simulation. For example, the drift velocity in
the presence of the field is given by the ensemble average of the component of the
velocity at the n th time step as
N
j
j
z z
v n t
N
v n t
1
where N is the number of simulated particles and j labels the particles in the ensemble.
This equation represents an estimator of the true velocity, which has a standard error
given by
N
s
,
where
2
is the variance which may be estimated from [Error: Reference source not
found]
N
j
z
j
z
v v
N N
N
1
2
2
2
Similarly, the distribution functions for electrons and holes may be tabulated by counting
the number of electrons in cells of k-space. From Eq. , we see that the error in estimated
average quantities decreases as the square root of the number of particles in the ensemble,
which necessitates the simulation of many particles. Typical ensemble sizes for good
statistics are in the range of 10
4
5
particles. Variance reduction techniques to
decrease the standard error given by Eq. may be applied to enhance statistically rare
events such as impact ionization or electron-hole recombination [Error: Reference source
not found].
6.1.4 Scattering Processes
Free carriers (electrons and holes) interact with the crystal and with each other
through a variety of scattering processes which relax the energy and momentum of the
particle. Based on first order, time-dependent perturbation theory, the transition rate from
an initial state k in band n to a final state k ’ in band m for the j th scattering mechanism is
given by Fermi’s Golden rule [
viii
]
k k
n k m k m k V r n k E E
j j
2
where V j
( r ) is the scattering potential of this process, E k
and E k ’
are the initial and final
state energies of the particle. The delta function results in conservation of energy for
long times after the collision is over, with the energy absorbed (upper sign) or emitted
(lower sign) during the process. Scattering rates calculated by Fermi’s Golden rule above
are typically used in Monte Carlo device simulation as well as simulation of ultrafast
processes. The total rate used to generate the free flight in Eq. , discussed in the previous
section, is then given by
k
k k
k k r k
,
2
, ,
2
,
m
j j
n m V n E E
Scattering Mechanisms
Defect Scattering Carrier-Carrier Scattering Lattice Scattering
Crystal
Defects
Impurity Alloy
Neutral Ionized
Intravalley Intervalley
Acoustic Optical Acoustic Optical
Deformation Nonpolar Polar
potential
Piezo-
electric
Scattering Mechanisms
Defect Scattering Carrier-Carrier Scattering Lattice Scattering
Crystal
Defects
Impurity Alloy
Neutral Ionized
Intravalley Intervalley
Acoustic Optical AcousticAcoustic OpticalOptical
Deformation NonpolarNonpolar PolarPolar
potential
Piezo-
electric
Figure 8. Scattering mechanisms in a typical semiconductor.
Figure 8 lists the scattering mechanisms one should in principle consider in a
typical Monte Carlo simulation. They are roughly divided into scattering due to crystal
defects, which is primarily elastic in nature, lattice scattering between electrons (holes)
and lattice vibrations or phonons, which is inelastic, and finally scattering between the
particles themselves, including both single particle and collective type excitations.
Phonon scattering involves different modes of vibration, either acoustic or optical, as well
as both transverse and longitudinal modes. Carriers may either emit or absorb quanta of
energy from the lattice, in the form of phonons, in individual scattering events. The
designation of inter- versus intra-valley scattering comes from the multi-valley band-
structure model mentioned in Section ??, and refers to whether the initial and final states
are in the same valley or in different valleys. The scattering rates
k k
n , ; m ,
j and
n , k
j
are calculated using time dependent perturbation theory using Fermi’s rule, Eqs.
and , and the calculated rates are then tabulated in a scattering table in order to select the
type of scattering and final state after scattering as discussed earlier.
6.1.5 Listing of the Bulk Monte Carlo Code for GaAs
The Monte Carlo program for the electrons for a three-valley semiconductor (such as GaAs in
this case) is initiated through the execution of the main() routine that first initializes various
parameters used in the code by calling the readin() subroutine and by the construction of the
scattering table by calling the sc_table() subroutine. The next step in the initialization procedure
is to initialize the carrier energies and wavevectors by calling the init() routine. Having initialized
the carriers, a check is performed by calling the histograms() subroutine that calculates the
energy and wavevector histograms of the carriers where one can see that carriers have been
initialized according to Maxwell-Boltzmann distribution. After the initialization procedure is
finished, a free-flight-scatter loop is performed until the end of the desired time window that is
specified in the subroutine readin(). Within each time step, the free_flight_scatter() subroutine
is called first, then the histograms() subroutine is called to check the change in the distribution
function under the application of uniform electric field and finally the write() subroutine is called
to record the values of the carrier wavevectors, energy, valley index and other relevant parameters
that can give us idea whether the time evolution of the variables is correct or not. A flow chart of
the main() program is summarized in Figure ??.
parameter(n_lev=1000, nele=1000)
C
C n_lev => # of energy levels in the scattering table
C nele => # total number of electrons that are simulated
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
common
&/time_1/dt,dtau,tot_time
character*40 file_name
integer nsim
C Read parameters form the input file
call readin(n_lev)
C Start the Monte Carlo Simulation
C Calculate the scattering table
call sc_table(n_lev)
C Initialize carriers
nsim = nele
call init(nsim)
file_name = 'initial_distribution'
call histograms(nsim,file_name)
time = 0.
j = 0
flag_write = 0.
do while(time.le.tot_time)
j = j + 1
time=dt*float(j)
call free_flight_scatter(n_lev,nsim)
file_name = 'current_distribution'
call histograms(nsim,file_name)
call write(nsim,j,time,flag_write)
flag_write = 1.
enddo! End of the time loop
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C WRITE AVERAGES IN FILES
C (these averages correspond to one time step)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine write(nsim,iter,time,flag_write)
parameter(n_val = 3)
common
&/ran_var/iso
&/pi/pi,two_pi
&/fund_const/q,h,kb,am0,eps_
&/dri/qh
&/ek/am(3),smh(3),hhm(3)
&/nonp/af(3),af2(3),af4(3)
&/variables/p(20000,7),ip(20000),energy(20000)
integer nsim,iter
real kb,time
character*40 file_1,file_
character*40 file_3,file_4,file_
real nvaly(n_val)
real velx_sum(n_val)
real vely_sum(n_val)
real velz_sum(n_val)
real velocity_x(n_val)
real velocity_y(n_val)
real velocity_z(n_val)
real sume(n_val)
C Open output files
file_1 = 'vx_time_averages'
file_2 = 'vy_time_averages'
file_3 = 'vz_time_averages'
file_4 = 'energy_time_averages'
file_5 = 'valley_occupation'
open(unit=1,file=file_1,status='unknown')
open(unit=2,file=file_2,status='unknown')
open(unit=3,file=file_3,status='unknown')
open(unit=4,file=file_4,status='unknown')
open(unit=5,file=file_5,status='unknown')
do i = 1,n_val
velx_sum(i) = 0.
vely_sum(i) = 0.
velz_sum(i) = 0.
sume(i) = 0.
nvaly(i) = 0
enddo
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C SAVE HISTOGRAMS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine histograms(nsim,file_name)
common
&/variables/p(20000,7),ip(20000),energy(20000)
integer nsim
real kx,ky,kz,e
character*30 file_name
open(unit=6,file=file_name,status='unknown')
write(6,*)'kx ky kz energy'
do i = 1, nsim
kx = p(i,1)
ky = p(i,2)
kz = p(i,3)
e = energy(i)
write(6,*)kx,ky,kz,e
enddo
close(6)
return
end
The readin() subroutine, which contains all the relevant parameters for modeling
electronic transport in GaAs bulk material system, is listed next. In our model, we have
assumed three conduction band valleys: , L and X as depicted in Figure ??. All
variables defined in the readin() subroutine have descriptive names and can be easily
identified.
k-vector
-valley
X-valley [100]
L-valley [111]
Conduction bands
Valence bands
Figure ??. Energy band model for GaAs. The L-valley is at the (111) point and there are 8 equivalent [111]
directions. Since these valleys are shared between Brillouin zones there are a total of four equivalent L
valleys. The X-valleys are at the [100] direction and since there are 6 equivalent [100] directions and the
valleys are shared between Brillouin zones, there are 3 equivalent X valleys.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C READ INPUT PARAMETERS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine readin(n_lev)
common
&/ran_var/iso
&/pi/pi,two_pi
&/fund_const/q,h,kb,am0,eps_
&/dri/qh
&/temp/tem,Vt
&/ek/am(3),smh(3),hhm(3)
&/valley_splitting/split_L_gamma,split_X_gamma
&/equiv_valleys/eq_valleys_gamma,eq_valleys_L,eq_valleys_X
&/nonp/af(3),af2(3),af4(3)
&/scatt_par/emax,de,w(10,3),tau_max(3),max_scatt_mech(3)
&/dielec_func/eps_high,eps_low
&/density/density,sound_velocity
&/time_1/dt,dtau,tot_time
&/force/fx,fy,fz
&/select_acouctic/acoustic_gamma,acoustic_L,acoustic_X
&/select_Coulomb/Coulomb_scattering
&/select_polar/polar_gamma,polar_L,polar_X
&/select_intervalley_1/intervalley_gamma_L
&/select_intervalley_2/intervalley_gamma_X
fz = 0.
C Define Masses for different valleys
rel_mass_gamma = 0.
rel_mass_L = 0.
rel_mass_X = 0.
C Define non-parabolicity factors
nonparabolicity_gamma = 0.
nonparabolicity_L = 0.
nonparabolicity_X = 0.
c nonparabolicity_gamma = 0
c nonparabolicity_L = 0
c nonparabolicity_X = 0
C Define valley splitting and equivalent valleys
split_L_gamma = 0.
split_X_gamma = 0.
eq_valleys_gamma = 1.
eq_valleys_L = 4.
eq_valleys_X = 3.
C Define low-ferquency and high-frequency dielectric constants
eps_high = 10.
eps_low = 12.
eps_high = eps_high*eps_
eps_low = eps_low*eps_
C Define crystal density and sound velocity
density = 5370.
sound_velocity = 5.22E
C Define parameters for the scattering table
emax=1.
de=emax/float(n_lev)
C Select scattering mechanisms
c Coulomb_scattering = 1
acoustic_gamma = 1
acoustic_L = 1
acoustic_X = 1
polar_gamma = 1
polar_L = 1
polar_X = 1
intervalley_gamma_L = 1
intervalley_gamma_X = 1
intervalley_L_gamma = 1
intervalley_L_L = 1
intervalley_L_X = 1
intervalley_X_gamma = 1
intervalley_X_L = 1
intervalley_X_X = 1
C Define coupling constants
sigma_gamma = 7.01! [eV]
sigma_L = 9.2! [eV]
sigma_X = 9.0! [eV]
polar_en_gamma = 0.03536! [eV]
polar_en_L = 0.03536! [eV]
polar_en_X = 0.03536! [eV]
DefPot_gamma_L = 1.8E10! [eV/m]
DefPot_gamma_X = 10.E10! [eV/m]
DefPot_L_gamma = 1.8E10! [eV/m]
DefPot_L_L = 5.E10! [eV/m]
DefPot_L_X = 1.E10! [eV/m]
DefPot_X_gamma = 10.E10! [eV/m]
DefPot_X_L = 1.E10! [eV/m]
DefPot_X_X = 10.E10! [eV/m]
phonon_gamma_L = 0.0278! [eV]
phonon_gamma_X = 0.0299! [eV]
phonon_L_gamma = 0.0278! [eV]
phonon_L_L = 0.029! [eV]
phonon_L_X = 0.0293! [eV]
phonon_X_gamma = 0.0299! [eV]
phonon_X_L = 0.0293! [eV]
phonon_X_X = 0.0299! [eV]
C Map parameters into internal variables used in the code
am(1) = rel_mass_gamma*am
am(2) = rel_mass_L*am
am(3) = rel_mass_X*am
af(1) = nonparabolicity_gamma
af(2) = nonparabolicity_L
af(3) = nonparabolicity_X
