InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Computer and Information Science » Numerical Analysis and Scientific Computing » "Theory and Applications of Monte Carlo Simulations", book edited by Victor (Wai Kin) Chan, ISBN 978-953-51-1012-5, Published: March 6, 2013 under CC BY 3.0 license. © The Author(s).

Chapter 5

Comparative Study of Various Self-Consistent Event Biasing Schemes for Monte Carlo Simulations of Nanoscale MOSFETs

By Shaikh Ahmed, Mihail Nedjalkov and Dragica Vasileska
DOI: 10.5772/53113

Article top


Three-dimensional Monte Carlo device simulator (MCDS 3-D) and the associated models/kernels used in this work.
Figure 1. Three-dimensional Monte Carlo device simulator (MCDS 3-D) and the associated models/kernels used in this work.
Rates for various scattering mechanisms as a function of electron energy used in the MCDS 3-D simulator.
Figure 2. Rates for various scattering mechanisms as a function of electron energy used in the MCDS 3-D simulator.
Variation of the quantum barrier field (QBF) as a function of distance from the Si/SiO2 interface and wavevector ky along the depth (left panel: low energy electrons, right panel: high energy electrons).
Figure 3. Variation of the quantum barrier field (QBF) as a function of distance from the Si/SiO2 interface and wavevector ky along the depth (left panel: low energy electrons, right panel: high energy electrons).
Typical plot of cumulative net charge at the source and the drain contacts.
Figure 4. Typical plot of cumulative net charge at the source and the drain contacts.
Weighting scheme used in the temperature-biasing method.
Figure 5. Weighting scheme used in the temperature-biasing method.
a) Weight distribution for T = 450K, and (b) Biasing decreases numerical particle number in the source/drain regions while increases the number in the channel region. Here 3n stands for numerical particles under bias = 3.
Figure 6. a) Weight distribution for T = 450K, and (b) Biasing decreases numerical particle number in the source/drain regions while increases the number in the channel region. Here 3n stands for numerical particles under bias = 3.
Particle split method and the distribution of particles in the device active region.
Figure 7. Particle split method and the distribution of particles in the device active region.
Enhancement of channel statistics: reduction of standard deviation in number of electrons in a particular cell in the channel region.
Figure 8. Enhancement of channel statistics: reduction of standard deviation in number of electrons in a particular cell in the channel region.
Biasing recovers precisely the self-consistent (a) average sheet density and (b) average kinetic energy of the electrons.
Figure 9. Biasing recovers precisely the self-consistent (a) average sheet density and (b) average kinetic energy of the electrons.
Comparison of channel currents obtained from (1) biased e-a rates (top panel), and (2) biased boundary distribution (bottom panel) methods from velocity consideration.
Figure 10. Comparison of channel currents obtained from (1) biased e-a rates (top panel), and (2) biased boundary distribution (bottom panel) methods from velocity consideration.
Terminal currents obtained from various methods by particle counting.
Figure 11. Terminal currents obtained from various methods by particle counting.
Numerical particle weight distribution in (a) e-a biasing, and (b) e-a/split biasing.
Figure 12. Numerical particle weight distribution in (a) e-a biasing, and (b) e-a/split biasing.

Comparative Study of Various Self-Consistent Event Biasing Schemes for Monte Carlo Simulations of Nanoscale MOSFETs

Shaikh Ahmed1, Mihail Nedjalkov2 and Dragica Vasileska1

1. Introduction

1.1. Semiclassical electron transport

Semiclassical Boltzmann transport has been the principal theory in the field of modeling and simulation of semiconductor technology since its early development. To date, most commercial device simulations including the full-band Monte Carlo (FBMC) method are based on the solution of the Boltzmann transport equation (BTE) and its simplified derivatives such as the hydrodynamic (HD) equations and the drift-diffusion (DD) model. The Boltzmann transport equation expresses the global variation of the non-equilibrium distribution function, f(r,k,t) , under the influence of various applied and built-in forces. In its most general form and for non-parabolic bands the BTE reads:


where v is the carrier group velocity. The terms on the left-hand side describe the change in the distribution function with respect to time, spatial gradients, and applied fields. The right-hand side represents the dissipation terms in the system, which account for the change of the distribution function due to various scattering mechanisms that balance the driving terms on the left. The Boltzmann equation is valid under the assumptions of semiclassical transport as follows [1]: (1) effective mass approximation; (2) First Born approximation for the collisions in the limit of small perturbation for the electron-phonon interaction and instantaneous collisions; (3) Scattering probability is independent of external forces; (4) Particle interactions are uncorrelated and forces are constant over distances comparable to the electron wave function; and (4) The phonons are usually treated as being in equilibrium, although the condition of non-equilibrium phonons may be included through an additional phonon transport equation.

Full analytical solutions of the BTE are possible only under very restrictive assumptions. The two classes of computational/numerical techniques that are used to solve the BTE are as follow [2][3][4][5]: (1) First, in deterministic methods, the BTE is discretized using a variety of methods (such as discrete ordinates, spherical harmonics, collision probabilities, nodal methods) and then solved numerically. However, direct numerical methods are limited by the complexity of the equation, which in the complete three-dimensional (3D) time-dependent form requires seven independent variables for time, space and momentum; (2) The second class of techniques, named Monte Carlo methods, construct a stochastic model in which the expected value of a certain random variable is equivalent to the value of a physical quantity to be determined [6][7][8][9]. The expected value is estimated by the average of many independent samples representing the random variable. Random numbers, following the distributions of the variable to be estimated, are used to construct these independent samples. There are two different ways to construct a stochastic model for Monte Carlo calculations. In the first case, the physical process is stochastic and the Monte Carlo calculation involves a computational simulation of the real physical process. This is achieved by tracing the trajectories of individual carriers as they are accelerated by the electric field and experience random scattering events. The particle movements between scattering events are usually described by the laws of classical mechanics, while the probabilities of the various scattering processes and the associated transition rates are derived from quantum mechanical calculations. The randomness of the events is treated in terms of computer generated random numbers, distributed in such a way as to reflect these probabilities. In the other case, a stochastic model is constructed artificially, such as the solution of deterministic equations by Monte Carlo [10]. Both the deterministic and the Monte Carlo stochastic methods have computational errors. Deterministic methods are computationally fast but less accurate; whereas Monte Carlo methods are computationally slow yet arbitrarily accurate. A great advantage of particle-based Monte Carlo methods is that it provides a unique and thorough insight into the underlying device physics and carrier transport phenomena [11]. This insight originates from the very nature of the method employed that relies on a detailed simulation of a large ensemble of discrete charge carriers in the device active region. Therefore, the full time-dependent distribution function and related macroscopic quantities of interest can be extracted directly from the simulation.

