Monte Carlo simulation for ultra-small MOS devices

Superlattices and Microstructures, Vol. 27, No. 2/3,2000 doi:10.1006/spmi.1999.0802 Available online at http://www.idealibrary.com on

Monte Carlo simulation for ultra-small MOS devices U. R AVAIOLI† , B. W INSTEAD , C. W ORDELMAN , A. K EPKEP Beckman Institute, University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A. (Received 9 November 1999) This paper discusses advanced needs of Monte Carlo simulation approaches for MOS silicon devices scaled below 0.1 µm channel length. For predictive simulation over a wide range of biases, it is necessary to provide the Monte Carlo procedure with tuning capabilities to adjust the mobility through calibration of the interface roughness scattering. This is accomplished by introducing a semi-empirical procedure with a physical elastic scattering rate and an inelastic rate with tunable strength. To resolve the role of hot carriers in relation to oxide interface damage, it is also important to realize fully bipolar MOS simulation, so that one can analyze the transport of impact-ionization generated carriers and secondary ionization. As the devices become quite small, three-dimensional simulation can be not only feasible, but also necessary to resolve the granularity of doping profiles and the complete carrier–carrier and carrier–ion interactions. Issues of device Monte Carlo implementation on parallel environments are discussed, and a practical approach for resolving the short-range forces of the charge–charge interaction in three dimensionals is described. Several examples and preliminary results are presented to illustrate the various issues. c 2000 Academic Press

Key words: Monte Carlo simulation, MOSFET, hot carriers.

1. Introduction The semiconductor industry is fast approaching the limit of 0.1 µm feature sizes for commercial integrated MOSFETs. It is reasonable to expect technology to evolve and reach perhaps 0.05 µm gate lengths with somewhat conventional MOS devices, but it is clearly perceived that between 0.05 and 0.01 µm (100 Å) the scaling of MOSFETs, as we know them today, will have necessarily to come to an abrupt stop due to intrinsic physical limitations. The drift-diffusion approach with tunable mobility models is still the workhorse for routine TCAD applications, well into the sub-micron device scale. However, despite the scaling of bias voltages to reduce high fields and hot carriers, high-end physical simulations such as particle Monte Carlo and other solutions of the Boltzmann transport equation (BTE) have become more and more necessary, to understand nonequilibrium and hot carrier phenomena and from time to time to calibrate models at lower levels in the hierarchy. With deep scaling we are also entering a range where granularity of the device will become a marked feature, since, at doping concentrations between 1018 and 1020 cm−3 , the ions will have average space separations approximately from 100 to 20 Å. In silicon devices we also need to make the additional consideration that, at room temperature, the mean-free path of carriers should stay below 100 Å, therefore † Author to whom correspondence should be addressed: Umberto Ravaioli, 3255 Beckman Institute, University of Illinois, 405 N. Mathews Avenue, Urbana, IL 61801, U.S.A. E-mail: [email protected]

0749–6036/00/020137 + 09 $35.00/0

c 2000 Academic Press

138

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

we expect size quantization effects and tunneling to be much more important than quantum coherence effects along the transport path. In the last stretch of MOSFET scaling, from 0.1 to 0.01 µm, transport physics at room temperature should be in a transition area, between semi-classical and quantum regimes, characterized by system granularity. Transport should have at the same time distinguishable semi-classical and quantum features, in a proportion that will vary as scaling is pushed further to smaller dimensions. An established simulation framework that can effectively deal with granularity issues is lacking, mainly because a full 3D simulation is required to include correctly carrier–ion and carrier–carrier interactions. In the practical range of devices that will be built in the foreseeable future, a particle Monte Carlo technique is the natural choice to resolve granularity and to ‘bridge’ between the semi-classical regime and the purely quantum regime of operation. While the usual semi-classical Monte Carlo procedure can be used to resolve motion and scattering, quantum features such as size quantization and tunneling can be prescribed by altering attributes of the simulated particles, in specific regions and in a self-consistent manner. In fact, Monte Carlo simulation techniques can go well beyond the typical range of applicability of the semi-classical BTE, and although a Monte Carlo simulation is computationally expensive, there is little penalty in the transition from 2D to 3D simulation. In this paper we will discuss a number of issues related to physical models and Monte Carlo implementation for practical simulation of ultra-small MOS structures.

