Monte Carlo methods for electron transport

The Boltzmann transport equation model has been the main tool used in the analysis of transport in semiconductors. The BTE equation is given by:

${\displaystyle {\frac {\partial f}{\partial t}}+{\frac {1}{\hbar }}\nabla _{k}E(k)\nabla _{r}f+{\frac {qF(r)}{\hbar }}\nabla _{k}f=\left[{\frac {\partial f}{\partial t}}\right]_{\mathrm {collision} }}$

${\displaystyle v={\frac {1}{\hbar }}\nabla _{k}E(k)}$

The distribution function, f, is a dimensionless function which is used to extract all observable of interest and gives a full depiction of electron distribution in both real and k-space. Further, it physically represents the probability of particle occupation of energy k at position r and time t. In addition, due to being a seven-dimensional integro-differential equation (six dimensions in the phase space and one in time) the solution to the BTE is cumbersome and can be solved in closed analytical form under very special restrictions. Numerically, solution to the BTE is employed using either a deterministic method or a stochastic method. Deterministic method solution is based on a grid-based numerical method such as the spherical harmonics approach, whereas the Monte Carlo is the stochastic approach used to solve the BTE.

The semiclassical Monte Carlo method is a statistical method used to yield exact solution to the Boltzmann transport equation which includes complex band structure and scattering processes. This approach is semiclassical for the reason that scattering mechanisms are treated quantum mechanically using the Fermi's Golden Rule, whereas the transport between scattering events is treated using the classical particle notion. The Monte Carlo model in essence tracks the particle trajectory at each free flight and chooses a corresponding scattering mechanism stochastically. Two of the great advantages of semiclassical Monte Carlo are its capability to provide accurate quantum mechanical treatment of various distinct scattering mechanisms within the scattering terms, and the absence of assumption about the form of carrier distribution in energy or k-space. The semiclassical equation describing the motion of an electron is

${\displaystyle {\frac {dr}{dt}}={\frac {1}{\hbar }}\nabla _{k}E(k)}$

${\displaystyle {\frac {dk}{dt}}={\frac {qF(r)}{\hbar }}}$

where F is the electric field, E(k) is the energy dispersion relation, and k is the momentum wave vector. To solve the above equation, one needs strong knowledge of the band structure (E(k)). The E(k) relation describes how the particle moves inside the device, in addition to depicting useful information necessary for transport such as the density of states (DOS) and the particle velocity. A Full-band E(K) relation can be obtained using the semi-empirical pseudopotential method.[1]

Both drift diffusion (DD) and the hydrodynamic (HD) models can be derived from the moments of the Boltzmann transport equation (BTE) using simplified approximation valid for long channel devices. The DD scheme is the most classical approach and usually solves the Poisson equation and the continuity equations for carriers considering the drift and diffusion components. In this approach, the charge transit time is assumed to be very large in comparison to the energy relaxation time.[2] On the other hand, the HD method solves the DD scheme with the energy balance equations obtained from the moments of BTE.[3][4] Thus, one may capture and calculate physical details such as carrier heating and the velocity overshoot effect. Needless to say, an accurate discretization method is required in HD simulation, since the governing equations are strongly coupled and one has to deal with larger number of variables compared to the DD scheme.

Average carrier velocity for an 80nm nmos comparing the various semiclassical simulation model (a) Vds= 0.3V (b) Vds= 0.6V

The accuracy of semiclassical models are compared based on the BTE by investigating how they treat the classical velocity overshoot problem, a key short channel effect (SCE) in transistor structures. Essentially, velocity overshoot is a nonlocal effects of scaled devices, which is related to the experimentally observed increase in current drive and transconductance.[5] As the channel length becomes smaller, the velocity is no longer saturated in the high field region, but it overshoots the predicted saturation velocity. The cause of this phenomenon is that the carrier transit time becomes comparable to the energy relaxation time, and therefore the mobile carriers do not have enough time to reach equilibrium with the applied electric field by scattering in the short channel devices.[6] The summary of simulation results (Illinois Tool: MOCA) with DD and HD model is shown in figure beside. In the figure (a), the case when the field is not high enough to cause the velocity overshoot effect in the whole channel region is shown. Note that at such limit, the data from the DD model fit well to the MC model in the non-overshoot region, but the HD model overestimate the velocity in that region. The velocity overshoot is observed only near the drain junction in the MC data and the HD model fits well in that region. From the MC data, it can be noticed that the velocity overshoot effect is abrupt in the high-field region, which is not properly included in the HD model. For high field conditions as shown in the figure (b) the velocity overshoot effect almost all over the channel and the HD results and the MC results are very close in the channel region.

