Generally Applicable Degradation Model for Silicon MOS Devices



The main cause of operational degradation in MOS devices is believed to be due to the buildup of charge at the Silicon-Oxide interface. This leads to reduced saturation currents and threshold voltage shifts in MOSFET devices. Physics-based models of the degradation process typically consider the breaking of Si-H bonds (depassivation) at the Silicon-Oxide interface to be the main cause of the operational degradation. A new general model of Si-H bond breaking has recently been included in Atlas, adding to the Silvaco TCAD portfolio of degradation models[1].

This article presents the theory of the new model, and describes its implementation in Atlas. The model is then applied to a simple MOSFET to illustrate the features of the model. Finally it is applied to model a realistic MOSFET for which experimental degradation data are available. It is able to simulate reasonably well the unusual behavior of the degradation as a function of stressing time.


General Framework Model

The model is based on a study of Si-H trap dynamics in which three bond breaking mechanisms are considered[2]. The first mechanism occurs at high electric field, which distorts the bond and reduces the amount of thermal energy needed to break the bond. The second mechanism involves a high energy (hot) carrier breaking the bond with a single interaction, and the third involves many lower energy (cold) carriers exciting a vibrational mode to higher and higher energies until the bond breaks. These two different carrier mediated processes are necessary in order to explain some aspects of Hot Carrier Degradation [3]. Along with that work, we refer to the hot carrier process as single-particle (SP) and the cold carrier process as multi-particle (MP). First we describe the single-particle process. The time evolution of the interface charge is assumed to be of the form


for electrons, where N(r,t) is negative interface charge density, and


for holes, where P(r,t) is positive interface charge density. The quantities and represent the saturated values of negative and positive interface charge density associated with the SP process, and the time is t in seconds. The reaction rate for this process at position r is given by


where f(E, r) is the anti-symmetric part of the carrier distribution function, g(E) is the density of states and ug is the group velocity. For electrons, the function σ (E, Esp) is defined for EEsp, where


where the Boltzmann energy KbT acts as an energy scale. This is known as a soft-threshold, as introduced by Keldysh in the context of impact ionization rate calculations. Therefore only electrons with an energy of more than Esp contribute to this integral. Analagously, the function is defined for holes as


Equation [3] is often referred to as an acceleration integral in the literature, although its units are s-1.

The MP process involves gradual excitation of the bending vibrational quantum states of the bond by less energetic carriers, followed by a thermal excitation from the highest bound state to the transport state of the Hydrogen. This thermal emission occurs over a barrier of height Eemi eV, with an attempt frequency of νemi Hz, giving an emission rate of


where T is the lattice temperature. There is also the reverse process for repassivation of the bond, where the hydrogen overcomes a barrier of height Epass to become bonded again. The overall repassivation rate is


The excitation of the bond by numerous cold carriers can be described by a set of coupled differential equations describing the occupation density of each level [2]. Entering these equations as parameters are Pu and Pd which are the probabilities of transition to the next higher vibrational state and the next lower vibrational state respectively. These are modelled by the expressions




where νphon is an attempt frequency and hω is the vibrational mode energy. The acceleration integral is


where f(E,r) is the anti-symmetric part of the carrier distribution function. The cross-section σ (E,Emp) is given by the expression


Because these processes depend on cold carriers, the threshold energies are less than the threshold energies in the SP process. After some mathematical manipulation and simplification, the density of traps created by the MP process is given by


for electrons, where N(r,t) is negative interface charge density, and


for holes, where P(r, t) is positive interface charge density. Nl is the number of bending mode vibrational levels in the Si-H bond. Analysis of equations (12) and (13) shows that the time evolution depends only on the emission rate. The saturation level depends on the unpassivated bond densities N and N , ratio of depassivation rate to passivation rate and the ratio of Pu to Pd, raised to the power of Nl. This last ratio will be very small in the absence of a significant acceleration integral, as it will be a Boltzmann factor with energy of approximately the binding energy of the ground state. From equations (8) and (9) it is seen that if K (MP)(r) is greater than the attempt frequency, the ratio of Pu to Pd is approximately one. The spatial distribution of traps depends, therefore, in a very non-linear manner on the acceleration integral.

The third component of the general framework model is a field-enhanced thermal degradation, which is modelled as


where Ktherm is an attempt frequency and Eb is the Field dependent Si-H bond energy. It has the same time dependence as the SP process and so it is simply added to K (SP)(r) in the calculation of defects after stressing time t.

Many of the model parameters can be set on the DEGRADATION, MATERIAL or MODELS statements. For example NTA.SP, NTA.MP, NTD.SP and NTD.MP on the DEGRADATION statement specify N , N , N , N respectively


Calculation of the Carrier Distribution Function

Equations (3) and (10) require the anti-symmetric part of the carrier distribution function. The capability to solve the Boltzmann Transport Equation (BTE) for the zeroth and first order terms in a Spherical Harmonic expansion of the carrier distribution function has recently been added to Atlas. The first order term is anti-symmetric and is used in equations (3) and (10). In a similar model Starkov et al [3] used Monte Carlo simulations to estimate the carrier distribution function. Reggiani et al [4] used an analytical formulation for the carrier distribution function, with parameters derived from the Spherical harmonic expansion solution to the Boltzmann transport equation. This approximation was made to improve calculation speed. The Atlas implementation of the BTE solver is sufficiently rapid that a further approximation of the carrier distribution function is not necessary. The BTE solver is based on the formulation of Ventura et al [5]. The equation for the zeroth order expansion, fo, of the carrier distribution function is