2. Interface roughness scattering Many tests have shown that the semi-classical transport model implemented in Monte Carlo simulators including the full band structure are reasonably accurate for silicon [1]. In MOS devices, however, interface roughness scattering at the Si/SiO2 boundary dominates the mobility. The interface morphology has an important role, and for any specific industrial process it would be necessary to know the resulting statistics of interface fluctuations amplitude and the correlation of the fluctuations, which enter physical scattering rate formulations. Few attempts have been made to determine precisely the statistical nature of the interface [2] and to implement physical scattering rates with detailed transverse wavefunctions in the channel [3]. The standard practice in semi-classical Monte Carlo simulation is to account for interface scattering processes by treating the particles reaching the oxide boundary either with a perfectly specular reflection or with a diffusive reflection, chosen in a random fashion. The statistical ratio between specular and diffusive processes can be adjusted to match a measured mobility. However, while performing tests with measured data, it is difficult to match closely current–voltage curves over the complete bias range with a single value for this ratio. In order to obtain a calibration procedure that avoids the need for demanding measurements of interface morphology and exceedingly complicated models that imply sharply increasing computational cost, we tried to implement a tunable semi-empirical model that mimics the interaction as much as possible. We applied two scattering mechanisms to particles on a very narrow layer under the interface. The first rate corresponds to an elastic mechanism, implemented like the physical mechanism discussed in [3] which best describes the behavior of low-energy carriers in low field conditions. This rate requires the knowledge of the correlation length and the amplitude of the roughness at the interface. The correlation length can be selected following typical values used in the literature, while the r.m.s. height of the roughness is taken as a free fitting parameter to calibrate the elastic scattering rate. At higher fields, it is necessary to provide a dissipation mechanism, which we implement and calibrate empirically in the absence of a known detailed model. The process is treated similarly to an optical phonon, which is inactive at low energies and gradually kicks in at higher energy values, by the use of a smooth function. The shape and strength of this function can be varied to fine tune the overall transport behavior and match a measured mobility. Monte Carlo simulation examples are shown in Fig. 1 for a wide MOSFET with conductive gate as employed in the classic measurement reported in [4]. The good agreement shows that it is relatively straightforward to match a

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

139

1000 900 800 700

Mobility (cm2 V–1 s–1)

600 500 400 300

200 Na = 3.9 × 1015 [4] Na = 7.2 × 1016 [4] Na = 3.0 × 1017 [4] Na = 7.7 × 1017 [4] MOCA (Na = 4 × 1015)

100

0.05

0.1

0.2

0.5

1

Effective field (MV cm–1) Fig. 1. Comparison of Monte Carlo simulations (×) calibrated with the measurement of mobility in the inversion layer [4].

wide range of conditions with the proposed tunable model. Work in progress is now devoted to perfect the calibration procedure for regular MOSFETs. A possible calibration procedure could start by selecting a reasonable guess for the r.m.s. roughness height parameter, to get a first acceptable fit for low field channel resistance. After that, the strength of the high-energy inelastic rate can be tuned to fit the saturation current and then a more precise fine-tuning of the r.m.s. roughness height parameter can be completed. For radically different interface processing or oxide growth conditions, different values of the correlation length of the roughness may have to be re-evaluated from measurements, first principle evaluations or other suitable means. One should keep in mind that when actual measurements are fitted as described above within a semiclassical Monte Carlo implementation, the interface quantum effects are also implicitly accounted for by the semi-empirical procedure. The main effects would be related to the discreteness of energy subband levels as a consequence of size quantization and the inter-subband scattering phenomena, which are not treated by a semi-classical transport approach. If the quantum effects were included separately in the overall simulation, the calibration procedure would only cover interface interaction effects, leading therefore to a different set of parameters. 2.1. Bipolar-mode simulation for reliability studies The unique feature of Monte Carlo simulation is to resolve the individual particle histories. While this implies a considerable computational cost, only a technique like Monte Carlo can be used to test the physics and statistics of physical mechanisms involving events, which are rare on the typical time scale of transport simulation. Device reliability is one such key area of study. The traditional empirical models link oxide in-

140

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