Band structure describes the relationship between energy(E) and wave vector(k). The band structure is used to compute the movement of carriers under the action of the electric field, scattering rate, and final state after the collision. Silicon band structure and its Brillouin zone are shown in figure below, but there is no analytical expression which satisfies entire Brillouin zone. By using some approximation, there are two analytical models for band structure, namely the parabolic and the non-parabolic modes.

For the concept of band structure, parabolic energy bands are generally assumed for simplicity. Electrons reside, at least when close to equilibrium, close to the minima of the E(k) relation. Then the E(k) relation can be extended in a Taylor series as

${\displaystyle E(k)=E(0)+\left.{\frac {\partial E(k)}{\partial k}}\right|_{\mathrm {k=0} }\cdot k+{\frac {1}{2}}{\frac {\partial ^{2}E(k)}{\partial k^{2}}}\cdot k^{2}}$

Because the first derivative vanishes at the band minimum, so the gradient of E(k) is zero at k = 0. Thus,

${\displaystyle E(k)={\frac {\hbar ^{2}k^{2}}{2m^{*}}}}$

which yields the definition of the effective mass tensor

${\displaystyle {\frac {1}{m^{*}}}={\frac {1}{\hbar ^{2}}}{\frac {\partial ^{2}E(k)}{\partial k^{2}}}}$

This expression is true for semiconductor which has isotropic effective mass, for instance GaAs. In case of silicon, conduction band minima does not lie at k = 0 and the effective mass depends on the crystallographic orientation of the minimum as

${\displaystyle E(k)={\frac {\hbar ^{2}}{2}}\left({\frac {k_{l}^{2}}{m_{l}^{*}}}+{\frac {2k_{t}^{2}}{m_{t}^{*}}}\right)}$

where ${\displaystyle m_{l}^{*},m_{t}^{*}}$ describe longitudinal and transverse effective mass, respectively.

For higher applied fields, carriers reside above the minimum and the dispersion relation, E(k), does not satisfy the simple parabolic expression described above. This non-parabolicity is generally described by

${\displaystyle E(1+\alpha E)={\frac {\hbar ^{2}k^{2}}{2m^{*}}}}$

where ${\displaystyle \alpha }$ is a coefficient of non-parabolicity given by

${\displaystyle \alpha ={\frac {(1-m^{*}/m_{0})^{2}}{E_{g}}}}$

where ${\displaystyle m_{0}}$ is the electron mass in vacuum, and Eg is the energy gap.[7]

For many applications, non-parabolic band structure provides reasonable approximation. However, in case of very high field transport, which requires the better physical model of the full band structure. For full band approach, numerically generated table of E(k) is used. Full band approach for Monte Carlo simulation was first used by Karl Hess at the University of Illinois at Urbana-Champaign. This approach is based on empirical pseudopotential method suggested by Cohen and Bergstresser [18]. Full band approach is computationally expensive, however, following the advancement of the computational power, it can be used as a more general approach.[8]

For this type of simulation, one carrier is injected and the motion is tracked in the domain, until it exits through contact. Another carrier is then injected and the process repeated to simulate an ensemble of trajectories. This approach is mostly useful to study bulk properties, like the steady state drift velocity as a function of field.

Instead of single carrier, a large ensemble of carriers is simulated at the same time. This procedure is obviously a good candidate for super-computation, since one may apply parallelization and vectorization. Also, it is now possible to perform ensemble averages directly. This approach is suitable for transient simulations.

This method couples the ensemble Monte Carlo procedure to Poisson's equation, and is the most suitable for device simulation. Typically, Poisson's equation is solved at fixed intervals to update the internal field, to reflect the internal redistribution of charge, due to the movement of carriers.

The probability that the electron will suffer its next collision during dt around t is given by

