Recent Advances in Device Simulation Using Standard Transport Models

In this paper we address a number of issues related with device simulation in 3-D, and point out a few deficiencies which still prevent 3-D device simulation to be widely accepted as a standard tool for device design and optimization in an engineering environment. More specifically, such deficiencies have to do with structure definition and mesh generation, as well as with the computational burden which is typically associated with discretization meshes featuring 50.000 nodes or more. Next, we address the problem of validating advanced physical models in 3-D by means of our device simulator HFIELDS-3D and use a flash-EEPROM cell manufactured at ST-Microelectronics as a test vehicle for the validation of hot-carrier injection into the floating gate and Fowler-Nordheim tunneling across the gate oxide.


INTRODUCTION
Numerical simulation of silicon devices has reached a mature development stage.Two-dimensional simula- tion codes are increasingly being used both for a bet- ter understanding of device performance and for design purposes; furthermore, commercial tools fea- turing a friendly user interface and sophisticated sim- ulation environments are currently made available to interested users.Three-dimensional simulation, instead, is still in an evolutionary phase, since a number of unsolved practical and theoretical problems still exist, which make it difficult to achieve an accurate prediction of the device performance.Such problems range from structure definition to mesh generation also, the size of the computational problem is often too large to be properly handled in an engineer- ing design environment.In addition, some physical models turn out to be very hard to quantitatively predict.
In this paper, we illustrate results achieved in the simulation of a flash-EEPROM cell designed for a 4  Megabit memory chip and, in doing so, we point out problems and inefficiencies arising from the incom- plete availability of the necessary software tools.More specifically, we discuss the problem of structure definition and mesh generation and, from the simula- tion of the programming and erasing characteristics of the cell, address the problem of validating the physi- cal models which are being used for the above func- tions, i.e. hot-electron injection into the floating gate, and Fowler-Nordheim tunneling across the gate oxide, respectively.

STRUCTURE DEFINITION AND MESH GENERATION IN 3-D
A major problem with structure definition in 3-D is the lack of reliable 3-D process-modeling tools.Generally speaking, such tools have a lesser predicting capability than device simulators, and this is especially true in 3-D, where process modeling is still in an early development stage.As a result, the device morpholgy is currently inferred from TEM cross sec- tions, which involves a rather cumbersome procedure.
Also, the impurity concentration is inferred from 1-D and 2-D process simulations empirically com- bined to get the 3-D profiles.However, due to the dif- ficulty of measuring experimental profiles in 3-D, the only feedback comes from the comparison between simulated and measured device characteristics.Quite often, a reverse-engineering process is needed, and the impurity profiles are adjusted in order to fit the turn-on characteristics.
Under such circumstances, the validation of sophisticated physical models, such as hot-electron injection into the gate oxide, turns out to be rather difficult, as the uncertainty in the device morphology and impu- rity concentration may lead to inaccurate electric fields.
Another area of importance is mesh generation.The spatial discretization of the semiconductor device equations has a substantial influence not only on the accuracy of the numerical results, but also on the con- vergence behavior.An appropriate spatial discretiza- tion must fulfill a number of requirements [1], including adaption to irregular geometries, accurate approximation of physical quantities, efficient alloca- tion of grid points and suitable element form-factors.
From the experience gained in the past, it has now become clear that no unique algorithm or discretiza- tion element can fulfill all the requirements which are needed in order to generate meshes suitable for differ- ent tasks, such as process and device modeling, simu- lation of micro-electro-mechanical sensors and actuators and so on.More recently, a new approach based on a a modular implementation of different algorithms has been devised [2].The latter allows new basic elements and properties to be incorporated within the code, without loss of generality.These modules can be selected depending on the grid requirements for the specific application.The poten- tial of the above approach raises the hope for a new generation of better and more reliable mesh generators.

