3.2 From Rates to Degradation Curves

The last section has been devoted to the possible physical explanations for the charge trapping process in NBTI. The involved mechanisms, such as quantum mechanical tunneling for instance, are characterized by their stochastic nature. This means that one must deal with probabilities instead of pre-determined transition times for the trapping events. Such problems can be best handled within the framework of homogeneous continuous-time Markov chain theory [131], which rests on the assumption that the transition rates do not depend on the past of the investigated system.

The continuous time Markov chain theory presumes a set of discrete states S ∈ {1,...,N } :

Π(t) = (π1(t),π2(t),...,πN (t))                (3.13)
The πi(t)  are the time-dependent occupation probabilities, which fulfill the normalization condition
   πi(t) = 1 .                       (3.14)
This equation states that the probability of finding the system in one of the states i ∈ S  equals unity. The time-dependent probability pij  for a transition from state i  to state j  can be formulated as
pij(Δt)  = P{X (t+ Δt) = j|X(t) = i} = rijΔt + O(Δt) ,         (3.15)
pii(Δt)  = P{X (t+ Δt) = i|X (t) = i} = 1-   rijΔt + O (Δt ) ,     (3.16)
where rij  denotes the transition rate belonging to pij  . The time evolution of the system can be described by the first-order partial differential equation
∂π (t)  =    ∑  π (t)r    -     ∑  π(t)r     .       (3.17)
 ti         j⁄=i j   ji        j⁄=i i  ij
            ◟---◝◜---◞         ◟--◝◜---◞
         transitions into state i transitions out of state i
The above equation is referred to as the Master or rate equation and controls the time-dependent occupation probabilities. The first term on the right-hand side represents the transitions Tj→i  while the second term stands for the opposite direction. With the definition
rii = -   j⁄=irij                      (3.18)
the rate equation (3.17) can rewritten as
∂tΠ(t) = Π (t) ⋅R  ,                    (3.19)
where R  denotes the matrix belonging to the elements rij  . The solution of equation (3.19) is given by
Π (t)  = Π(0)⋅exp(RΔt )                  (3.20)
and can be expanded into a truncated Taylor series.
Π(t)  ≈ Π(0)⋅(1+ RΔt )                  (3.21)
Assuming detailed balance, the time derivative of Π(t)  vanishes and the equilibrium occupation probabilities Πeq  must satisfy
Πeq ⋅R = 0 .                        (3.22)

In the case of charge trapping, a single defect is represented by one vector Π(t)  . When considering a large ensemble of defects, the stochastic behavior vanishes and πi(t)  and Π(t)  can be replaced by their corresponding expectation values fi(t)  and F(t)  , respectively. Then the equations (3.19) and (3.22) read

∂tF (t)  = F (t)⋅R                      (3.23)
Feq ⋅R = 0 ,                       (3.24)

When applying Markov theory to charge trapping in NBTI, each state must be assigned to a certain configuration and a certain charge state of a defect. The rates rij  linking the states i  and j  must be related to certain defect transitions, such as charge transfer reactions or thermally activated rearrangements of the defect structure. The trapping dynamics are then governed by equation (3.23), which must be solved as a function of time. The details of the applied numerical procedures are outlined in Fig. 3.2.


Figure 3.2: A flow chart diagram for the calculation of stress/relaxation curves. At first, the occupations for equilibrium conditions (VG = Veq  , T = Teq  ) are calculated in steps A, B, and C. During one stress or relaxation cycle, the band edges (D) and the transition rates (E) must be determined at first. Then the corresponding the time-dependent occupations (F) and the resulting threshold voltage shift (G,F) are evaluated at each time step tk  . This complete cycle must be performed for each stress or relaxation phase.

Before the real device is subject to stress, the device is assumed to be in equilibrium. The corresponding band diagram is computed by a P/SP-solver (A) for the equilibrium conditions VG = Veq  and T = Teq  . Based on this information, the transition rates req,ij,m  (B) can be evaluated for each defect m  . With the rates req,ij,m  at hand, the equilibrium occupation probabilities feq,i,m  are calculated using the equation (3.24) and must be subsequently stored (C). When stress sets in, the gate voltage VG  and the temperature T  are changed to Vs  and Ts  , respectively. Since this alters the band diagram and in consequence the transition rates rs,ij,m  , the S/SP-solver step (D) and the calculation of the rates (E) must be repeated. Then the implicit Euler method is employed for the numerical time integration of equation (3.23). This iteration scheme must be continued until the end of the stress time ts  has been reached. At each tk  , the change in the defect occupancies Δfi,m(tk+1 )  (G) and the corresponding threshold voltage shift ΔVth(tk+1)  (H) are computed. With the beginning of the relaxation phase, VG  and T  are modified again and the whole iteration loop including the steps D, E, F, G, and H must be repeated. This can be continued for several stress/relaxation cycles with different stress/relaxation conditions. The steps E, F, and G of numerical procedure in Fig. 3.2 must be repeated for each individual defect, where the calculation of the rates requires most of the computation time.

For the evaluation of ΔVth  , the charge sheet approximation is employed.

             (       )
ΔV    = q0tox- 1 - xm- Δf                    (3.25)
  th,m    εrε0      tox     i0,m
xm  and tox  are the trap depth and the thickness of the dielectric, respectively. The above expression gives the threshold voltage shift ΔVth,m  due to a single trap m  in the ‘charged’ state i0  . It is noted that the term in the parentheses of equation (3.25) accounts for the trap depth dependence of ΔVth,m  according to which traps closer to the interface have a larger impact on the threshold voltage. However, there also exists a strong dependence on the lateral location of the traps, totally neglected in the above equation. Recall in this context that the Vth  steps in the TDDS measurements have also been caused by single charged traps. The wide variations of these steps have been ascribed to the position of the traps relative to the current percolation path and follow an exponential distribution according to an investigation carried out by Kazcer et al.  [50]. This issue has been intensively studied under the name random dopant fluctuations [13213313450] and gained large interest due to its serious influence on the lifetime projection. For instance, the charging of one defect can even produce step heights in the threshold voltage beyond 30mV  , which already violates typical lifetime criteria. By contrast, the charge sheet approximation rests on the assumption that a trapped charge is distributed over the whole plane parallel to the interface. Therefore, variations in the threshold voltage shift due to random dopant fluctuations remain unconsidered in this approximation.

Small-area devices often contain only a handful of defects. Then the steps B, C, E, F, G, and H can be performed at computationally feasible costs. But since the calculation time increases with the number of defects, the simulation of large-area devices can become time-consuming. In order to reduce computation time for these simulations, a certain number of defects with similar properties are grouped together and replaced by one representative trap. This approach implies that the channel area covers several defects with similar properties including the trap level and the spatial location among others1. Furthermore, it is important to note here that the value of trap density is often quite inaccurate and can differ by some orders of magnitude. Therefore, the number of traps must be treated as a variable in the simulations. This means that the calculated degradation curves can be scaled to the experimental data in order to achieve reasonable agreement of the simulations with the measurements.