terface damage to substrate current measured in MOSFETs under stress-bias, and simple physical models of damage involve injection of both holes and electrons across the oxide interface. However, as device sizes and biases are continuously scaled down, conditions are now being reached where little or no impact ionization would take place in normal operation. At this point one has to question the fundamental basis of the traditional approaches and pay more attention to detailed physical model of the carrier–interface interaction. The standard way to passivate dangling bonds at the surface is to place hydrogen at the interface. When hydrogen atoms are progressively kicked out of their sites by interaction processes, charge builds up at the interface with eventual degradation of device threshold and transconductance. The discovery of the hydrogen-deuterium isotope effect at the oxide interface [5] has provided an important reference to examine the physics of degradation. Atomic-level calculations are beginning to indicate that a much richer spectrum of bound states is possible for hydrogen atoms at the oxide interfaces [6]. In addition, one should expect it possible to desorb hydrogen atoms from the dangling sites they occupy, by a series of lower energy electron impacts, since hydrogen requires a considerable time to relax energy to the lattice from excited bound states. Deuterium is instead hard to displace because it very efficiently relaxes energy to the lattice phonons. It is important to link atomic physical models with a detailed transport formalism to complete the connection between theory and experiments, in particular to prove the multipleimpact model of hydrogen desorption. The transport model should provide not only distribution functions in the channel, but also a precise statistics of particle impacts on the interface, including particle energy, momentum and angle of incidence at the time of impact. This is possible only with Monte Carlo simulation. We have therefore devoted considerable effort to developing a combined electron–hole simulation of MOSFETs where the physical basis of damaged could be evaluated consistently under various conditions of stress and regular operation bias. The simulator solves for individual holes in the substrate as well as for impact-ionization generated holes, which can be tagged and followed during the simulation, to understand their history and possible role in damage events. Two similar programs for full band electron and hole simulations are run synchronously under a single- or multiple-processor environment, feeding at each time step the corresponding space charge distribution to a Poisson solver that updates the electric field and potential. With this software architecture, the overall simulation can be easily managed and parts of the code can be reused for other applications. All the interface impact events are recorded to build up a statistical database that can then be used to test theories of damage and to attempt a quantitative correlation with available measurements. An example of collected data is shown in Figs 2 and 3, which were obtained for a 0.05 µm n-channel device with complete structure details posted on the Web page of the ‘Well-Tempered-MOSFET’ Project [7]. In these pictures we show the distribution of electron and hole impacts onto the oxide interface, collected into bands of energy for clarity. The results are obtained from long simulation runs, which are then scaled to provide an estimate of impacts for each mesh along the interface, per micron of width in the third direction and per year. This way of presenting data is practical to attempt the correlation of microscopic simulations with damage phenomena that may involve hours to days in stress-bias conditions, or months to years at normal operation biases. The same procedure can be iterated to actually map the possible evolution of the damage over long times, by sequencing quasi-steady-state Monte Carlo simulations for the same device. Surface charge can be progressively adjusted to map the damaged predicted by the interface model, based on the transport Monte Carlo data obtained at the previous step. The computational cost of these procedures is in general considerable, but the cost decreases as devices are scaled, and with the progressive reduction of biases a detailed simulation may simply become the only possible analysis option. For the 0.025 µm device posted in [7], our Monte Carlo simulations do not record a single-impact ionization event at the suggested operation biases, which are low enough to keep the vast majority of electrons well below the ionization threshold. Even if one could manage to operate this device at a higher stress bias for impact ionization to occur, any extrapolation of damage trends from a gate current measurement would become at best questionable. Considerable work

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

141

0 0.05 Xecross;Elec; n – wtm 50 nm; Vd = 2.0 V; Vg = 1.0 V

0.1

Number of grids /µm /y

1025 0.0 0.5 1.0 1.5 2.0 2.5 1020

1015 0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

x (µm) Fig. 2. Statistical data collected from Monte Carlo simulations for electron impacts on the oxide interface as a function of energy in a 0.05 µm channel MOS device as specified in [7] with drain bias Vd = 2.0 V and gate bias Vg = 1.0 V. The legend associates with each symbol the lower bound of the represented energy band, in units of electron volts.

is now in progress to prototype the coupling between Monte Carlo simulation and atomic-level interface simulation. 2.2. 3D Monte Carlo simulation Device scaling and the corresponding reduction in the number of particles in a structure make a full 3D Monte Carlo simulation approach much more affordable. The extension is not only desirable for the sake of completeness, but also to resolve charge–charge interaction, which is not yet fully understood to the extent that it affects the details of the high-energy distribution function. This may have considerable bearing on quantitative assessment of damage models. A 2D simulation cannot naturally resolve the complete shortrange carrier–carrier and carrier–ion Coulombian interactions, and the only available approach is to add to the solution of Poisson equation (providing the long-range interaction) equivalent scattering models that inevitably have to become questionable when devices are extremely scaled. While the shortcomings of all standard ionized-impurity scattering approaches have been known for a long time, since failures may be observed at high doping levels or low temperatures, still one could compensate for those problems with empirical adjustments from mobility measurements. However, in ultra-small devices the notion of smooth doping profiles departs severely from reality, since doping becomes granular and individual ions should be resolved. Because of statistical fluctuations, devices with nominally identical doping profiles might exhibit considerable differences in their actual behavior due to differences in ions placement, even if exactly identical in total number. As an example, device threshold is one of the important parameters that would dramatically suffer from granular doping profile fluctuation [8]. When considering carrier–carrier scattering rates, it is not