PHYSICAL MODELS
Within the context of standard simulation tools, some physical effects turn out to be very difficult to model this is especially true for physical mechanisms char- acterized by an energy threshold, such as impact ioni- zation and hot-carrier injection into the oxide.In order to accurately predict such effects, the knowl- edge of the energy-distribution function of both types of carriers is a necessary prerequisite.This involves the solution of the Boltzmann Transport Equation (BTE) which is typically done using the Monte Carlo technique.If such information is not available, a rela- tionship must be established between impact ioniza- tion and electric field or carrier temperature, if the latter is available.
Within homogeneous materials, and under a uni- form electric field, there is a one-to-one correspondence between electric field, carrier temperature and energy distribution function.In deep submicron devices, however, the device length is comparable with the energy-relaxation distance: hence the energy distribution is never in equlibrium with the local field.
Hence, any impact-ionization model based on the local value of the electric field turns out to be grossly inadequate and typically leads to an overestimation of the impact-ionization effect.Likewise, hot-carrier injection into the oxide requires a non-local model to account for the dynamics of high-energy electrons.
Better results can be achieved using the hydrody- namic transport model.Within this model, a simplified energy distribution function is associated with the carrier temperature.The occurrence of the investi- gated effect is then inferred from the population of the tail of the distribution function.Such an approach is obviously not general enough to warrant its applica- bility to a variety of devices under different operating conditions.However, for a given class of devices, we expect the above assumption to be reasonably correct.