where E is energy in eV, g(E) is the density of states in m-3-1, F is field in V/m, τ(E) is a scattering lifetime in seconds, ug is the group velocity in m/s, cop is optical phonon scattering coefficient in m3J/s, Nop is the optical phonon occupation number and the optical phonon energy is ω in eV. N is the optical phonon occupation number plus one, simplified as follows


where Kb is Boltzmanns constant and T1 is the lattice temperature. The first order expansion, f1, is then obtained from


The lifetime (E) is derived from the carrier scattering mechanisms. Scattering mechanisms which are included by default are optical phonon scattering, acoustic phonon scattering and ionized impurity scattering. Impact ionization scattering can also be included if required. Quantities such as carrier density, drift velocity and energy can be calculated from the carrier distribution functions. For example, the drift velocities as a function of homogeneous field are shown in Figure 1 for electrons and Figure 2 for holes. Results are shown for three different values of dopant concentration.

Figure 1. Homogeneous velocity field curves for electrons.


Figure 2. Homogeneous velocity field curves for holes.


To initialize Atlas for solving the BTE, the flags BTE.PP.E for electrons and BTE.PP.H for holes must be set on the MODELS statement. After the BTE has been solved for a specific bias set, Atlas includes the acceleration integrals when it saves the structure to file. See the Atlas manual[1] for more details on the BTE solver.


Implementation of the General Framework Model

An Atlas device is biased to the stressing configuration using the drift-diffusion or energy-balance models. A SOLVE statement with the flags DEVDEG.GF.E for electrons and DEVDEG.GF.H for holes will solve the Boltzmann transport equation. Up to 10 degradation times can be simulated using the parameters TD1 .. TD10 on the SOLVE statement. The interface charge densities are calculated using equations (1), (2), (12) and (13) for each requested degradation time, and the results are written to a structure file. For example, the Atlas statement

SOLVE DEVDEG.GF.E TD1=1.0e-2 TD2=1.0e-1 TD3=1.0 TD4=10.0 TD5=1.0e2 OUTFILE=simstd.str

will result in files


being written out, each having an interface charge density corresponding to the simulated degradation time.


Example: Simple MOSFET

The first example is for the MOSFET structure shown in Figure 3. Each of the three different degradation models are looked at in turn for the case of electrons in this device. The MP Keldysh cross-section σemp,0 was set to zero and the SP Keldysh cross-section was set to be 1.0 × 10-22 cm2, with a threshold energy of 2.2 eV. The saturated dangling bond density Nasp was set to 4 × 1012 cm-2. With a gate bias of 2 V, a BTE solution was obtained at a drain voltage of 2 V and also at a drain Voltage of 4 V. The acceleration integral for the SP process is shown in Figure 4 for 2V Drain bias and Figure 5 for 4V Drain bias. At 4V drain bias it is many orders of magnitude larger than at 2V drain bias. In Figure 6 the first order component of the electron distribution function is plotted, at the node where the SP acceleration integral is a maximum. The electron distribution function at 4V drain bias is much larger at higher energies than the equivalent distribution at 2V drain bias. The energy threshold in the calculation of acceleration integral is 2.2 eV, and clearly the electron distribution function at 4V drain bias is much larger above this energy.


Figure 3. Example structure.


Figure 4. SP Acceleration integral at 2V drain bias (logarithmic scale).


Figure 5. SP Acceleration integral at 4V drain bias (logarithmic scale).


Figure 6. First order component of electron distribution function.

The simulation was performed with degradation times of 10 milliseconds, 100 milliseconds, 1 second, 10 seconds and 100 seconds. The threshold Voltage after each simulation time was calculated from the Gate bias required to achieve a specified drain current, and the threshold Voltage shifts calculated. At 2V drain bias there was negligible threshold voltage shift, and so the calculation was performed with drain biases of 3V and 4V, with the resulting shifts being shown in Figure 7.

Figure 7. Threshold Voltage shifts due to SP process as a function of stressing time.


In order to study the MP process in isolation the SP Keldysh cross-section σesp,0 was set to zero and the MP Keldysh cross-section σemp,0 was set to be 1.0 × 10-13 cm2, with default values for other parameters, including a threshold energy of 1 eV. The default parameters give a value of Pemi of approximately 0.036 /second, and so at 100 seconds the time evolution will be essentially complete. The saturated dangling bond density N was set to 1 × 1013cm-2. With a gate bias of 2 V, a BTE solution was obtained at a drain voltage of 2 V and also at a drain Voltage of 4 V. The MP Acceleration integral is shown in Figures 8 and 9 for these two bias points. There is less difference between the two cases than for the SP process, due to the lower threshold energy. From Figure 6, it is shown that the distribution functions are very similar up to about 1.5 eV, and consequently give similar contributions to the MP acceleration integrals in this range. Because of the lower threshold energy and higher value of cross-section the values of the MP acceleration integral are much higher than the SP acceleration integral under the same conditions. High values are required to give a sizeable value of interface charge density, and in Figure 10 it is seen that the maximum of interface charge density is at the same position as the maximum acceleration integral.