In the last decade, with the continued downsizing of device dimensions into the nanoscale regime, many new and challenging questions have emerged concerning the physics and characterization of these small devices. However, contrary to the technological advances, present state-of-the-art in device simulation is lacking in the ability to treat the new challenges pertaining to device scaling [12]. First, for devices where gradual channel approximation (analytical modeling) cannot be used due to the two-dimensional nature of the electrostatic potential and the electric fields driving the carriers from source to drain, drift-diffusion models have been exploited. These models are valid, in general, for large devices in which the fields are not that high so that there is no degradation of the mobility due to the electric field. The validity of the drift-diffusion models can be extended to take into account the velocity saturation effect with the introduction of field-dependent mobility and diffusion coefficients. When velocity overshoot becomes important, drift diffusion model is no longer valid and hydrodynamic model may be considered. The hydrodynamic model has been the workhorse for technology development and several high-end commercial device simulators have appeared including Silvaco, Synopsys, Crosslight, etc. The advantage of the hydrodynamic model is that it allows quick simulation runs. However, standard hydrodynamic models do not provide a sufficiently accurate description since they neglect significant contributions from the tail of the phase space distribution function in the active region of the device. Also, the velocity overshoot in hydrodynamic models depends upon the choice of the energy relaxation time. The smaller is the device, the larger is the deviation when using the same set of energy relaxation times. In addition, the energy relaxation times are material, device geometry and doping dependent parameters, so their determination ahead of time is not possible. To avoid the problem of the proper choice of the energy relaxation times, a direct solution of the Boltzmann Transport Equation using the Monte Carlo method is the best method of choice. To-date, most semiclassical device simulations have been based on stochastic Monte Carlo solution methods, which involve the simulation of particle trajectories rather than the direct solution of partial differential equations. An additional problem arises in small structures, where one must begin to worry about the effective size of the carriers themselves. In the nanoscale regime, the charge transport is expected to be dominated by quantum effects throughout the active region. Quantum effects in the surface potential will have a profound impact on both the amount of charge, which can be induced by the gate electrode through the gate oxide, and the profile of the channel charge in the direction perpendicular to the surface (the transverse direction). Also, because of the two-dimensional confinement of carriers in the channel region, the mobility (or microscopically speaking, the carrier scattering) will be different from the three-dimensional (bulk) case. A well-known approach to study the impact of spatial confinement on mobility behavior is based on the self-consistent solution of the 3D Poisson–1D Schrödinger–3D Monte Carlo, which demands enormous computational resources as it requires storage of position dependent scattering tables that describe carrier transition between various subbands. In the smallest size devices, for a full quantum mechanical description, various quantum formalisms based on density matrices [13], Wigner functions [14], Feynman path integrals [15], and non-equilibrium Green’s functions (NEGF) [16] have been developed and proposed with varying success to address these issues. The Green’s functions approach is the most exact, but at the same time appears, from the historical literature perspective, as the most difficult of all. In contrast to, for example, the Wigner function approach (which is Markovian in time), the Green's functions method allows one to consider simultaneously correlations in space and time or space and energy, both of which are expected to be important in nanoscale devices. However, today, although the NEGF transport formalism has been well-established, accurate and full three-dimensional modeling of scattering-dominated transport in realistically-sized semiconductor devices using the NEGF approach is prohibitively expensive and the computational burden needed for its actual implementation are perceived as a great challenge. A successful utilization of the Green's function approach available commercially is the NEMO (Nano-Electronics Modeling) simulator [17], which is effectively 1-D and is primarily applicable to resonant tunneling diodes. On the other hand, within the semiclassical Boltzmann transport formalism, for modeling quantum effects in nanostructures, recently, different forms of quantum effective potentials have been proposed [18][19][20][21][22], and when used in conjunction with the BTE, provide satisfactory level of accuracy. The idea of quantum potentials is quite old and originates from the hydrodynamic formulation of quantum mechanics, first introduced by de Broglie and Madelung [23][24], and later developed by Bohm [25]. The work presented here focuses mainly on the modeling and simulations of nanoscale MOSFET structures within a quantum-corrected Monte Carlo transport framework.

2. The Monte Carlo method

Monte Carlo (MC) method was originally used and devised by Fermi, Von Neumann and Stanislaw Ulam to solve the BTE for transport of neutrons in the fissile material of the atomic bomb during the Manhattan Project of World War II [26]. Since these pioneering times in the mid 1940's, the popularity and use of the MC method has grown with the increasing availability of faster and cheaper digital computers. Its application to the specific problems of high-field electron transport in semiconductors is first due to Kurosawa [27] in 1966. Shortly afterwards the Malvern, UK, group [28] provided the first wide application of the method to the problem of the Gunn effect in GaAs. Applications to Si and Ge flourished in the 1970s, with an extensive work performed at the University of Modena, Italy. In the mid-1970s, a physical model of silicon was developed, capable of explaining major macroscopic transport characteristics. The band structure models were represented by simple analytic expressions accounting for nonparabolicity and anisotropy. The review articles by C. Jacoboni and L. Reggiani [29] and by Peter J. Price [30] provide a comprehensive and deep historical and technical perspective.