COMPARISON WITH EXPERIMENTS
In this section we are going to show simulation results of a flash-EEPROM cell which is a staked-gate device designed for the SGS-Thomson's 4 Megabit flash chip.This cell is intrinsically a three-dimen- sional device due to the morphology of the field-oxide region, shown in the TEM Fig. 1 micro- graph of figure 1, which represents a cross section normal to the current flow.For this reason, a FIGURE TEM micrograph of the EEPROM cross section two-dimensional simulation of this device cannot predict the cell behavior, and reasonable results can be achieved only after a tuning of the effective device width and after an experimental determination of the coupling coefficient between the control gate and the floating gate.The latter, however, is a function of the specific operating conditions.Hence, an accurate pre- diction of the device performance in 2-D is very diffi- cult to achieve.
Using a 3-D simulator, all capacitive couplings between the floating gate and the accessible contacts are automatically taken into account, and the scaling of the device current with respect to channel width can be investigated.Moreover, the space distribution of the generation rate due to impact ionization and band-to-band tunneling around the drain and source diffusions, in programming or erasing conditions, can be carefully studied.Therefore, the simulations are expected to predict the actual device behavior.The device cross-section parallel to the current Fig. 2 flow is sketched in figure 2.
To investigate such effects, cells with different widths were considered.It is worth noting that, in order to minimize process modifications which would introduce additional uncertainties in the geometry, the fabrication of cells with different channel widths has been attained by modifying only one process mask, FIGURE 2 Device cross-section along the channel length namely the one controlling the size of the active area.As a consequence, the dimension of the floating gate was kept constant, thus yielding a variable wing dimension.This means that the floating and con- trol-gate areas are almost the same in cells having dif- ferent widths, resulting in an unchanged interpoly capacitance, i.e. the capacitance between the control gate and the floating gate.In order to optimize the trade-off between accuracy and computational work- load, significant effort has been devoted to the generation of a suitable structure for the simulation.Rather than attempting to reproduce the shape of the bird's beak by a process simulator, the morphology of the cell was extracted from the TEM micrograph in figure to act as a template.A commercial solid modeller and mesh generator were used to generate the two-dimensional grid, which was then replicated in the third dimension.As usual, symmetry conditions were taken advantage of by considering only half of the EEPROM cell.
The fine-tuning of the cell geometry and doping profiles for the source, drain, and channel regions was pursued by reverse-engineering of the transfer and output characteristics of the cell for three different channel widths.So doing, a validation of the 3-D structure thus generated was achieved.The agreement between simulations and measurements shown in figures 3 and 4 is acceptable, in spite of a slight devia- tion in threshold voltage and saturation behavior.It is worth mentioning that no physical parameter was adjusted to fit the data.The former deviation can be ascribed to inaccuracies in the surface-channel doping and to an insufficient mesh refinement in the same region, while the deviations in the output current could be due to a limitation of the surface-mobility model in the high doping region.In any case, local changes of the threshold voltage on wafer are of the same order of magnitude.
In order to reproduce the writing and erasing char- acteristics of the cells, i.e. the curves representing the threshold-voltage shift as a function of time, it is nec- essary to determine the relation between the thresh- old-voltage shift AV and the floating-gate charge Qfg as evaluated by HFIELDS-3D.A typical assumption is that the stored charge shifts the threshold voltage according to the relation [3,4] AVt--Qfg (1) where Cpp is the interpoly capacitance.For an unwrit- ten cell it is Qfg 0, hence AV 0. During a writing operation, the building-up of negative charge within   the floating gate progressively increases the thresh- old-voltage.During an erasing operation, instead, the threshold-voltage progressively decreases due to the extraction of electrons from the floating gate.
By simulating the transfer characteristics with differ- ent amounts of charge stored in the floating gate, the linear relation between Qfg and AV expressed by Eq.
(1) is confirmed and the interpoly capacitance relative to half cell can be estimated as 1.8 x 10 -15 F for all channel widths, as shown in figure 5.
A non-linear scaling of the drain current with channel width W is experimentally observed the narrower cell carries a higher current density, as shown in fig-  reported against the control-gate voltage.This can be explained as follows.Considering for instance a transfer characteristic at V s O, Vb= O, Qfg 0, and using a simplified capacitive equivalent circuit for the cell, the floating-gate coupling to other contacts can be represented by the following linear equation [4,5]   gfg ggcg n t-(Zdgd (2) where (Zg Cp/CT is the coupling ratio between con- trol and floating gates, o d CadCT is the floating gate to drain coupling ratio, and C T is the total capacitance (C T Cpp + C d + C + CI).Eq. ( 2) simplifies to Vfg OgVcg Cpp Vcg (3) CT due to the fact that cx d is typically one order of magni- tude smaller than Og and Vd 0.1 V. Simulations show that increasing the channel width at constant gcg decreases the potential of the floating gate, leading to a higher threshold voltage.As expected, this means that, while Cpp is unchanged, C T increases, due to the increase in the planar region of the oxide.The increase in the floating-gate potential at decreasing channel widths strongly affects the erasing behavior of the cell, as shown in section 6.

PROGRAMMING OPERATION
The cell under investigation is programmed by hot-electron injection at the drain end of the channel.
Hot-eletron injection requires a large lateral electric field in silicon for the carriers to attain a high enough energy to overcome the barrier.Once they are injected into the floating gate, they are confined by the poten- tial barrier of approximately 3 eV at the polysili- con-oxide interface.
Figure 7 shows the threshold voltage shift vs pro- gramming time for different values of the drain volt- age.As more and more negative charge is injected into the floating gate, the voltage of the latter decreases until it reaches the drain voltage.At this point, the field in the oxide reverses, and slows down further charge injection.This effect is revealed by the sudden slope change visible in the curves of figure 7.
The information contained in the above measure- ments is converted to a curve relating the injection current to the floating-gate voltage, using the follow- ing expression: and FIGURE 7 Experimental threshold-voltage shift against programming time for different values of the drain voltage.Vcg 11 V 'g (Zg(Vcg AVt) --}-(ZdV d A typical example of the gate current vs floating-gate voltage is shown in figure 8. Clearly, the corrent drops off as soon as the floating gate voltage approaches the drain voltage, due to the reverse oxide field which opposes electron injection. The hot-electron injection model incorporated in HFIELDS-3D accounts to some extent for both the non-Maxwellian form of the electron-energy distribu- tion and for the non-local nature of carrier heating.Fundamentally, it is based on the expression for the  energy distribution of hot electrons f(E) derived in [6]   under some simplifying assumptions for electron transport (i.e., carriers belong to a single nonparabolic conduction band and their dominant energy-loss mechanism is the emission of optical phonons) I(E) Cexp IFI1.5 (6) where E is the electron energy, F represents the elec- tric field, C is a weak function of F, and ) =const.In order to account for non-local effects, in [7] the actual electric field F is replaced by an effective field F e depending linearly on the electron temperature T n.In case of carrier heating (i.e., T n > TI where T/ is the lattice temperature), the hot-electron injection current density can be written as Jinj --qnA j['E+AF E lS exp ( ZE3 e.5 ) dE (7)   where E B denotes the height of the barrier at the Si-SiO 2 interface; AEB is the increase in the barrier height due to the oxide-field reversal which may occur near the drain if the drain voltage exceeds the floating-gate voltage; n is the electron concentration, and the values of the constants A and ) are deter- mined by fitting experimental data.Eq. ( 7) is then multiplied by a damping term, taking into account the electron scattering in the image force potential well into the SiO 2 layer [8].The model has been imple- mented as a post processor, due to the fact that the gate current is normally orders of magnitude lower than the other currents of interest.The above model was then improved and extended to holes [9] in anal- ogy with a substrate current model based on carrier temperature 10].
Figure 9 shows a comparison between simulation and experiment.An obvious shortcoming of the model is that the simulated threshold-voltage shift is sensibly larger than the measured one at large values of time, since the former curve does not account for the sudden change in the Fig. 10 slope which characterizes the experimental data.This is corroborated even fur- ther by figure 10, where simulated and experimental current data are reported in a linear plot versus the floating-gate voltage.Figure 10 exhibits a smoother drop of the simulated current as the floating-gate volt- age decreases.
Such a shortcoming has been overcome by assum- ing a modified functional form for the electron-distri-  bution function [11,12].A linear combination of two exponentials has been posited f (E) C [exp ( ZaE3 Tn.5 ) + Coexp ( -bE3 (8)   in order to overcome the assumpion that the electrons belong to a single nonparabolic conduction band.The resulting gate current model has the decided advan- tage to fit correctly the observed shoulder in the measured electron gate current at the expense of a larger number of fitting parametes.This is shown in figure 11, where the experimental and simulated gate currents turn out to be in a much better agreement than before 11 ].The computational workload involved in the simu- lation of the programming characteristics of realistic three-dimensional structures is rather heavy, primarily due to the size of the mesh and the transient nature of the effect under investigation.This seems to be one of the most important hurdles to the acceptance of three-dimensional device simulators in an engineering environment.

ERASING OPERATION
Cell erasing is performed by keeping the drain contact floating, whilst quickly ramping a moderate positive bias on the source contact and a large negative bias on the gate contact.In this situation, a large number of carriers tunnel through the Poly-Si-SiO 2 potential barrier.For emisson into a dielectric, using a free-electron gas model for polycrystalline silicon, and the WKB approximation for the tunneling proba- bility, and taking the effect of the image force on the barrier into account, the following expression for the Fowler-Nordheim tunneling current density is obtained 3, 4] q3mF2 rtckT 8rthm;x t2(y) sin(rtckT) exp ( 4v/2m,;x)3 ) 3qhF v(y) (9) where  FIGURE 11 Sirnulated and measured gate current against floating-gate voltage for W 0.7 gm, Vcg 11 V and Vd 6.0 V F is the electric field, m2x is the effective mass of the electron in the forbidden gap of the oxide, repre- sents the Poly-Si-SiO 2 barrier height, and is the normalized image-force barrier lowering (y A)/).Both v(y) and t(y) are slowly varying func- tions of y their expressions can be found in [13].The other quantitites have their usual meaning.A simplifi- cation is often made b.ecause -2 varies only between 1.00 and 0.81 and the correction due to the multiplicative term rtckT/sin(rtckT) is--for the ordinary tem- peratures; the following expression is obtained q3mF2 ( a v/2mxd3 ) JFN 8rthm;x---exp 3qhF v(y) (10) Only if the image-force correction is completely neglected, it is possible to derive from Eq. ( 10 where ) is expressed in eV.The tunneling current flows almost uniformly across the thin oxide, and no important edge effect was detected.For this reason, one can argue that the injection through the oxide depends on the channel width via the injection area and the coupling coefficient ag.The former effect increases the tunnel current at increasing channel widths, whereas the latter causes a variation in the opposite direction.In figure 12, measured erasing curves at different channel widths are represented: the narrower the cell, the faster the discharging process.Hence, the geometrical effect mentioned above has a strong impact on the erasing process.Note that the largest cell is markedly slower than the others.Coming now to the simulations, it was found that on a given structure it is possible to fit the erasing curve at different source voltages, as shown in figures 13, 14 and 15.However it is not possible to do so with a unique set of fitting parameters (i.e., with the same value of ) for all the curves at different channel widths, as shown in figure 16.Since there are no 0.0 w 0.5 gm w =0.7 lxm -0.5 --" w 0.9 Is   The two simulated curves refer to @ 2.83 eV and 2.95 eV physical reasons to believe that the barrier height is different in cells having different widths and, on the other hand, the values obtained for are different and lower than expected (i.e., it is necessary to enhance the field emission into the oxide), such a lack of agreement seems to indicate the existence of oxide thinning at the end of the bird's beak of the cell for the smaller cells.This effect leads to a local enhancement of the tunneling field which, due to the exponential dependence of the tunneling current upon the field, may cause a strong enhancement of the erasing char- acteristics.
In principle, it would be possible to use a unique set of values for AFN and BFN by introducing a smaller effective-oxide thickness.Once the parameters are fixed for the largest cell, they are expected to scale as AFNY 2 and BEN for the narrower cells, where "f=tox/teff=Fefl/Fox is determined by fitting the experimental curves.Of course, the above reasoning is based on the assumption that the oxide thinning affects more sensibly the narrower cells than the largest one, as it seems.
When considering the amount of oxide thinning, it turns out that the above discrepancies are fully explained by assuming a thinning of 0.5 nm, i.e. about one atomic layer.Thus, the extreme sensitiviy of the tunneling current upon the oxide thickness makes it very hard to accurately predict the erasing behavior of the cell.Moreover, the uniformity of the process can hardly be controlled within the above dimensional limits, so that the simulation is within the spreading of the experimental results.

