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, , 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) 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 : (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 . 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 . 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 . 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 . 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 , Wigner functions , Feynman path integrals , and non-equilibrium Green’s functions (NEGF)  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 , 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 , 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 , and later developed by Bohm . 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 . 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  in 1966. Shortly afterwards the Malvern, UK, group  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  and by Peter J. Price  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. ):
Solution of the field (Poisson) equation: determine forces on electrons.
Simulate the electron dynamics (free-flight and scattering):
Accelerate the electrons.
Determine whether the electron experiences a scattering/collision.
If electron scatters, select the scattering mechanism.
Update the electron position, energy and wavevector.
Charge assignment to the numerical nodes/grids.
Compute measurable quantities such as average energy, average velocity, and mobility of the particle ensemble.
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 
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 , Tomizawa , and Kunikiyo et al.  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.  and .
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 . 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 . 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 .
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.
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. .
Quantum mechanical size-quantization effects have been accounted for via a parameter-free effective potential scheme . 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 : (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 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, 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.
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 . 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) , has been recently extended for the Single Particle Monte Carlo technique (stationary problem) . 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 of generic physical quantities a such as carrier density and velocity given by
Where and Equation (8) denotes the integration over the phase space and time , and 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 . The latter is determined by the initial condition in evolution problems  or, in the case of stationary transport, by the boundary conditions . The second equality in (8) follows from the relationship between an integral equation and its adjoint equation. It shows that the mean value is determined by and by the solution of the adjoint Boltzmann equation
where S is the usual scattering rate from lattice imperfections, is the total out-scattering rate, is the device domain indicator, which is discussed later, is the Heaviside function and the trajectories, initialized by , are formulated with the help of the electrical force F and the velocity v as
If both, initial and boundary conditions are taken into account, it can be shown that becomes
While is defined only at the initial time , the function is defined only at the device boundary Γ and for values of k such that the corresponding velocity inwards the domain, D. is the velocity component normal to Γ so that a velocity-weighted distribution drives the particle flux, injected into the device at times . in (13) governs both the transient and the stationary behavior of a device. The latter is established in the long time limit, provided that is time independent. Usually 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 :
Consider the second term in (14) augmented with the help of two probabilities and to become expectation  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). is evaluated according to the numerical Monte Carlo theory as follows. P0 and P are used to construct numerical trajectories: (i) selects the initial point of the trajectory. (ii) selects the next trajectory point Q provided that is given. The fraction W2 in front of A, called weight, is a product of weight factors , and evaluated at the corresponding points , selected by application of . The sample mean of N realizations of the r.v., calculated over N trajectories , estimates the mean value :
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 . Next, one establishes the link between (16) and the EMC technique, which is due to particular choice of the initial, , and transition, , densities. , 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 is then the domain indicator which takes values 1 (one) if the trajectory belongs to D and 0 (zero) otherwise. The choice of 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 for each estimator is obtained from and 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 for each of the estimators becomes Ni and NJ respectively, and can be eliminated by the choice and . 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. is the number of such particles at time , and 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 where the position of the biased particles is accounted in .
The Boltzmann equation for Coulomb carriers becomes nonlinear via the interaction component 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, . The particle system is evolved in time intervals fs. At the end of each time step, at say time , the charge density is calculated and assigned to the corresponding grid points. One uses the relation between Cl and the distribution function , which is estimated with the help of (18) by introducing a mesh in the wavevector space, , as
The charge density Cl is used to find the solution of the Poisson equation, which provides an update for the electric force . The latter governs the trajectories evolving the particles in the next time interval . 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 obtained from the above formula is entirely different from . 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, presents the initial condition for the next time step. Numerical particles, having distribution and weights Wn present a biased initial condition for this step. The initial weight will be updated in the time interval 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.
4.2.1. Biasing the initial/boundary temperature
Denoting the equilibrium distribution,
where, is the individual carrier energy and , one chooses a biased distribution
which corresponds to higher temperature (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 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.
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 , otherwise with . 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 cm-3, and oxide thickness of 0.8 nm. Similar device has already been fabricated by Intel . 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).
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.
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 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.
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.