The Monte Carlo method is the most popular method used to solve the Boltzmann transport equation without any approximation made to the distribution function. In the Monte Carlo method, particles are used to represent electrons (holes) within the device. For bulk simulations, the momentum and energy of the particles (electrons and/or holes) are continuously updated. For actual device simulations, the real space position of the particle is updated as well. As time evolves, the updated momentum (and corresponding energy) is calculated from the various forces applied on the particle for that time step. The general concept of a Monte Carlo simulation is that the electrons (holes) are accelerated by an electric field until they cover a predetermined free-flight time. At the end of this, a scattering mechanism (due to impurities, lattice vibrations, etc.) is randomly chosen based on its relative frequency of occurrence, which randomize the momentum and energy of charge particles in time. The basic steps in a typical Monte-Carlo device simulation scheme include (see Refs. [31][32][33]):

  1. Initialization.

  2. Solution of the field (Poisson) equation: determine forces on electrons.

  3. Simulate the electron dynamics (free-flight and scattering):

    1. Accelerate the electrons.

    2. Determine whether the electron experiences a scattering/collision.

    3. If electron scatters, select the scattering mechanism.

    4. Update the electron position, energy and wavevector.

  1. Charge assignment to the numerical nodes/grids.

  2. Compute measurable quantities such as average energy, average velocity, and mobility of the particle ensemble.

  3. Repeat steps 2-5 for each iteration.

The initialization includes the calculation/setting of various material parameters such as the electronic band structure, scattering rates and (particularly for device simulation) the definition of the computational domain and initial electron distributions. In order to determine device characteristics, such as drain current in a MOSFET, it is necessary to solve the carrier transport equations using the device parameters (including doping concentration and type, channel length, channel width, and oxide thickness) and boundary conditions (such as terminal voltages) for the particular device being simulated. The device is divided into grids/cells by a numerical discretization technique (usually finite difference method). In the actual implementation, for each dimension, maximum number of node points is set as a program parameter and the relevant array sizes are optimally (often dynamically) allocated. Provision may be kept for using nonuniform mesh spacing along different axial dimensions of the device. The source and drain contact charges are then calculated by multiplying the doping density with the volume of a single cell. Then, the Poisson equation is solved for the applied gate bias only (keeping source/drain/substrate bias equal to zero) and the resulting potential distribution is used to populate/generate carriers in each of the cells in the active region through the use of appropriate equilibrium distribution function (usually Fermi-Dirac). It is always convenient to begin the simulation in thermodynamic equilibrium where the solution is known. When an electron is introduced in the device active region, the real space coordinates of the electron are determined based upon the position of the node where it is being added. Under the assumption of charge-neutrality, the total number of free carriers within the device must equal the total number of ionized dopants. In the non-equilibrium simulation, the remaining boundary biases are applied at different contacts/terminals (source/drain/substrate). The program control is then transferred to the Monte Carlo iterative transport/dynamics kernel. Here, for each pre-selected time step (typically fractions of femtosecond) the carriers undergo the free-flight-scatter sequence. During a time interval, Δt, each electron is accelerated according to Newton’s second law of motion. For a semiconductor with ellipsoidal constant energy surfaces, such as silicon, the effective mass is a tensor quantity. Using Hering-Vogt transformation [31]

ki=mmiki  i=x,y,z

where mi is the effective mass of the particle in the ith (x, y, or z) direction, one can make a transformation from k-space to k*-space in which the constant energy surfaces are spherical. In the original k-space, due to the mass anisotropy, the constant energy surfaces are ellipsoidal. To determine the time between scattering events (free-flight time), the electron energy at the time of scattering must be known, since the scattering rate is a function of the electron energy. The rates for various scattering mechanisms included into the Monte Carlo model are pre-tabulated (in number of collisions per second) as a function of energy, which are usually calculated quantum mechanically using perturbation theory. Once the type of scattering terminating the free flight is selected, the final energy and momentum (as well as band or subband) of the particle due to this type of scattering must be selected. For elastic scattering processes such as ionized impurity scattering, the energy before and after scattering is the same. In contrast, for the interaction between electrons and the vibrational modes of the lattice described as quasi-particles known as phonons, electrons exchange finite amounts of energy with the lattice in terms of emission and absorption of phonons. The final k-vector is then chosen randomly. After a new k-vector is chosen, a new free-flight time is determined. The electron is then accelerated for the remainder of the time interval or until it encounters another scattering event.

To simulate an actual device such as a MOSFET, the device boundaries must be treated properly. The Ohmic contacts are often assumed to be perfect absorbers, so carriers that reach them simply exit the device. At the end of each time step, thermal electrons are injected from the contacts to maintain space-charge neutrality therein. The noncontacted free surfaces are treated as reflecting boundaries. For field-effect transistors, roughness at the surface of the channel can cause scattering. A simple approach is to treat some fraction of the encounters with the surface as specular scattering events and the remainder as diffusive scattering events. The specific fraction is usually selected to match transport measurements, such as the low-field mobility as a function of the transverse electric field. A carrier deletion scheme is also implemented at this stage. When completed, randomly distributed charges from the Monte Carlo simulations are assigned to the nearest node points using an appropriate charge-assignment method. A proper charge assignment scheme must ensure proper coupling (match) between the charged particles and the actual Coulomb forces acting on the particles, and thereby maintain zero self-force and a good spatial accuracy of the forces. Once the updated grid charge distribution is known, Poisson equation is then solved to determine the resulting potential distribution within the device. It is important to note that the electron concentration (in an n-MOSFET simulation) is not updated based upon the potential at the node point. In contrary, the hole concentration at each node point is calculated using the updated node potential with the assumption that the hole quasi-Fermi level is equal to that of the electrons (quasi-static approximation). The error introduced with the assumption that the quasi-Fermi levels for electrons and holes are equal is small since the hole concentration is negligible in the active region of the (n-MOSFET) device. The force on each electron is then interpolated from the nearest node points. The resultant electric field is then used to drive the carriers during the free-flight in the next time step. The whole process is repeated for several thousand iterations until the steady state is achieved. In summary, the Monte Carlo algorithm consists of generating random free-flight times for each particle, choosing the type of scattering occurring at the end of the free-flight, changing the final energy and momentum of the particle after scattering, and then repeating the procedure for the next free flight. At any time, sampling the particle motion throughout the simulation allows for the statistical estimation (by averaging over the particles within each slab of the active region) of physically interesting and key device quantities such as the single particle distribution function, average carrier density, velocity, mobility, and energy versus position. Lundstrom [32], Tomizawa [33], and Kunikiyo et al. [34] discuss the application of this approach to the simulation of two-dimensional transistors. A detailed description of a 3D Monte Carlo device simulator can be found in Refs. [31] and [35].

