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:
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
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
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
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
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.
To simulate an actual
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
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.
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
4.1. Theory of Event Biasing
The Ensemble Monte Carlo (EMC) technique is designed to evaluate averaged values of generic physical quantities
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
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,
A recursive replacement of equation (9) into itself gives rise to the von-Neumann expansion, where the solution
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.
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
The initial probability for each estimator is obtained from and respectively, with the help of two normalization factors: the number of initial carriers
Equation (18) accounts that only trajectories which belong to
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
The charge density
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
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
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 (
The biasing scheme is illustrated in Figure 5. Increasing
4.2.2. Particle Split
Particle weight is controlled by choosing a desired weight
A choice of
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
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
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
Regarding the validation of the biasing techniques,
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.