${\displaystyle p(t)\,dt=P[k(t)]\exp[-\int _{0}^{t}P[k(t')]\,dt']\,dt}$

where P[k(t)]dt is the probability that an electron in the state k suffers a collision during the time dt. Because of the complexity of the integral at the exponent, it is impractical to generate stochastic free flights with the distribution of the equation above. In order to overcome this difficulty, people use a fictitious “self-scattering” scheme. By doing this, the total scattering rate, including this self-scattering, is constant and equal to, say, ${\displaystyle \Gamma }$. By random selection, if self-scattering is selected, k′ after the collision is the same as k and the carrier continues its flight without perturbation. Introducing a constant ${\displaystyle P(k)=\tau _{0}^{-1}}$, the above equation reduces to

${\displaystyle p(t)={\frac {1}{\tau _{0}}}\exp(-t/\tau _{0}).}$

Random numbers r can be used very simply to generate stochastic free flights, which duration will then be given by ${\displaystyle t_{r}=-\tau _{0}\ln(r)}$. The computer time used for self-scattering is more than compensated for by the simplification of the calculation of the free-flight duration.[9] To enhance the speed of free flight time calculation, several schemes such as “Constant Technique”, and “Piecewise Technique” are used to minimize the self-scattering events.

Important charge transport properties of semiconductor devices such as the deviance from Ohm's law and the saturation of carriers mobility are a direct consequence of scattering mechanisms. It is thus of great importance for a semiconductor device simulation to capture the physics of such mechanisms. The semiconductor Monte Carlo simulation, in this scope, is a very powerful tool for the ease and the precision with which an almost exhaustive array of scattering mechanisms can be included. The duration of the free flights is determined from the scattering rates. At the end of each flight, the appropriate scattering mechanism must be chosen in order to determine the final energy of the scattered carrier, or equivalently, its new momentum and scattering angle. In this sense, one will distinguish two broad types of scattering mechanisms which naturally derive form the classic kinetic theory of collision between two bodies:

Elastic scattering, where the energy of the particle is conserved after being scattered. Elastic scattering will hence only change the direction of the particle's momentum. Impurity scattering and surface scattering are, with a fair approximation, two good examples of elastic scattering processes.

Inelastic scattering, where energy is transferred between the scattered particle and the scattering center. Electronphonon interactions are essentially inelastic since a phonon of definite energy is either emitted or absorbed by the scattered particle. Before characterizing scattering mechanisms in greater mathematical details, it is important to note that when running semiconductor Monte Carlo simulations, one has to deal mainly with the following types of scattering events:[9]

Acoustic Phonon: The charge carrier exchanges energy with an acoustic mode of the vibration of atoms in the crystal lattice. Acoustic Phonons mainly arise from thermal excitation of the crystal lattice.

Polar Optical: The charge carrier exchanges energy with one of the polar optical modes of the crystal lattice. These modes are not present in covalent semiconductors. Optical phonons arise from the vibration against each other of atoms of different types when there is more than one atom in the smallest unit cell, and are usually excited by light.

Non-Polar Optical: Energy is exchanged with an optical mode. Non-polar optical phonons must generally be considered in covalent semiconductors and the L-valley of GaAs.

Equivalent Intervalley Phonon: Due to the interaction with a phonon, the charge carrier transitions from initial states to final states which belong to different but equivalent valleys. Typically, this type of scattering mechanism describes the transition of an electron from one X-valley to another X-valley, or from one L-valley to another L-valley.[10]

Non Equivalent Intervalley Phonon: Involves the transition of a charge carrier between valleys of different types.

Piezoelectric Phonon: For low temperatures.

Ionized Impurity: Reflects the deviation of a particle from it ballistic trajectory due to Coulomb interaction with an ionized impurity in the crystal lattice. Because the mass of an electron is relatively small in comparison to the one of an impurity, the Coulomb cross section decreases rapidly with the difference of the modulus of momentum between the initial and final state.[9] Therefore, impurity scattering events are mostly considered for intravalley scattering, intraband scattering and, to a minor extent, interband scattering.

Carrier-Carrier: (electron-electron, hole-hole and electron-hole interactions). When carrier concentration is high, this type of scattering reflects the electrostatic interaction between charge carriers. This problem becomes very quickly computationally intensive with an increasing number of particles in an ensemble simulation. In this scope, Particle-Particle–Particle-Mesh (P3M) algorithms, which distinguish short range and long range interaction of a particle with its surrounding charge gas, have proved efficient in including carrier-carrier interaction in the semiconductor Monte Carlo simulation.[11] Very often, the charge of the carriers is assigned to a grid using a Cloud-in-Cell method, where part of the charge of a given particle is assigned to a given number of closest grid points with a certain weight factor.

Plasmon: Reflects the effect of the collective oscillation of the charge carriers on a given particle.

A computationally efficient approach to including scattering in Monte Carlo simulation consists in storing the scattering rates of the individual mechanisms in tables. Given the different scattering rates for a precise particle state, one may then randomly select the scattering process at the end of the free flight. These scattering rates are very often derived using the Born approximation, in which a scattering event is merely a transition between two momentum states of the carrier involved. As discussed in section II-I, the quantum many-body problem arising from the interaction of a carrier with its surrounding environment (phonons, electrons, holes, plasmons, impurities,...) can be reduced to a two-body problem using the quasiparticle approximation, which separates the carrier of interest from the rest of the crystal.[9] Within these approximations, Fermi's Golden Rule gives, to the first order, the transition probability per unit time for a scattering mechanism from a state ${\displaystyle |k\rangle }$ to a state ${\displaystyle |k'\rangle }$:

${\displaystyle S(k,k')={\frac {2\pi }{\hbar }}\left|\langle k|H'|k'\rangle \right|^{2}\cdot \delta (E-E')}$

where H' is the perturbation Hamiltonian representing the collision and E and E′ are respectively the initial and final energies of the system constituted of both the carrier and the electron and phonon gas. The Dirac ${\displaystyle \delta }$-function stands for the conservation of energy. In addition, the term ${\displaystyle \langle k|H'|k'\rangle }$, generally referred to as the matrix element, mathematically represents an inner product of the initial and final wave functions of the carrier:[12]

${\displaystyle \langle k|H'|k'\rangle ={\frac {1}{Vol}}\int _{\mathrm {Vol} }\psi _{k}(r)H'\psi _{k'}^{*}(r)\,dr}$