3. Description of the Monte Carlo device simulator

The major computational components of the 3-D Monte Carlo device simulator (MCDS 3-D) used in this work are shown in Figure 1. Regarding the Monte Carlo transport kernel, intravalley scattering is limited to acoustic phonons. For the intervalley scattering, both g- and f-phonon processes have been included. It is important to note that, by group symmetry considerations, the zeroth-order low-energy f- and g-phonon processes are forbidden. Nevertheless, three zeroth-order f-phonons and three zeroth-order g-phonons with various energies are usually assumed [28]. This selection rule has been taken into account and two high-energy f- and g-phonons and two low-energy f- and g-phonons have been considered. The high-energy phonon scattering processes are included via the usual zeroth-order interaction term, and the two low-energy phonons are treated via a first-order process [36]. The first-order process is not really important for low-energy electrons but has a significant contribution for high-energy electrons. The low-energy phonons are important in achieving a smooth velocity saturation curve, especially at low temperatures. The phonon energies and coupling constants in this model are determined so that the experimental temperature-dependent mobility and velocity-field characteristics are consistently recovered [37].

Figure 2 shows the rates for various scattering mechanisms as a function of electron energy used in the MCDS 3-D simulator. At present, impact ionization and surface-roughness scattering are not included in the model. Impact ionization is neglected, as, for the drain biases used in the simulation, electron energy is typically insufficient to create electron-hole pairs. Also, since it is not quite clear how surface roughness scattering can be modeled when carriers are displaced from the interface due to the quantum confinement effects, it is believed that its inclusion is most likely to obscure the quantum confinement effects. Also, band-to-band tunneling and generation and recombination mechanisms have not been included in the simulations.


Figure 1.

Three-dimensional Monte Carlo device simulator (MCDS 3-D) and the associated models/kernels used in this work.


Figure 2.

Rates for various scattering mechanisms as a function of electron energy used in the MCDS 3-D simulator.

The Incomplete Lower-Upper (ILU) decomposition method has been employed for the solution of the 3D Poisson equation. To treat full Coulomb (electron-ion and electron-electron) interactions properly, the simulator implements two real-space molecular dynamics (MD) schemes: the particle-particle-particle-mesh (P3M) method and the corrected Coulomb approach. The effective force on an electron is computed as a combination of the short-range molecular dynamics force and the long-range Poisson force. The implementation details of these models and methodologies have been discussed in Ref. [35].


Figure 3.

Variation of the quantum barrier field (QBF) as a function of distance from the Si/SiO2 interface and wavevector ky along the depth (left panel: low energy electrons, right panel: high energy electrons).

Quantum mechanical size-quantization effects have been accounted for via a parameter-free effective potential scheme [38]. The approach is based on a perturbation theory around thermodynamic equilibrium and derived from the idea that the semiclassical Boltzmann equation with the quantum corrected potential and the Wigner equation should possess the same steady state. It leads to an effective potential/field, which takes into account the discontinuity at the Si/SiO2 barrier interface due to the difference in the semiconductor and the oxide affinities. It possesses no fitting parameters, as the size of the electron (wavepacket) is determined from its energy. The resultant quantum potential is, in general, two-degrees smoother than the original Coulomb and barrier potentials of the device, i.e. possesses two more classical derivatives, which essentially eliminates the problem of the statistical noise. The calculated quantum barrier field (QBF) for low-energy (left column) and high-energy (right column) electrons are shown in Figure 3 having the following salient features: (1) QBF decays almost exponentially with distance from the Si/SiO2 interface proper; (2) QBF increases with increasing the wavevector of the carriers along the normal direction; (3) The contour plots clearly reveal the fact that the electrons with lower momentum feel the quantum field far from the interface proper, and (4) A similar trend is also observed with the variation in electron energy—electrons with higher energy can reach the vicinity on the interface, thus, behaving as classical point-like particles.

The device output current is determined using two different yet consistent methods [31]: (1) In the first approach, the charges (electrons) entering and exiting each terminal/contact are tracked with time. The net number of charges over a period of the simulation experiment is then used to calculate the terminal current. The net charge crossing a terminal boundary is determined by


where q is the electronic charge, nabsorbed is the number of particles that are absorbed by the contact, ninjected is the number of particles that have been injected at the contact, Ey is the vertical field at the contact. The second term in (3) on the right-hand-side accounts for the displacement current due to the changing field at the contact. Eq. (3) assumes that the contact is at the top of the device and that the fields in the x and z direction are negligible. The net charge calculations are made for both the source and the drain contacts. Figure 4 shows a typical plot of net terminal charges for a 2 ps simulation run. The terminal current is measured from the slope of Q(t)versus time plot. In steady state, the current can be calculated by