CONCLUSION
In this paper we examine some of the problems related with device simulation in 3-D, and address the major deficiencies which still prevent a wide accept- ance of 3-D simulation in an engineering design envi- ronment for process and device optimization.More specifically, we identify structure definition, mesh generation and CPU time as the most important areas where a substantial effort is needed in order to make Biographies Giorgio Baccarani received his Dr. Ing. degree in Electrical Engineering in 1967 and his Dr. degree in Physics in 1969, both from the University of Bologna, Italy, where he is currently Full Professor of "Dig- ital-System Electronics".His research interests include numerical-device simulation, integrated-circuit design, analog and digital architectures for image processing, character and voice recognition, neural networks and fuzzy systems.He is author or coauthor of over 100 papers published on international scien- tific journals, and editor of two books.He is senior member of IEEE, and corresponding member of the Academy of Sciences in Bologna.
Massimo Rudan received a degree in Electrical Engineering in 1973 and a degree in Physics in 1976, both from the University of Bologna, Italy.He joined the Department of Electronics (DEIS) of the Univer- sity of Bologna in 1975, where he investigated the physical properties of the MOS structures and the problems of analytical modelling of semiconductor devices.Since 1983 he has been working in a group involved in numerical analysis of semiconductor devices, acting as Task leader in a number of EEC-supported Projects in the area of CAD for VLSI.In 1986 he has been a visiting scientist, on a one-year assignment, at the IBM Thomas J. Watson Research Center at Yorktown Heights, NY.In 1990, he was appointed Full Professor of Microelectronics at the University of Bologna.
Martino Lorenzini received his degree in Electri- cal Engineering from the University of Bologna, Italy,  in 1993.Currently, he is working towards his Ph.D. degree at the Department of Electronics (DEIS) of the same University in the field of the numerical simula-