142

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

0 0.05 Xecross;Hole; n – wtm 50 nm; Vd = 2.0 V; Vg = 1.0 V

0.1

1020

Number of grids /µm /y

1018

0.0 0.5 1.0 1.5 2.0

1016

1014

1012

1010 0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

x (µm) Fig. 3. Statistical data collected from Monte Carlo simulations for hole impacts on the oxide interface as a function of energy in a 0.05 µm channel MOS device as specified in [7] with drain bias Vd = 2.0 V and gate bias Vg = 1.0 V. The legend associates with each symbol the lower bound of the represented energy band, in units of electron volts.

clear how applicable the conditions are under which they are derived or the screening models on which they depend, once the extreme granular regime is approached. In order to eventually make 3D simulation more practical, some of our efforts are focused on developing a modular simulation framework adaptable to parallel super-clusters based on typical low-cost PC-processors. At present, target systems are the large clusters available within the National Computational Science Alliance (the new phase of NCSA) with platforms based on NT and Linux operating systems. In order to effectively map MOSFET devices to the processors, we partition the simulation domain by slicing with parallel planes running from source to drain, so that each slice would be fairly naturally balanced with an approximately equal number of particles. By treating each device slice with a processor, the main communication work consists of handling plane crossings between neighboring processors, with no ambiguous situation of trajectories in the vicinity of corners which would arise from additional partition of the slices themselves. If both holes and electrons are simulated, two processors could be assigned to each slice, depending also on the availability of memory at each node to store the necessary band structure and scattering tables. The overall parallelization flow is therefore fairly simple, with the various processors updating at each iteration the particle positions, performing scattering events and exchanging particles crossing over slice boundaries. At the end of each time step the fields are updated by solving Poisson equation, for instance, or a more complete coulombian interaction solver including short-range effects. We have been working on approaches for the complete inclusion of carrier–carrier and carrier–ion shortrange interactions, although all the computational applications have so far been demonstrated only on single processors. Our goal has been to implement an efficient method, which preserves the flexibility of nonuniform rectangular meshes used for the long-range Poisson solution, necessary for practical applications, and that

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

× 107

143

Electron velocity versus channel position

2.5

2D 3D 3D ci 3D cicc

vx (cm s–1)

2

1.5

1

0.5

0

–0.5 30

40

50

60

70

80

90

100

110

120

x (nm) Fig. 4. Monte Carlo results for average velocities along the channel of a simple 0.05 µm MOS structure using: a 2D approach, a standard 3D approach with smooth doping profiles and ionized-impurity scattering rates, a 3D approach including doping granularity with explicit carrier–ion interactions (ci), and a 3D approach with both carrier–ion and carrier–carrier interactions (cicc).

minimizes the computational cost. We have extended to device simulation a technique known as particle– particle–particle–mesh (P3 M) technique [9], where forces are obtained from a local molecular dynamics calculation performed on a short-range domain around each particle (particle–particle force), coupled with the overlapping mesh where the Poisson equation is solved (particle–mesh force). It is very useful to retain the Poisson equation for the long-range force calculation, because boundaries and field line distortions at dielectric interfaces are naturally included while these would be quite problematic issues for a molecular dynamics calculation extending to the entire simulation domain. Since we choose the radius of the domain for short-range calculation to span at least two meshes in all directions [10], the two evaluation methods provide force components that when added contain a certain degree of double counting, as a result of the overlap of the two solution domains. The short-range Coulombian force between two particles i and j is evaluated as FiCoul = j

qi q j (ri − r j ) 4πε|ri − r j |3

but a correction must be applied for elimination of double counting, to obtain a short-range force Fisrj to be added to the Poisson equation result, where Fisrj = FiCoul − Ri j . The excess force Ri j is defined reference force j in [9]. Although an exact evaluation of the reference force is not possible in the general case, one can obtain a reasonable approximation by assuming that each particle under consideration is associated with a space charge distribution, peaked at the center of the short-range sphere and going to zero at the sphere boundary. For any appropriate choice of this charge distribution in the short-range domain one can obtain an explicit reference force result to correct the Coulombian force evaluated by the molecular dynamics procedure. The

144

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