where nnet is the net number of particles (electrons) exiting the contact over a fixed period of time, Δt. It is important to note that in steady state the source current must be equal to the drain current. Therefore, for proper current measurements only those portions of the curves where the two slopes are equal, is considered. Note that the method is quite noisy, due to the discrete nature of the electrons; (2) In the second method, the channel current is calculated from the sum of the electron velocities in a portion of the channel region of the device. The electron (channel) current density through a cross-section of the device is given by

where vd is the average electron drift velocity and n is the carrier concentration. If there are a total of N particles in a differential volume, dV=dLdA,the current found by integrating (5) over the cross-sectional area, dA, is



where vx(n) is the velocity along the channel of the nth electron. The device is divided into several sections along the x-axis, and the number of electrons and their corresponding velocity is added for each section after each free-flight. The total x-velocity in each section is then averaged over several time-steps to determine the current for that section. Total device current can be determined from the average of several sections, which gives a much smoother and noise-free result compared to counting terminal charges. By breaking the device into sections, individual section currents can be compared to verify that there is conservation of particles (constant current) throughout the device. Note that the sections in the vicinity of the source and drain regions may have a large y-component in their velocity and should be excluded from the current calculations. Finally, by using several sections in the channel, the average energy and velocity of electrons along the channel can be observed to ensure proper physical characteristics.


Figure 4.

Typical plot of cumulative net charge at the source and the drain contacts.

4. Statistical enhancement: The self-consistent event biasing scheme

Statistical enhancement in Monte Carlo simulations aims at reduction of the time necessary for computation of the desired device characteristics. Enhancement algorithms are especially useful when the device behavior is governed by rare events in the transport process. Such events are inherent for sub-threshold regime of device operation, simulations of effects due to discrete dopant distribution as well as tunneling phenomena. Virtually all Monte Carlo device simulators with statistical enhancement use population control techniques [39]. They are based on the heuristic idea for splitting of the particles entering a given phase space region Ω of interest. The alternative idea―to enrich the statistics in Ω by biasing the probabilities associated with the transport of classical carriers―gives rise to the event-biasing approach. The approach, first proposed for the Ensemble Monte Carlo technique (time-dependent problem) [40], has been recently extended for the Single Particle Monte Carlo technique (stationary problem) [41]. In the next section, the basic steps of derivation of the approach in presence of both initial and boundary conditions has been discussed. Utilized is the linearity of the transport problem, where Coulomb forces between the carriers are initially neglected. The generalization of the approach for Hartree carriers has been established in the iterative procedure of coupling with the Poisson equation. Self-consistent simulation results are presented and discussed in the last section.

4.1. Theory of Event Biasing

The Ensemble Monte Carlo (EMC) technique is designed to evaluate averaged values A of generic physical quantities a such as carrier density and velocity given by


Where Q=(k,r,t)and Equation (8) denotes the integration over the phase space and time t(0,), and A=aθΩδ(tτ) introduces the indicator θΩ of the phase space domain, where the mean value is evaluated at time τ. Equation (8) is the usual expression for a statistical mean value, augmented by a time integral with the purpose to conveniently approach the formal theory of integral equations. It has been shown that the Boltzmann equation can be formulated as a Fredholm integral equation of a second kind with a free term f0. The latter is determined by the initial condition in evolution problems [38][42] or, in the case of stationary transport, by the boundary conditions [42]. The second equality in (8) follows from the relationship between an integral equation and its adjoint equation. It shows that the mean value A is determined by f0 and by the solution of the adjoint Boltzmann equation

K=S(k,k,r)ettλ(K(y),R(y))dyθD(r)δ(r r)θ(tt)

where S is the usual scattering rate from lattice imperfections, λ is the total out-scattering rate, θD is the device domain indicator, which is discussed later, θ is the Heaviside function and the trajectories, initialized by (k,r,t), are formulated with the help of the electrical force F and the velocity v as




If both, initial fi and boundary fb conditions are taken into account, it can be shown that f0 becomes


While fi is defined only at the initial time t=0, the function fb is defined only at the device boundary Γ and for values of k such that the corresponding velocity inwards the domain, D. v is the velocity component normal to Γ so that a velocity-weighted distribution drives the particle flux, injected into the device at times tbt. f0 in (13) governs both the transient and the stationary behavior of a device. The latter is established in the long time limit, provided that fb is time independent. Usually fb is assumed to be the equilibrium distribution function.

A recursive replacement of equation (9) into itself gives rise to the von-Neumann expansion, where the solution g is presented as a sum of the consecutive iterations of the kernel on A. If replaced in (8), the expansion gives rise to the following series for A:


Consider the second term in (14) augmented with the help of two probabilities P0 and P to become expectation [41] value of a random variable (r.v.):


It takes values determined by the second half with a probability given by the product in the first half in equation (15). A2 is evaluated according to the numerical Monte Carlo theory as follows. P0 and P are used to construct numerical trajectories: (i) P0(Q) selects the initial point Q of the trajectory. (ii) P(Q,Q) selects the next trajectory point Q provided that Q is given. The fraction W2 in front of A, called weight, is a product of weight factors f0P0, and KP evaluated at the corresponding points Q0Q1Q2, selected by application of P0PP. The sample mean of N realizations of the r.v., calculated over N trajectories (QQ1Q2)n,n=1...N, estimates the mean value A2:


The iterative character of the multiple integral (15) has been used to introduce a consecutive procedure for construction of the trajectories. It can be shown that a single trajectory, obtained by successive applications of P, contributes to the estimators of all terms in (14) simultaneously i.e. the procedure is generalized in (16) for a direct evaluation of A. Next, one establishes the link between (16) and the EMC technique, which is due to particular choice of the initial, P0B, and transition, PB, densities. PB, which can be deduced from (10), is a product of the conditional probabilities for free-flight and scattering, associated with the evolution of the real carriers. The ratio K/PB is then the domain indicator θD which takes values 1 (one) if the trajectory belongs to D and 0 (zero) otherwise. The choice of P0B is complicated by the presence of both initial and boundary terms in (13). They decompose (16) into two terms which are evaluated separately as


