MOS/sup 2/: an efficient MOnte Carlo Simulator for MOS devices

IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN, VOL. 7, NO. 2, FEBRUARY 1988

259

MOS2: An Efficient A4Onte Carlo Simulator for MOS Devices ENRICO SANGIORGI, BRUNO RICCO,

MEMBER, IEEE, A N D

Abstract-An efficient Monte Carlo device simulator has been developed as a postprocessor of a two-dimensional numerical analyzer based on the drift-diffusion model. The MC package described in this work performs Monte Carlo analysis of real VLSI MOSFET’s in a minicomputer environment, overcoming some existing theoretical and practical problems. In particular the particle free-flight time distribution is obtained with a new algorithm leading to a CPU time saving of a t least one order of magnitude compared with the traditional approach. To describe r a r e electron configurations, such as the high energy tails of the distributions and the particle dynamics in the presence of large retarding fields, a multiple repetition scheme was implemented. Selected applications a r e presented to illustrate the simulator’s capabilities.

Average carrier density in the real space. Net number of carrier collected at an electrode. Thermal equilibrium number of optical phonons. 7) at a Probability that an electron in given time will be scattered between t and t + dt. Electron charge. Random number evenly distributed in the interval [ 0, 1 3. Carrier position in space. Domain of rare states. Time. Free flight time. Lattice temperature. Carrier velocity vector. Longitudinal sound velocity. Control volume around node i. Number of charge units of the impurity center. Non-parabolicity parameter. Impurity screening length. Scattering rate. Self-scattering rate. Total scattering rate including self-scattering . Root-mean-square deviation of the Si-Si02 interface from flatness. Effective dielectric constant at the Si-Si02 interface. Dielectric constant. Carrier mobility. Crystal density. Electric potential. Angular frequency Ensemble average at a given time.

(c,

LIST OF SYMBOLS

e D

z,7)

D(

Domain of common states. Carrier diffusion coefficient. Domain of all the possible states in the momentum and real space. Deformation potential constant for opticalphonon scattering with electron. Electric field Electron acoustic-phonon deformation potential. Reduced Planck constant. Normalization factor used to relate driftdiffusion and Monte Carlo quantities. Current. Carrier wave vector. Boltzmann constant. Correlation length of roughness. Electron effective mass. Electron conductivity effective mass. Carrier density in momentum and real space. Carrier density before scattering. Depletion charge per unit area ( at/cm2). Impurity concentration. Inversion charge per unit area ( at/cm2).

Manuscript received May 12, 1987; revised July 30, 1987. The review of this paper was arranged by R . W. Dutton. E. Sangiorgi was with the University of Bologna, Bologna, Italy. He is now with the Department of Physics, University of Udine, Udine, Italy. B. Ricco and F. Venturi are with the Department of Electronics, University of Bologna, Bologna, Italy. IEEE Log Number 87 17573.

FRANCO VENTURI

I. INTRODUCTION HE MOST COMMON approach for semiconductor device modeling is the drift-diffusion model (DD) given by a set of partial differential equations (PDE) (poisson’s plus current continuity and transport) [I], [2] to be solved for the electric potential (+), electron and hole concentrations ( n , p ) , respectively, or an equivalent

T

0278-0070/88/0200-0259$01.OO O 1988 IEEE

IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN, VOL. I, NO. 2, FEBRUARY 1988

260

set of variables (for instance and pseudopotentials). This model is derived from Boltzmann’s transport equation under certain approximations and its validity is, therefore, restricted to cases where the carriers are in equilibrium with the local electric field [3]; furthermore, it makes use of macroscopic parameters, such as the carrier mobility p and the diffusion coefficient D , which can be determined only experimentally. In modern VLSI devices, however, the internal electric fields are continuously increasing due to a shrinking of the geometrical dimensions unaccompanied by a suitable reduction in voltages; consequently high-field, nonlocal phenomena that cannot be described by means of the DD model progressively have become more and more important. An approach to these phenomena is to consider also the energy and momentum balance equations [4]-[6] in the relaxation time approximation. The resulting model ( D D energy and momentum equations ), often called the hydrodynamic model, although reaching beyond the simple drift-diffusion theory, is still macroscopic, besides being critically dependent on parameters such as mobility, diffusivity and relaxation times whith a dependence on field and/or carrier energy that must be determined empirically. Because of its inherent macroscopic nature the hydrodinamic model cannot describe the electron energy distribution with the detail that is essential to study any phenomenon characterized by an energy threshold, such as gate and substrate currents in silicon MOSFET’s and negative differential resistivity in GaAs MESFET’s. The most complete approach to device modeling features Boltzmann’s transport equation without the relaxation time approximation and, of course, Poisson’s equation. Unfortunately the exact solution of the resulting problem is, in general, not available. A numerical method for solving Boltzmann’s transport equations is the Monte Carlo technique (MC) that first was applied to the charge transport problem by Kurosawa [7] and later used to solve more complete and realistic cases in one [8], [9] and two dimensions [lo], [111. The MC approach to charge transport consists of following in space and time the history of each individual particle. A complete set of information, including carrier energy and velocity distributions, is obtained from which the results provided by any of the macroscopic models can easily be derived by means of suitable averages. The major drawbacks to MC approach come from the intensive CPU requirements and difficulties in handling real multidimensional devices. In this paper, a 2-D MOSFET Monte Carlo Simulator ( MOS2) that solves some of the mentioned problems will be presented. The program has been developed on a MICROVAX’) minicomputer and is therefore suitable for extensive use in technology development of VLSI devices where high field effects that cannot be described by conventional device analyzers are crucially important.

+

‘MICROVAX is a registered trademark of Digital Equipment Corporation.

The paper is organized in the following way: in Section I1 the Monte Carlo method is briefly described; in Sections 111 and IV the main features of our program will be illustrated. Section V contains selected applications and comparisons with experimental data; finally some conclusions are drawn in Section VI. 11. THE MONTECARLO(MC) METHOD The MC method applied to charge transport in semiconductors [ 121 consists of the simulation of one or more carriers subject to external forces (electric field) and given scattering mechanisms. The carrier free-flight duration (i.e., the time between two successive collisions) and the scattering event at the end of each flight are selected stochastically, according to given probabilities describing the microscopic processes. Therefore, the method relies on the generation of a sequence of random numbers with given distribution probability, usually obtained from a computer generated sequence of random numbers evenly distributed between 0 and 1. Each scattering mechanism is then described by a microscopic model that relates the electron parameters (energy, momentum, * * ) before and after the collision by appropriate scattering probabilities so that after collision a particle can start a new free flight with a well-defined new state, and the process is continuously repeated. If the purpose of the simulation is to describe a steady state, it is sufficient to follow only one carrier at a time across the device, and then repeat the procedure until the results are statistically significant. On the contrary, when dealing with non stationary processes, a large number of particles has to be simultaneously simulated. These two schemes are usually referred to as single- and multiparticle, respectively. In our case, we have so far limited our analysis to stationary device characteristics implementing the former of these two approaches-in a simulator that is described in the following section.

-

OF MOS2 111. DESCRIPTION A. MOS2 as a Post-Processor MOS2 is a 2-D single-particle Monte Carlo simulator for the transport of electrons in MOSFET’s. To save computer time the MC package has been conceived as a postprocessor of a traditional 2-D numerical simulator [13] based on the drift-diffusion model. The DD model provides the electric field used within the simulation domain of the MC simulator. Expensive iteration between Poisson’s equation and the MC transport description usually required to make the whole calculation self-consistent [141 is approximated with a simple two-step procedure illustrated in the flowchart in Fig. 1. First a complete simulation of the device is performed with the 2-D drift-diffusion package obtaining the potential ($), and the electron and hole concentrations (n, p ) at each node of the discretization mesh. Second, the space domain for the MC simulation is chosen as a sub-region of the initial device based on a criterion described in Section 111-F, the

SANGIORGI et al. : M0S’-EFFICIENT DOPING

MONTE CARLO SIMULATOR

TOPOLOGY

26 1

flight; first, unacceptable errors in the simulation of regions with very inhomogeneous electric fields are avoided; second, the trajectory computation, following the particle element by element, takes care of the assigning and identifying of the locations of the particles during the simulation with no extra effort, even using a nonuniform triangular mesh; iv) the results of the MC simulation can be averaged so that at each grid node we can compare the results given by the MC with those obtained with the DD simulator.

CONCENTRATION

CHECK FOR

2-D POTENTIAL

ELECTRON CONCENTRATION

M C RESULTS

Fig.

Flowchart of the two-step simulation proce

re.

2-D potential distribution is used to compute the electric field driving the electrons, and the MC calculation is performed. Finally, as the present scheme is not self-consistent, tests are performed to assess the agreement between the results of the MC and of the DD analysis where comparison is possible. In particular this is done for the electron concentration that is the crucial ingredient of Poisson’s equation not reconsidered after the MC calculation. The comparison of DD and MC results has been extensively used in the development of MOS’ to make sure of the validity of the MC simulator. For most of the cases of interest the previous check shows very good agreement because the potential is mainly determined by the bulk of the electron distribution, well described by the DD model. Only for very short devices an increasing disagreement between the two models appears, thus suggesting the need of a more sophisticated self-consistent scheme. To find the field acting on the electrons, a partition of the triangular mesh is performed following the ‘perpendicular bisector’ method and defining a control volume ( V , ) for each node. The potential is then supposed linear throughout the volume V, and a 2-D step-graded electric field distribution is obtained with a constant value around each node. This approach has several advantages since: i) it allows to make use of sophisticated algorithms available in the DD simulator to generate a nonplanar triangular grid with adaptive refinements in the regions of interest (junctions, high field regions, inversion layers, etc.) that are indispensable to deal with real, nonplanar devices; ii) within each V, because the electric field is constant, the electron motion equation +

-+

A k = - qE (1) can be integrated analitically ; the global electron trajectories are then made up of parabolic segments corresponding to the paths within each V,; iii) such a way to compute trajectories requires several subtle numerical problems to be solved but it presents advantages on simpler methods [ 141 where the initial electric field is kept constant during the free

B. The Self-Scattering Scheme If y ( k‘( t ) , 3 ( t l )dt is the probability that an electron with wave vector k and position 3 suffers a collision during the time dt, then the probability 6 ( t ) dt that it will dt after the previous be scattered between time t and t scattering event is given by

+

6 ( t ) dt = y ( k ’ ( t ) ,7 ( t ) )

. exp

[- ir

1

y ( k ‘ ( ~ )7, ( ~ ) d7 ) dt.

(2)

Consequently in MC simulations a stochastic free-flight time sequence with the same distribution as (2) must be generated. Starting from a random number r evenly distributed between 0 and 1, the flight time tf will be given by

~I@(T) d7

=

r

(3)

Solving (3) is difficult and too time consuming for the extremely repetitive use required in MC simulations. Consequently, the following simple method has been proposed [15], [16] to ac_hieve the same result: if r is the maximum value of y ( k ( t ) , 7 ( t )) during the whole simulation, a new artificial scattering mechanism, the so called self-scattering, is introduced to make the total scattering probability per unit time equal to F and independent of the electron wave vector and energy. When a selfscattering occurs, the electron state is assumed not to be changed by the collision, so that the electron trajectory continues as if no scattering occurred,The self-scattering probability per unit time for a st_ate ( k ( t ) , 3 ( t ) )is then automatically defined as r - y ( k ( t ) , 7 ( t )). Under these conditions (2) can be rewritten simply as

6 ( t ) dt

=

dt

(4)

By means of the direct technique [12] the flight time duration with the distribution of (4)is given by 1

tf = --In

r

(r)

(5)

and for each free flight, tf is a simple function of the selected random number r .

262

IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN, VOL. I, NO. 2, FEBRUARY 1988

C. Choice of the Scattering Mechanism During each flight the electron dynamizs is governed by (1), so that the energy E, wave vector k , and position 7 are known at any instant and the probability of each scattering mechanism can be evaluated at the end of the flight. The type of collision terminating the free flight is then chosen according to the ratio between its probability y i and r.

D. Average Quantities In order to obtain results from a Monte Carlo steadystate simulation comparable to those of the DD approach, suitable averages must be computed. In our program we used the synchronous-ensemble method [ 171, [ 121 according to which the average value of the quantity A ( 3) can be evaluated as

The carrier concentration at the nodej, instead, is given by (11) where the sum is made over all the flights ended within the jth volume. Of course, the normalization constant H must be used in calculating the terminal currents, as given by

I,

=

qHz1

( 12)

where X1is the net number of carriers collected at thejth electrode. Because the actual simulation is performed only in two dimensions (i.e., the third component of the applied electric field is equal to zero), the currents are calculated per unit of device width.

E. The Variable r Scheme The self-scattering scheme, although very simple to imk plement,_can be very expensive in terms of CPU time where C is a normalization constpt, and n ( i , 7) the when y ( k ( t ) , 3 ( t ) ) during the simulation varies over a electron+distribution function in ( k , 7). It can be shown wide range. This case arises when the difference between that n ( k, 7) is given by the maximum and minimum values of scattering rate is very large because (5) wi! produce sequences of short (7) free-flights even when y ( k ( t ) , 7 ( t ) ) is near the minimum, namely for low electron energies (for example in a where nb( k , 3 ) is the distribution function before scat- low field region) where the actual free-flight is continutering and is propoGiona1 to the probability that an elec- ously intempted to perform self-scatterings. In MOSFET's this situation is unavoidable because the tron+ is found in ( k , 7) immediately before scattering, y ( k, 7) the scattering probability of an electron in ( k , domain of interest includes the source/drain low field re7 ) , and T~ an appropriate constant. Equation (6) then can gions (indispensable for correct boundary conditions) , and a channel where a strong electric field substantially inbe written as creases the camer energy. Even worse, because of obvious concentration and dimension differences the over( A ( ; ) ) =- i whelming majority of the simulated electrons is within the y;' S/D regions, and an enormous amount of CPU time is i wasted in simulating their uninteresting self-scattering and, in the case of total constant scattering rate (self-scatevents. Quantitatively the scattering rates can be as high tering scheme), simply as as lOI5 s-' near the end of the channel and as low as 10" s-' in the low field regions, so that an electron travelling (9) within the source region suffers a true collision only every io4 self-scatterings. where the sum covers all N electron free-flights and Ab, A more efficient proposal by Hockney [18] and Mogindicates the value of A before the ith scattering event. lestue [14] is that of choosing instead of a unique value With this method, average values and distribution func- of I', a number of constant levels ( r,) each applicable in tions of energy, momentum, velocity, particle space dis- the state domain where the scattering rate does not exceed tribution, etc., are calculated and progressively refined r r .When an electron starts a flight in the ith domain, the during the simulation. In particular, the camer concentra- tf is+calculated using r, in (5); if at the end of the flight tion can be evaluated from the particle space distribution y ( k ( t f ) , 7 ( t f ) )is greater than I',, then r, is increased by means of an appropriate normalization factor taking (for example replacing it with + ) and the flight is reinto account the relatively small number of simulated par- peated, possibly many times, until the final value of scatticles. In MOS2, this normalization factor H is obtained tering rate is less than the ultimate I', . However when this by imposing the condition that both the drift-diffusion and method was implemented in MOS2, the results were inthe Monte Carlo simulators produce the same total mobile consistent both with those obtained from the drift-diffucharge within the simulation domain. Thus sion simulator and also with those from Monte Carlo schemes with constant r (in good agreement with one other), even in simple nf-n-n+ 1-D structures. In partic-t

c

SANGIORGI er

U/:

MOS'-EFFICIENT

263

MONTE CARLO SIMULATOR

ular, the current calculated in any section of the device as the flux of the mobile charge was not solenoidal and did not coincide with that given by (12). In the spirit of the synchronous-ensemble method [ 171, it can be noticed [ 121 that, the weighting factor used in (8) must be a unique function of k' at the end of the free flight whereas, _when using the multiple r scheme, the same final state ( k ( t ) , ? ( r ) ) can be weighted differently depending on the electron previous history. This inconsistency has been overcome in the following way. As already mentioned, in order to reduce the number of self-scatterings while using a simple relationship between the random number and the free-flight duration time, it is very useful to deal with a step-shaped probability distribution. Such a distribution can be obtained by splitting th_e domain of the particle wave vector and positio? D( k ( t ) , 7 ( t ) ) into a number of sub-regions D,(k ( t ) , 7 ( t ) ) and defining a self-scattering rate yII which makes the total scattering rate r, constant within each of them. The particle trajectory, made up of portions comjng from each V,, can be evaluated in the domain D( k ( t ) , 7 ( t ) )obtaining the sequence of the sub-regions (01, . . . , D,} crossed by the electron and the corresponding sequence of maximum scattering rates (includwhere we allow ocing self-scattering) { rl, . . . , r,?} currences such as D,= D, (r, = r,) indicating that the electron has come back to a previous state. For a given

and the time of free flight then is given by

In MOS2 each D,(z ( t ) , 7 ( t ) ) is made by the-? in the same volume around a mesh node and by the k corresponding to a specified energy range { E,, E, + } ( E is a function of I k I only). Notice that the particle trajectory is independent on the scattering rate that only limits the length of the electron path along it; thus the sequence { t , , . . . , t, } can be evaluated without any assumption about the scattering rate. Equation (17) can, therefore, be implemented in the following manner: for each free-flight, first a random number r between 0 and 1 is chosen and rlog = -In ( r ) is evaluated ( rloeis always 2 0 ) ; then the time interval to reach the next sub-region ( A t ) is computed and the corresponding scattering rate r extracted from a table. This procedure is repeated, adding the products A t F l , and is terminated just before C A t , r , = rprev (Fig. 2) exceeds rlog.Finally the time spent in the last subregion is computed, dividing the difference between Y, t, if j > i ), the real scattering rate (without self-scattering) as a functhe probability 6 that a particle scattered at the time to = tion of time, evaluated only for t = t f . The dotted, step0 will scatter again between t and t + dt is given by (see graded line instead represents the values of the discretized (2)): total scattering rate I?, (including self-scattering) in each 1-1 sub-region crossed by the particle; finally the solid line s ( t )dt = r, exp [ - ( t - rl)rl(t, - t,-l)r,]dr shows how the flight time rf is computed for a given -In J=I ( r ) summing the factors At, r,. As for the average quantities this algorithm is obviously (13) consistent with (8) where yI must be replaced by the total where t, - < t < t , , and rJ is the time of exit from thejth scattering rate r, associated with the subregion where the domain DJ.Thus (3) allows us to compute the relationship ith flight terminates. between random numbers and free-flight times. In fact, deriving (13) F. The Choice of Boundary Conditions In order to save computer time the MC simulation domain should be restricted to the relevant channel region. and, because 6 (O)/r,= 1 and 6 ( t ) / r ,always is con- With this procedure, however, the electric field at the boundaries is significant and the choice of boundary continuous, ditions (BC) (i.e., the parameters of the entering elec1 trons) is not trivial. I? practice one should know the elec6 ( 7 )d7 = 1 - - 6 ( t ) . (15) tron distribution in k that is a priori unknown, and we rI have found that the use of a special (guessed) distribution Then using (3), (13), and (14) we have gives unsatisfactory results. In particular with the MC re1 - 1 gional approach [ 191 carrier concentrations and currents (t, - t,-l)I', = 1 - r (16) were found to be inconsistent with those of the DD simexp - ( t - t , ) r ,J=1 ulator, and with experimental results determined on simand, since r is evenly distributed between 0 and 1 , we can ple test structures. A solution to this problem is to place boundaries only substitute ( 1 - r ) with r on the right hand side of (16)

c

i:

[

1

264

IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN, VOL. I, NO. 2, FEBRUARY I RANDOM NUMBER O