In a crystal lattice, the wavefunctions ${\displaystyle \psi _{k}(r)}$ and ${\displaystyle \psi _{k'}(r)}$ are simply Bloch waves. When it is possible, analytic expression of the Matrix elements are commonly found by Fourier expanding the Hamiltonian H', as in the case of Impurity scattering [13] or acoustic phonon scattering.[14] In the important case of a transition from an energy state E to an energy state E' due to a phonon of wave vector q and frequency ${\displaystyle \omega _{q}}$, the energy and momentum change is:

${\displaystyle E'-E=E(k')-E(k)\pm \hbar \omega _{q}\,}$

${\displaystyle k'-k\pm q={\begin{cases}0&{\text{ }}\\R&{\text{Umklapp-process}}\end{cases}}}$

where R is a reciprocal lattice vector. Umklapp processes (or U-processes) change the momentum of the particle after scattering and are therefore limiting the conduction in semiconductor crystals. Physically, U-processes occur when the final momentum of the particle points out of the first Brillouin zone. Once one knows the scattering probability per unit time from a state k to a state k', it is interesting to determine the scattering rate for a given scattering process. The scattering rate gives the probability per unit time to scatter from a state k to any other state in the reciprocal space. Therefore, the scattering rate is

${\displaystyle \lambda (k)=\sum _{k'}S(k,k')}$

which can be readily used to determine the free flight time and the scattering process as discussed in section 3-3. It is important to note that this scattering rate will be dependent on the band structure of the material (the dependence arises from the matrix elements).

At the end of a free flight, a scattering mode and angle must be randomly chosen. In order to determine the scattering mechanism, one has to consider all the scattering rates ${\displaystyle \lambda _{1},\lambda _{2},...,\lambda _{n}}$ of the mechanisms relevant to the simulation as well as the total scattering rate at the time of scattering ${\displaystyle \lambda _{tot}(t_{sc})=\sum _{i}\lambda _{i}.}$ Selecting a scattering mechanism then simply results in generating a uniformly distributed random number 0 < r < 1 and referring to the following rules

${\displaystyle {\begin{aligned}r&<{\frac {\lambda _{1}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}1\\r&<{\frac {\lambda _{1}+\lambda _{2}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}2\\&{}\ \vdots \\r&<{\frac {\sum _{i=0}^{n}\lambda _{i}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}n\end{aligned}}}$

A computationally efficient approach to selecting the scattering mechanism consists in adding a “void” scattering mechanism so that ${\displaystyle \lambda _{\mathrm {tot} }}$ remains constant over time. If a particle is scattered according to this mechanism, it will keep its ballistic trajectory after scattering takes place. In order to choose a new trajectory, one must first derive the energy (or momentum) of the particle after scattering

${\displaystyle E(k')=E(k)\pm \hbar \omega _{q}\pm \Delta E_{C}\,}$

where the term ${\displaystyle \hbar \omega _{q}}$ accounts for phonon emission or absorption and the term ${\displaystyle \Delta E_{C}}$ is non-null for inter-valley scattering. The final energy (and the band structure) directly yield the modulus of the new momentum k'. At this point one only needs to choose a new direction (or angle) for the scattered particle. In some simple cases as phonon scattering and a parabolic dispersion relation, the scattering angle is random and evenly distributed on the sphere of radius k'. Using spherical coordinates, the process of choosing the angle is equivalent to randomly picking two angles ${\displaystyle \theta }$and ${\displaystyle \psi }$. If the angle is distributed with a distribution ${\displaystyle p(\theta ,\psi )}$, then for a uniform distribution of angles, the probability to pick a point of the sphere is

${\displaystyle p(\theta ,\psi )\,d\theta d\psi ={\frac {\sin \theta \,d\theta \,d\psi }{4\pi }}}$

It is possible, in this case, to separate the two variables. Integrating over ${\displaystyle \psi }$ then over ${\displaystyle \theta }$, one finds

${\displaystyle p(\theta )={\frac {\sin \theta }{2}}}$

${\displaystyle p(\psi )={\frac {1}{2\pi }}}$

The two spherical angles can then be chosen, in the uniform case, by generating two random numbers 0 < r₁, r₂ < 1 such that

${\displaystyle r_{1}=\int _{0}^{\psi }p(\psi ')\,d\psi '={\frac {\psi }{2\pi }}}$

${\displaystyle r_{2}=\int _{0}^{\theta }p(\theta ')\,d\theta '={\frac {1-\cos \theta }{2}}}$

The Wigner transport equation forms the bases for the Wigner-based quantum correction.

${\displaystyle {\frac {\partial f}{\partial t}}+r\cdot \nabla _{r}f-{\frac {1}{\hbar }}\nabla _{r}V\cdot \nabla _{k}f+\sum _{\alpha =1}^{\infty }{\frac {(-1)^{\alpha +1}}{\hbar 4^{\alpha }(2\alpha +1)!}}\times (\nabla _{r}\nabla _{k})^{2\alpha +1}Vf=\left({\frac {\partial f}{\partial t}}\right)_{c}}$

where, k is the crystal momentum, V is the classical potential, the term on the RHS is the effect of collision, the fourth term on the LHS represents non-local quantum mechanical effects. The standard Boltzmann Transport Equation is obtained when the non-local terms on the LHS disappear in the limit of slow spatial variations. The simplified (for ${\displaystyle \alpha =0}$) quantum corrected BTE then becomes

${\displaystyle {\frac {\partial f}{\partial t}}+r\cdot \nabla _{r}f-{\frac {1}{\hbar }}\nabla _{r}V\cdot \nabla _{k}f=\left({\frac {\partial f}{\partial t}}\right)_{c}}$

where the quantum potential is contained in the term ${\displaystyle V_{\omega }}$ (must be an error: ${\displaystyle V_{\omega }}$ was never mentioned).

This method for quantum correction was developed by Feynman and Hibbs in 1965. In this method the effective potential is derived by calculating the contribution to the path integral of a particle's quantum fluctuations around its classical path. This calculation is undertaken by a variational method using a trial potential to first order. The effective classical potential in the average point on each path then becomes

${\displaystyle V_{\mathrm {eff} }(x)={\frac {1}{\sqrt {2\pi a}}}\int _{-\infty }^{\infty }V(x')e^{-{\frac {(x'-x)^{2}}{2a^{2}}}}dx'}$

${\displaystyle a^{2}={\frac {\hbar ^{2}}{12m^{*}k_{B}T}}}$

This approach involves periodical solving of a Schrödinger equation in a simulation with the input being the self-consistent electrostatic potential. The exact energy levels and wavefunctions relating to the electrostatic potential solution are employed to calculate the quantum potential. The quantum correction obtained on the bases of this method can be visualised by the following equation

${\displaystyle V_{\mathrm {schr} }(z)=-k_{B}T\cdot \log(n_{q}(z))-V_{p}(z)+V_{0}}$

where V_schr is the quantum correction potential, z is the direction perpendicular to the interface, n_q is the quantum density from the Schrödinger equation which is equivalent to the converged Monte Carlo concentration, V_p is the potential from the Poisson solution, V₀ is the arbitrary reference potential far away from the quantum region such that the correction goes to null in the region of semi-classical behavior. Even though the above-mentioned potentials for quantum correction differ in their method of calculation and their basic assumptions, yet when it comes to their inclusion into Monte Carlo simulation they are all incorporated the same way.