The initial probability P0B for each estimator is obtained from fi and vfb respectively, with the help of two normalization factors: the number of initial carriers Ni and the total number NJ of the injected particles into the device. The ratio f0/P0B for each of the estimators becomes Ni and NJ respectively, and can be eliminated by the choice N1=Ni and N2=NJ. The two sums can be merged back to give


Equation (18) accounts that only trajectories which belong to D give contributions. As only the endpoint of such trajectories matters for the estimator, we speak about particles inside the device. Nτ is the number of such particles at time τ, and θΩ(n) is 1 or 0 if the n-th particle is inside or outside Ω. All particles have weight unity and evolve as real Boltzmann carriers and the EMC technique for transport problems posed by initial and boundary conditions is recovered. A choice of alternative probabilities is called event biasing. Biased can be the probabilities for initial and/or boundary distributions, free-flight duration, type of scattering and the selection of the after-scattering state direction. It can be shown that (18) is generalized to A=n=1NτbWnθΩ(n)an where the position of the Nτb biased particles is accounted in θΩ.

The Boltzmann equation for Coulomb carriers becomes nonlinear via the interaction component F(f)(r,t) of the electric force. As the results of the previous section are based on the linearity of the integral equations involved, it is no more possible to apply the steps used to derive the event biasing. The solution is sought in the iterative procedure of coupling of the EMC technique with the Poisson equation. The latter is discretized, as stated earlier, by a decomposition of the device region into mesh cells, ψl. The particle system is evolved in time intervals Δt0.1 fs. At the end of each time step, at say time τ, the charge density eC(rl,τ) is calculated and assigned to the corresponding grid points. One uses the relation between Cl and the distribution function fl,m=f(rl,km,τ), which is estimated with the help of (18) by introducing a mesh ϕm in the wavevector space, (Ωl,m=ψlϕm), as


The charge density Cl is used to find the solution of the Poisson equation, which provides an update for the electric force F(r,t). The latter governs the trajectories evolving the particles in the next time interval (τ,τ+Δt). Between the steps of solving the Poisson equation the electric field is frozen so that event biasing can be applied. Assume that at time τ the particles emerge with weights Wn. Due to the event biasing the behavior of the biased particles differs from that of the EMC particles. The distribution function of the biased particles fl.mnum obtained from the above formula is entirely different from fl,m. Nevertheless, as seen from (15), any biasing does not change the values of the physical averages. The Boltzmann distribution function is recovered by using the weights Wn as


Accordingly, the correct F is provided by the Poisson equation. As the evolution is Markovian, fl,m presents the initial condition for the next time step. Numerical particles, having distribution fl.mnum and weights Wn present a biased initial condition for this step. The initial weight will be updated in the time interval (τ,τ+Δt) by the weight factors according to the chosen biased evolution. It follows that the particle weights survive between the successive iteration steps, which completes the proof of the self-consistent biasing scheme.

4.2. Event biasing methods employed

Three different event-biasing techniques have been employed in this work. The chosen biasing techniques aim at increasing the number of numerical particles in the channel, but keep the total particle number in the device constant (for example, equal to 105). Particles, which enrich the high energy domain of the distribution on the expense of obtaining weights less than unity readily overcome the source potential barrier. Particle number in the low energy domain is less than the conventional ones, with higher weights and remain longer in the S/D regions. The methods are discussed in the following.


Figure 5.

Weighting scheme used in the temperature-biasing method.


Figure 6.

a) Weight distribution for T = 450K, and (b) Biasing decreases numerical particle number in the source/drain regions while increases the number in the channel region. Here 3n stands for numerical particles under bias = 3.

4.2.1. Biasing the initial/boundary temperature

Denoting the equilibrium distribution,


where, ε is the individual carrier energy and ε¯=1.5kBT, one chooses a biased distribution


which corresponds to higher temperature Tb=bias×T (bias > 1). Having higher kinetic energy, the numerical particles readily overcome the source potential barrier and enrich the statistics in the channel. The weight distribution is governed by the formula


The biasing scheme is illustrated in Figure 5. Increasing Tb increases the spread of the weight (and hence the distribution function) further away from unity which may lead to an increased variance of the physical averages, obtained by the mean of heavy and light particles. Thus finding the appropriate bias is a matter of compromise between the need for more particles in the channel (high temperature) and keeping the spread of the weight low (low temperature). The particle weight distribution for a particular choice of bias = 1.5 (corresponding to a temperature of 450K) is shown in Figure 6(a). That biasing increases the numerical particle number in the channel region by decreasing the same in the S/D regions in illustrated in Figure 6(b) with bias=3. Also noticeable is the fact that different values of bias do not change the actual number of physical electrons throughout the device region.

4.2.2. Particle Split

Particle weight is controlled by choosing a desired weight w1 of the numerical particles with kinetic energy below given level ε1 and weight w2 with kinetic energy above ε1. fb is obtained from feq as follows:


w2 is obtained as a function of w1 and ε1 from the condition for normalization of fb:


A choice of w1 > 1 effectively reduces the number of particles below ε1 as compared to the unbiased case. In both the above cases of biasing heavy particles which enter the channel perturb the statistics accumulated by a set of light weight particles (Figure 7). It is thus desirable to apply the technique of particle splitting in parallel to the temperature biasing in order to minimize the spread of the weight.


Figure 7.

Particle split method and the distribution of particles in the device active region.

4.2.3. Biasing Phonon Scattering (e-a)

Artificial carrier heating can be achieved by biasing the phonon scattering rates. For a given scattering mechanism, the probability for phonon absorption is increased at the expense of phonon emission, controlled by a parameter w1,


If in the course of the simulation, a phonon absorption is selected, the particle weight is updated by a multiplication with λabs/λabsb, otherwise with λem/λemb. The distribution of the flight time is not affected, because the sum of emission and absorption rate is not changed.

