## 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,

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]):

Initialization.

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 [31]

where *m*_{i} is the effective mass of the particle in the *i*^{th} (*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.

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 (P^{3}M) 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].

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/SiO_{2} 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/SiO_{2} 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, *n*_{absorbed} is the number of particles that are absorbed by the contact, *n*_{injected} is the number of particles that have been injected at the contact, *E*_{y} 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

where *n*_{net} 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 *v*_{d} is the average electron drift velocity and *n* is the carrier concentration. If there are a total of *N* particles in a differential volume, *dA*, is

or,

where *v*_{x}(n) is the velocity along the channel of the *n*^{th} 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 [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* such as carrier density and velocity given by

Where

where *S* is the usual scattering rate from lattice imperfections,

and

If both, initial

While *D*.

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

It takes values determined by the second half with a probability given by the product in the first half in equation (15). *P*_{0} and *P* are used to construct numerical trajectories: (i) *Q* provided that *W*_{2} in front of *A*, called *weight*, is a product of weight factors *N* realizations of the r.v., calculated over *N* trajectories

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 *D* and 0 (zero) otherwise. The choice of

The initial probability *N*_{i} and the total number *N*_{J} of the injected particles into the device. The ratio *N*_{i} and *N*_{J} respectively, and can be eliminated by the choice

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*-th particle is inside or outside *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

The Boltzmann equation for Coulomb carriers becomes nonlinear via the interaction component *C*_{l} and the distribution function

The charge density *C*_{l} is used to find the solution of the Poisson equation, which provides an update for the electric force *frozen* so that event biasing can be applied. Assume that at time *W*_{n}. 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 *W*_{n} as

Accordingly, the correct F is provided by the Poisson equation. As the evolution is Markovian, *W*_{n} present a biased initial condition for this step. The initial weight will be updated in the time interval *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 10^{5}). 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,

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 *T*_{b} 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* 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 *w*1 of the numerical particles with kinetic energy below given level *ε*1 and weight *w*2 with kinetic energy above *ε*1. *f*^{b} is obtained from *f*_{eq} as follows:

*w*2 is obtained as a function of *w*1 and *ε*1 from the condition for normalization of *f*^{b}:

A choice of *w*1 > 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 *w*1,

If in the course of the simulation, a phonon absorption is selected, the particle weight is updated by a multiplication with

### 4.3. Results from biasing experiments

The MOSFET chosen for the simulation experiments has gate length of 15 nm, channel doping of ^{-3}, and oxide thickness of 0.8 nm. Similar device has already been fabricated by Intel [43]. The applied potential *V*_{G} = 0.375V, *V*_{D} = 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 (10^{5}) 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 *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 *w*2 particles are additionally split, which recovers the *conventional* split technique.

## 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.

## Acknowledgements

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.