FIGURE 3
FIGURE 3 Simulated and experimental turn-on characteristics for the three EEPROM cells with different channel widths.Vds O. V

FIGURE 4
FIGURE 4 Simulated and experimental output characteristics for the three EEPROM cells with different channel widths.Vcg 5 V

FIGURE 5 FIGURE 6
FIGURE 5 Simulated turn-on characteristics for different values of the floating-gate charge

FIGURE 8
FIGURE 8 Gate current against floating-gate voltage as extracted from the experimental data of figure7.Vcg11 V

FIGURE 9 FIGURE 10
FIGURE 9 Experimental and simulated threshold-voltage shift against programming time for W= 0.7 gm, Vcg 11 V and V d4.5 V

10 FIGURE 12
FIGURE 12 Measured threshold-voltage shift vs erasing time for the three different channel widths.Vcg -8 V, V 4 V, I d 0 A

FIGURE 13
FIGURE13 Measured and simulated threshold-voltage shift vs erasing time for three different values of the source voltage.W 0.9 gm, gcg ---8 V, (I)= 2.95 eV

FIGURE 14 10 FIGURE 15 10 FIGURE 16
FIGURE 14 Measured and simulated threshold-voltage shift vs erasing time for three different values of the source voltage.W 0.7 gm, Vcg =-8 V, 2.83 eV