4.3. Results from biasing experiments

The MOSFET chosen for the simulation experiments has gate length of 15 nm, channel doping of 2×1019cm-3, and oxide thickness of 0.8 nm. Similar device has already been fabricated by Intel [43]. The applied potential VG = 0.375V, VD = 0.1V correspond to a subthreshold regime at lattice temperature T = 300 K.

Figure 8 shows the standard deviation in electron number as a function of evolution time in a certain cell within the channel region. One can clearly see the improvement achieved from the use of event biasing techniques in the simulation. The standard deviation is least for method (a) where the initial/boundary temperature of the electrons is biased (at T = 450 K).


Figure 8.

Enhancement of channel statistics: reduction of standard deviation in number of electrons in a particular cell in the channel region.


Figure 9.

Biasing recovers precisely the self-consistent (a) average sheet density and (b) average kinetic energy of the electrons.

Regarding the validation of the biasing techniques, first, the consistency of the biasing techniques in the thermodynamic limit of a very large number (105) of simulated particles is investigated. Both Boltzmann (conventional) and biased stochastic processes must give the same evolution of the physical averages. Figure 9 shows that the biased experiments recover precisely the physical averages of electron sheet density and electron energy along the channel of the device.

Secondly, the convergence of the cumulative averages for the channel and terminal currents obtained from the velocity and particle counting, respectively, is investigated. Biasing the phonon scattering rates (e-a) is applied in the half of the source region near the barrier in a 4 nm depth. Figure 10 (top panel) shows the biased channel current as compared to the EMC (conventional) result for 30 ps evolution time. The 5% error region (straight lines) around the mean value is entered 2.5 ps earlier and the convergence is better. The channel current from the boundary temperature biasing (T = 450K) is shown in the bottom panel of Figure 10. The temperature-biased curve shows a superior behavior.


Figure 10.

Comparison of channel currents obtained from (1) biased e-a rates (top panel), and (2) biased boundary distribution (bottom panel) methods from velocity consideration.

The corresponding terminal currents (shown in Figure 11) are much noisier and show long-time correlations due to the inter-particle interactions. The e-a biased curve is very unstable and enters the 5% error region in Figure 11 (top panel), after 15 ps evolution. One can associate this behavior with numerical error. The poor statistics is due to the appearance of very heavy (W>2) particles in the source as can be seen in Figure 12(a). To check this, it is sufficient to apply the conventional particle splitting coupled with the e-a biasing. The result is presented by the dotted curve in Figure 11. The behavior is significantly improved with the expense of a 30% increase of the simulated particles (the corresponding weight distribution is shown in Figure 12(b)). The terminal current corresponding to the biasing of the temperature of the injected particles again shows a superior behavior. This is due to an improved weight control. The weight, determined during the injection remains constant in the evolution. Its maximal value for T = 450 K is exactly 1.5 as depicted in Figure 6(a). Furthermore the probability for interaction with the impurities, which dominates the source/drain regions, drops for the majority of the particles due to their high energies. The conventional splitting technique cannot achieve such superiority. The terminal current from particle split technique is shown by the dotted curve on Figure 11 (bottom panel). The behavior of the curve resembles the EMC counterpart. An improvement is expected if the w2 particles are additionally split, which recovers the conventional split technique.


Figure 11.

Terminal currents obtained from various methods by particle counting.


Figure 12.

Numerical particle weight distribution in (a) e-a biasing, and (b) e-a/split biasing.

5. Conclusion

In conclusion, three event biasing techniques for particle-based Monte Carlo simulations of semiconductor devices have been derived in presence of both initial and boundary conditions and generalized for self-consistent simulations. All these approaches are confirmed and validated by simulations of a MOSFET with a channel length of 15 nm in the subthreshold regime of operation. A bias technique, particularly useful for small devices, is obtained by injection of hot carriers from the boundaries. The coupling with the Poisson equation requires a precise statistics in the source/drain regions. It is shown that a combination of event biasing and population control approaches is advantageous for this purpose.


Shaikh Ahmed would like to thank Sharnali Islam for her assistance in editing the text. This work was partially supported by the National Science Foundation (Grant No. ECCS-1102192). Computational resource/time supported by the ORAU/ORNL High-Performance Computing Grant 2009 is acknowledged. Mihail Nedjalkov would like to acknowledge the Austrian Science Fund (FWF) P21685-N22 and the Bulgarian Science Fund DTK 02/44.