effect of the correction is to provide a smooth transition, from the short-range domain to the long-range domain where the Poisson equation solution dominates. This method has been tested by evaluating the low field mobility of silicon for a wide range of dopiness in a 3D resistor-like structure where the actual doping ions were placed at randomly chosen locations while no ionized-impurity scattering rate was invoked [10]. Similar results were shown [11] with a technique where the forces obtained from the Poisson equation are corrected using pre-tabulated error tables obtained from molecular dynamics calculations, with the only apparent drawback of being limited so far to uniform mesh. A prototype MOSFET simulator has been implemented to fully extend the P3 M technique to device applications. In Fig. 3 we show an example of recent preliminary results for average electron velocity in a simple structure with 0.05 µm channel length and 1 V applied to gate and drain. We compare several simulation approaches, including a 2D model and a 3D model with only the standard Poisson solution and ionized-impurity scattering rate, a 3D model with granular carrier–ion interaction, and a 3D model with both carrier–ion and carrier–carrier interactions. There are some significant differences in the position of the peak velocity and the velocity relaxation moving towards the drain (the channel is located between 50 and 100 nm in the plot). Such results are now being examined to gain more insight into the detailed role of doping granularity and electron–electron interaction with device scaling.

3. Conclusions This paper has presented an overview of activities in Monte Carlo device simulation, addressing the specific needs of ultra-small MOS devices. While the bulk transport physics is correct in Monte Carlo, in situations where a close match with measured device I –V curves is desired, calibration capabilities for the mobility must be added since first principle models of interface interaction are too costly and not completely understood. We propose a semi-empirical procedure that mimics the complete transport distortion at the interface with two scattering rates. Interface scattering and quantum effect at the interface are intimately linked and they are both taken into account by the calibration procedure. For ultra-scaled devices one expects that traditional approaches of reliability predictions based on substrate current measurements will eventually fail, due to the lack of impact ionization at scaled bias while it is probable that lower energy degradation phenomena indicated by atomistic calculations will gain importance. Bipolar Monte Carlo simulation can be used as a tool to understand this transition and to provide a direct link between atomistic interface calculation and experimental observations. Finally, 3D simulation is becoming a necessity for ultra-small devices, because of the importance of resolving doping granularity and carrier–carrier interaction issues. Super-clusters of PC-based processors promise to be an economical way to exploit parallelization schemes for 3D simulation. We have shown a practical approach to account for the complete Coulombian interaction in a Monte Carlo simulator, without appreciably altering the existing architecture of the codes, since local molecular dynamics calculation for short-range forces can be coupled to the existing Poisson equation framework. The last area that needs to be addressed in the near future is the inclusion of quantum effects in the Monte Carlo procedure. Deep scaling is expected to bring new effects that have rarely been considered in devices, such as current enhancement due to tunneling through the source/drain barrier. We believe that the most efficient way is to adopt an approach in which the potential is locally modified to account for quantum effects, while preserving the semi-classical character of the particle trajectories. This can be accomplished by basing the transport on quantum-corrected versions of Boltzmann equation. Such an approach avoids the complication of coupling semi-classical transport to a local quantum system, which still presents considerable numerical and theoretical challenges. Work is in progress to include a quantum correction scheme to the full band Monte Carlo framework discussed in this paper. Acknowledgements—This work was partly supported by the National Science Foundation, through the Distributed Center for Advanced Electronics Simulation, grant ECS-9802730, the Semiconductor Research Corporation, contract 98-SJ-406, and an equipment grant by the IBM Shared University Research Progam.

Superlattices and Microstructures, Vol. 27, No. 2/3,2000

145

References [1] K. Hess (ed.), Monte Carlo Device Simulation: Full Band and Beyond (Kluwer Academic Publishers, Norwell, MA, 1991). [2] S. M. Goodnick, D. K. Ferry, C. W. Wilmsen, Z. Liliental, D. Fathy, and O. L. Krivanek, Phys. Rev. B32, 8171 (1985). [3] S. Yamakawa, H. Ueno, K. Taniguchi, C. Hamaguchi, K. Miyatsuji, K. Masaki, and U. Ravaioli, J. Appl. Phys. 79, 911 (1996). [4] S. Takagi and A. Toriumi, IEDM Technical Digest p. 711 (1992). [5] J. W. Lyding, K. Hess, and I. C. Kizilyalli, Appl. Phys. Lett. 68, 2526 (1996). [6] B. Tuttle, Phys. Rev. B60, 2631 (1999). [7] World Wide Web page address: http://www-mtl.mit.edu:80/Well/. [8] A. Asenov, Electron Device Lett. 45, 2505 (1998). [9] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (McGraw-Hill, New York, 1981). [10] C. Wordelman and U. Ravaioli, Integration of a particle-particle-particle-mesh algorithm with the ensemble Monte Carlo method for the simulation of ultra-small semiconductor devices, to appear in IEEE Trans. Electron Devices. [11] W. J. Gross, D. Vasileska, and D. K. Ferry, IEEE Electron Device Lett. 20, 463 (1999).