Figure 8. MP Acceleration integral at 2V drain bias (logarithmic scale).


Figure 9. MP Acceleration integral at 4V drain bias (logarithmic scale).


Figure 10. Trapped interface charge density from MP process.

The simulation was performed with the same degradation times as before, and for drain biases of 3 V and 4V. The resulting shifts in threshold voltage are shown in Figure 11. At a drain bias of 2 V the shifts were negligible.

The final mechanism to consider is the field-enhanced thermal degradation. The cross-sections σesp,0 and σemp,0 were set to zero, and the rate Ktherm was changed from its default value of 0 to be 1 × 1012s-1. The saturated dangling bond density Nasp was set to 4 × 1012cm-2. The gate was biased to 12 V with a drain bias of 0.01 V and the simulation carried out for the degradation lifetimes as above. The interface charge density, shown in Figure 12 after 100 s of stressing, is much more uniform than in the case of either SP or MP process degradation. The threshold voltage shifts typical of this process are shown in Figure 13.

Figure 11. Threshold Voltage shifts due to MP process as a function of stressing time.


Figure 12. Trapped interface charge density from thermal process.


Figure 13. Threshold Voltage shifts due to field-enhanced thermal process as a function of stressing time.


Fitting to Experimental Data

In an actual MOS device all of the three aforementioned degradation mechanisms may be important. In this section results from the model are used to analyze experimental data for a p-channel MOSFET. The device structure is shown in Figure 14 and the stressing biases are gate bias set to -2.1 Volts and drain bias set to -5.5 Volts, with all other contacts grounded. At these biases impact ionization is significant and so degradation by both electrons and holes is important. Impact ionization scattering can be included in the Boltzmann transport solver by specifying BTE.IMPACT on the MATERIAL statement. The main degradation metric used was the change in current in the linear regime, at a fixed gate bias. This shows an enhancement at short stressing time and a decrease at longer simulation time, as shown in Figure 15. One possible interpretation is that some interface acceptor traps are created on a very short time scale with a larger contribution from interface donor traps occurring over a longer stressing timescale.


Figure 14. p-MOSFET with net doping shown.


Figure 15. Experimental degradation data (Linear drain current).


The device stressing was simulated by using the Boltmann transport equation solver for both electrons and holes with simulation times of 1,2,5,10,20,50,100,200,500 and 800 minutes respectively. The trap densities were set as follows


and the emission and passivation parameters were set as follows


which result in a lifetime associated with the MP processes of approximately 50 seconds in the simulation. Other MP process parameters were


resulting in an MP electron integral having a maximum value of over 1013/s, and a saturated acceptor charge density along a 0.06 microns length of the device. This gives the initial enhancement in the current as the negative interface charge is created, which persists until approximately 10 minutes.

Figure 16. SP process acceleration integral.

The time evolution of the donor traps depends on Khf (SP)(r) and this quantity is shown in Figure 16. The Keldysh parameters used were


and as can be seen from the figure this produces a maximum value of Khf (SP)(r) of about 15 s-1. The maximum value is away from the interface, and on the interface the maximum value is of the order of 1 s-1, but with a significant part of the interface having values down to 10-5s-1, which match the maximum timescale of the degradation stressing. Figure 17 shows the evolution with stressing time of the interface charge, along a part of the interface. The positive interface charge generated then reduces the drain current at -5 V, with the current reducing with increased stress time. The percentage change in current is shown as a function of stressing time in Figure 18. This simulation shows good qualitative agreement with experiment.

Figure 17. Stress time evolution of donor interface charge.


Figure 18. Simulated degradation data (Linear drain current).



  1. Atlas User’s Manual, Silvaco, (2014).
  2. C. Guerin, V. Huard and A. Bravaix,’ General framework about defect creation at the Si/SiO2 interface’, J.Appl.Phys., Vol. 105, (2009), 114513.
  3. I. Starkov,S. Tyaginov, H. Enichlmair, J.Cervenka, C.Jungemann, S. Carniello, J.M.Park, H.Ceric and T.Grasser,’Hot-carrier degradation caused interface state profile - Simulation versus experiment’, J.Vac. Sci. Technol B, Vol. 29, 01AB09-1/8 (2011).
  4. S.Reggiani,G.Barone,S.Poli,E.Gnani,A.Gnudi,G.Baccarani, M-Y.,Chuang, W.Tian, R.Wise,’TCAD Simulation of Hot-Carrier and Thermal Degradation in STI-LDMOS Transistors’,IEEE Trans. Elec. Dev. Vol. 60, No.2 , (2013), pp.691-698.
  5. D.Ventura, A.Gnudi, G.Baccarani, ‘A deterministic approach to the solution of the BTE in semiconductors’, Rivista del Nuovo Cimento, Vol. 18, No. 6 pp. 1-32, (1995).