1 - Iafrate G. Introduction. In: Ferry DK., Jacoboni C. (ed.) Quantum Transport in Semiconductors. Plenum Press: New York and London; 1992. pxv-xxi.
2 - Ferry DK. Semiconductor Transport. Taylor & Francis: London, UK; 2000.
3 - Lewis EE., Miller, WF. Computational Methods of Neutron Transport, American Nuclear Society, Inc.: La Grange Park, Illinois; 1993.
4 - Majumdar A. Development of a Multiple Perturbation Monte Carlo Method for Eigenvalue Problems and Implementation on Parallel Processors. Ph.D. dissertation: The University of Michigan; 1996.
5 - Duderstadt JJ., Martin WR. Transport Theory. John Wiley & Sons: New York; 1979.
6 - Hammersley JM. et al. Monte Carlo Methods. John Wiley & Sons: New York; 1964.
7 - Lux I., Koblinger L. Monte Carlo Particle Transport Methods: Neutron and Photon Calculations. CRC Press: Florida; 1991.
8 - Rubinstein RY. Simulation and the Monte Carlo Method. John Wiley & Sons: New York; 1981.
9 - Spanier J., Gelbard EM. Monte Carlo Principles and Neutron Transport Problems. Addison-Wesley: Reading, Massachusetts; 1969.
10 - Kosina H., Nedjalkov M., Selberherr S. Theory of the Monte Carlo method for semiconductor device simulation. IEEE Transactions on Electron Devices 2000; 47 1898–1908.
11 - Fjeldly TA., Shur MS. Simulation and Modeling of Compound Semiconductor Devices. In: Shur MS. (ed.) Compound Semiconductor Electronics: The Age of Maturity. World Scientific; 1996. 317–364.
12 - Vasileska D., Mamaluy D., Khan H. R., Raleva K., Goodnick SM. Semiconductor Device Modeling. Journal of Computational and Theoretical Nanoscience 2008; 5 1–32.
13 - Ferry DK., Goodnick SM. Transport in Nanostructures. Cambridge University Press; 1997.
14 - Wigner E. On the quantum correction for thermodynamic equilibrium. Phys. Rev 1932; 40 749–759.
15 - Feynman R., Kleinert H. Effective classical partition functions. Phys. Rev. A 1986; 34 5080–5084.
16 - Datta S. Quantum Transport: Atom to Transistor. Cambridge University Press; 2005.
17 - Lake R., Klimeck G., Bowen RC., Jovanovic D. Single and multiband modeling of quantum electron transport through layered semiconductor devices J. Appl. Phys. 1997; 81 7845.
18 - Iafrate GJ., Grubin HL., Ferry DK. Utilization of quantum distribution functions for ultra-submicron device transport. Journal de Physique 1981; 42 (Colloq. 7) 307–312.
19 - Ferry DK., Zhou JR. Form of the quantum potential for use in hydrodynamic equations for semiconductor device modeling. Phys. Rev. B 1993; 48 7944–7950.
20 - Ferry DK. Effective potential and the onset of quantization in ultrasmall MOSFETs. Superlattices and Microstructures 2000; 28 419–423.
21 - Ringhofer C., Gardner CL. Smooth QHD model simulation of the resonant tunneling diodes. VLSI Design 1998; 8 143–146.
22 - Vasileska D., Khan H., Ahmed SS. Modeling Coulomb effects in nanoscale devices. Journal of Computational and Theoretical Nanoscience 2008; 5(9) 1793–1827.
23 - Broglie L. De. Sur la possibilité de relier les phénomènes d'interference et de diffraction á la théorie des quanta de luminère. C. R. Acad. Sci. Paris 1926; 183 447–448.
24 - Madelung E. Quantatheorie in hydrodynamischer form. Z. Phys. 1926; 40 322–326.
25 - Bohm D. A suggested interpretation of the quantum theory in terms of hidden variables I and II. Phys. Rev. 1952; 85 166–193.
26 - Fischetti M., Laux S. Monte Carlo study of electron transport in silicon inversion layers. Physical Review B 1993; 48 2244–22743.
27 - Kurosawa. J. Phys. Soc. Jpn. 1966; 21 424.
28 - Fawcett W., Boardman DA., Swain S. J. Phys. Chem. Solids 1970; 31 1963.
29 - Jacoboni C., Reggiani L. The Monte Carlo Method for the Solution of Charge Transport in Semiconductors with Applications to Covalent Materials. Rev. Modern Phys 1983; 55 645–705.
30 - Price PJ. Semiconductors and Semimetals 1979; 14 249–308.
31 - Gross WJ. Three-dimensional particle-based simulations of deep sub-micron MOSFET devices. Ph. D. dissertation: Arizona State University; 1999.
32 - Lundstrom M. Fundamentals of Carrier Transport. Cambridge University Press; 2000.
33 - Tomizawa K. Numerical Simulation of Submicron Semiconductor Devices. Artech House: Boston; 1993.
34 - Kunikuyo T., Takenaka M., Kamakura Y., Yamaji M., Mizuno H., Morifuji M. A Monte Carlo simulation of anisotropic electron transport in silicon using full band structure and anisotropic impact-ionization model. Journal of Applied Physics 1994; 75 297–312.
35 - Ahmed SS. Quantum and Coulomb Effects in Nanoscale Devices. Ph.D. dissertation: Arizona State University; 2005.
36 - Ferry DK. First-Order Optical and Intervalley Scattering in Semiconductors. Phys. Rev. B 1976; 14 1605–1609.
37 - Gross WJ., Vasileska D., Ferry DK. 3D Simulations of Ultra-Small MOSFETs with Real-Space Treatment of the Electron-Electron and Electron-Ion Interactions. VLSI Design 2000; 10 437–452.
38 - Ahmed SS., Ringhofer C., Vasileska D. Parameter-Free Effective Potential Method for Use in Particle-Based Device Simulations. IEEE Transactions on Nanotechnology 2005; 4 465–471.
39 - Wordelman C., Kwan T., Snell C. Comparison of statistical enhancement methods for Monte Carlo semiconductor simulations. IEEE Transactions on Computer-Aided Design 1998; 17(12) 1230–1235.
40 - Rota L., Jacoboni C., Poli P. Weighted Ensemble Monte-Carlo. Solid-State Electronics 1989; 32(12) 1417–1421.
41 - Kosina H., Nedjalkov M., Selberherr S. The stationary Monte-Carlo method for device simulaton–Part I and Part II. Journal of Applied Physics 2003; 93(6) 3553–3571.
42 - Jacoboni C. Principles of Quantum Transport. In: Ferry D., Jacoboni C. (ed.) Quantum Transport in Semiconductors. Plenum Press: New York and London; 1992. p1–15.
43 - Jacoboni C. A New Approach to Monte Carlo Simulation. PRIEDM, IEEE Electron Devices Society 1989; 469–472.
44 - Chau R., Boyanov B., Doyle B., Doczy M., Datta S., Hareland S., Jin B., Kavalieros J., Metz M. Silicon nano-transistors for logic applications. Proceedings of the 4th Int. Symp. on Nanostructures and Mesoscopic Systems, 17-21 February 2003, Tempe Mission Palms Hotel, Temepe, Arizona, USA.