The Physics of Non–Equilibrium Reliability Phenomena

Chapter 4 Theoretical Approaches & Models

Electrochemical reactions triggered by a non–equilibrium carrier ensemble fundamentally define the degradation mechanisms discussed in this work. In particular heated carriers can provoke Si–H bond excitation at the Si/SiO2 interface or influence the dynamics of defects in the SiO${}_2$ via charge trapping. This Chapter addresses the key concepts introduced in Chapter 2 and establishes the theoretical descriptions relevant for the model applications in the last Chapter 5. The methods presented here are not limited to the Si/SiO2 material system and have been used in a much broader field to understand the interactions within emerging material combinations and related phenomena.

4.1 Hot–Carriers & Interface Defects

This part of the thesis focuses on the theoretical formulation of transition rates ${\mathrm{\Gamma }}_{i,f}$ between individual vibrational eigenstates of the Si–H ground state potential $V(q)$ ultimately leading to bond breakage and the creation of interface defects. The cornerstone is a consistent description of the interaction between a Si–H bond and its environment, including energetic, hot charge carriers. It will be shown that two major components determine the kinetics: Dynamics related to vibrational relaxation and inelastic scattering. The latter one can be split into dipole–induced transitions and excitations due to a transient population of a resonance electronic state.

4.1.1 System–Bath Interactions

Ideally, all dissipative and inelastic couplings leading to changes of the states in both subsystems¹ are taken into account. However, due to the macroscopic nature of the bath (the Si/SiO2 environment), an explicit treatment of the system’s (the Si–H bond) dynamics is desired, while the bath entity is implicitly described. Therefore, focussing solely on the dynamics and properties of the system, it is convenient to separate the total system into a system and bath entity

$$\begin{array}{}\text{(4.1)}& \hat{H}={\hat{H}}_{\mathrm{S}}+{\hat{H}}_{\mathrm{B}}+{\hat{V}}_{\mathrm{SB}}={\hat{H}}_0+{\hat{V}}_{\mathrm{SB}},\end{array}$$

where ${\hat{V}}_{\mathrm{SB}}$ is the system-bath coupling operator. Open–system–density–matrix theory is a well suited approach for the problem at hand where the evolution of the system is described by the reduced density operator ${\hat{\rho }}_{\mathrm{S}}$, satisfying the open-system Liouville-von Neumann (LvN) equation [268]

$$\begin{array}{}\text{(4.2)}& i\hslash \frac{\partial }{\partial t}{\hat{\rho }}_{\mathrm{S}}=i\hslash \mathcal{L}{\hat{\rho }}_{\mathrm{S}}(t)=[{\hat{H}}_0,{\hat{\rho }}_{\mathrm{S}}(t)]+{\int }_0^{t}d\tau {\mathcal{L}}_{\mathrm{D}}(t,\tau ){\hat{\rho }}_{\mathrm{S}}(\tau ).\end{array}$$

Within this formulation the system–bath interaction is included in the dissipative Liouvillian (superoperator) ${\mathcal{L}}_{\mathrm{D}}$. However, for all practical applications the dissipative part of (4.2) needs to be approximated. Thereby, the Lindblad semigroup functional formulation [269] provides the most general type of a Markovian and time–homogeneous master equation preserving the laws of quantum mechanics². Within this approach the bath modes have been already traced out resulting in equations of motion solely for the system part of the problem. The reduced density operator ${\hat{\rho }}_{\mathrm{S}}$ can be written as

$$\begin{array}{}\text{(4.3)}& \frac{\partial }{\partial t}\widehat{{\rho }_{\mathrm{S}}}(t)={\mathcal{L}}_{\mathrm{S}}{\hat{\rho }}_{\mathrm{S}}+{\mathcal{L}}_{\mathrm{D}}{\hat{\rho }}_{\mathrm{S}}=-\frac{i}{\hslash }[{\hat{H}}_{\mathrm{S}},{\hat{\rho }}_{\mathrm{S}}]+{\mathcal{L}}_{\mathrm{D}}{\hat{\rho }}_{\mathrm{S}},\end{array}$$

where ${\mathcal{L}}_{\mathrm{S}}$ and ${\mathcal{L}}_{\mathrm{D}}$ denote the system Hamiltonian and the dissipative superoperator. Furthermore, utilizing the eigenstate representation of ${\hat{\rho }}_{\mathrm{S}}$, such as the eigenfunctions $|{\varphi }_n⟩$ of the system Hamiltonian, ${\rho }_{nm}=⟨{\varphi }_n|\hat{\rho }|{\varphi }_m⟩$, allows the interpretation of the diagonal elements ${\rho }_{nn}$ as the level populations $P_n$ and the off–diagonal represents the coherences between two states.

The superoperator ${\mathcal{L}}_{\mathrm{D}}$ within the Lindbladian is given by

$$\begin{array}{}\text{(4.4)}& {\mathcal{L}}_{\mathrm{D}}{\hat{\rho }}_{\mathrm{S}}=\sum _k({\hat{C}}_k{\hat{\rho }}_{\mathrm{S}}{\hat{C}}_k^{†}-\frac{1}{2}{[{\hat{C}}_k^{†}{\hat{C}}_k,{\hat{\rho }}_{\mathrm{S}}]}_{+}),\end{array}$$

with ${\hat{C}}_k$ being the Lindblad operators, defining the nature of $k$ dissipative channels. They account for population transfer between eigenstates of the ground state PEC $V(q)$ of the system due to vibrational relaxation and inelastic scattering. The Lindblad operators can be expressed as

$$\begin{array}{}\text{(4.5)}& {\hat{C}}_k=\sqrt{{\mathrm{\Gamma }}_{i,f}}|{\varphi }_f⟩⟨{\varphi }_i|,\end{array}$$

which is the projection of an initial state onto a final state weighted with a rate ${\mathrm{\Gamma }}_{i,f}$. The individual transition rate depends on the type of perturbation as well as on the bond complex and its properties. A general formulation is given by FGR

$$\begin{array}{}\text{(4.6)}& {\mathrm{\Gamma }}_{i,f}=\frac{2\pi }{\hslash }\sum _{I,F}{f}_I(T)[1-f_F(T)]\cdot |⟨{\varphi }_f{\mathrm{\Phi }}_F|\hat{V}|{\varphi }_i{\mathrm{\Phi }}_F⟩{|}^2\delta (E_{i,I}-E_{f,F}).\end{array}$$

The summations run over all collective initial $I$ and final states $F$, weighted with the energy distribution function $f$ of carriers. Furthermore, the coupling operator $\hat{V}$ specifies the nature of the perturbation, and thus the interaction strength of $|{\varphi }_i⟩$ and $|{\varphi }_f⟩$.

In the following (4.6) is used to calculate transfer rates of localized Si–H eigenstates due to various mechanisms. Contributions originating from bond vibration-phonon coupling as well as from inelastic electron tunneling are included. The latter can either arise due to dipole scattering or resonance scattering. Eventually, the total rate ${\mathrm{\Gamma }}_{i,f}$ is the sum of all individual contributions:

$$\begin{array}{}\text{(4.7)}& {\mathrm{\Gamma }}_{i,f}^{\mathrm{tot}}={\mathrm{\Gamma }}_{i,f}^{\mathrm{vib}}+{\mathrm{\Gamma }}_{i,f}^{\mathrm{dip},\mathrm{inel}}+{\mathrm{\Gamma }}_{i,f}^{\mathrm{res},\mathrm{inel}}\end{array}$$

4.1.2 Model Potential

As was shown in Section 3.2 and in [MJJ3], Si–H bond breakage at the Si/SiO2 interface follows a trajectory which is a combination of the bond’s bending and stretching mode. The one–dimensional PEC, where the collective ionic motion has been mapped onto a single reaction coordinate, can be assessed directly from the DFT results. For all practical applications³, however, an analytic fit of the following form has been used:

$$\begin{array}{}\text{(4.8)}& V_{\mathrm{g}}(q)=v_0{q}^4-v_2{q}^2+\sigma \sqrt{\frac{v_0}{2v_2}}q,\end{array}$$

which represents an asymmetric double–well potential, see Fig. 4.1. The parameters $v_0$, $v_2$ and $\sigma$ are given in Table 4.1a. Overall the (reference) potential supports two minima at $q=\text{−}3.2a_0\sqrt{m}$ and $q=3.0a_0\sqrt{m}$ with a barrier of $2.7eV$ at $q=0$ and exhibits 32 bound states (below the barrier), localized either in the left or right well. Selected eigenstates of the ground state potential $V(q)$ are summarized in Table 4.1b.

Furthermore, also the characteristics of the excited PECs are captured by the same functional form as the ground state potential, see Sec. 3.4. The negatively (positively) charged profile possesses a substantially lowered transition barrier of $\sim 1.7eV$ ($\sim 2.1eV$) and a shifted equilibrium position by $\sim 0.9a_0\sqrt{m}$ ($\sim 0.4a_0\sqrt{m}$), see Fig. 4.1 and Tab. 4.1a. Note that the positive potential $V^{+}(q)$ additionally exhibits a shifted transition point.

The analytic fits presented here do not aim to provide the most accurate description of the DFT results, but to give a qualitatively correct description of the ground and excited state potential energy surfaces. Due to the amorphous nature of the Si/SiO2 system and the interface itself, an intrinsic distribution of the parameters such as the transition barrier and the minimum position as well as the resonance energy is to be expected, see Table 4.1b.

4.1.3 Dissipative Relaxation

The basic interaction between the Si–H bond and its environment can be described by the coupling of the bond potential to a bath of oscillators. The resulting vibrational relaxation rates ${\mathrm{\Gamma }}_{i,f}^{\mathrm{vib}}$, and thus the lifetimes ${\tau }_i$, strongly influence the bond dynamics. As already mentioned above, the reaction coordinate of Si–H breakage at a Si/SiO${}_2$ interface is a combination of the Si–H bending and Si–H stretching modes. Experimental and theoretical studies, however, have focused on the vibrational relaxation of various Si–H modes on a silicon surface. These studies revealed lifetimes on the order of ns and ps for the stretching and bending mode [270–273], respectively. However, for the problem at hand – a mixture of various Si–H modes and a substantially different phonon spectrum – these values can only be used to derive an educated guess for the vibrational relaxation.

Investigating the detailed coupling between the fundamental frequency of the ground state PEC $V(q)$ and the substrate is, therefore, an essential component in the presented framework⁴. Following a perturbation theory approach allows to calculate the system–bath coupling elements and ultimately the transition rates using Fermi’s Golden Rule, see [271, 272, 274]. By introducing a new set of coordinates one can separate the $3N$–dimensional coordinates $\mathbf{r}$ of the total system into system and bath degrees of freedom (DOF). The new DOF $\mathbf{q}$ and $\mathbf{Q}$ can be constructed using an orthogonal transformation

$$\begin{array}{}\text{(4.9)}& R_i=\sum _{j=1}^{3N}{O}_{i,j}{r}_j,\mathrm{with}\sum _{i=1}^{3N}{O}_{i,k}{O}_{i,l}={\delta }_{k,l},\end{array}$$

where $q_i=R_i(i=1\dots M)$ describes the system coordinates and $Q_i=R_{i+M}(i=1\dots 3N-M)$ corresponds to the bath modes. The total Hamiltonian can be decomposed in the usual system–bath form:

$$\begin{array}{}\text{(4.10)}& \hat{H}={\hat{H}}_S+{\hat{H}}_{SB}+{\hat{H}}_B,\end{array}$$

with the individual contributions of the system, bath and the system–bath coupling given by:

$$\begin{array}{}\text{(4.11)}& \begin{array}{rl}{\hat{H}}_S& =-\frac{{\hslash }^2}{2}\sum _{i=1}^M\frac{{\partial }^2}{\partial q_i^2}+V(\mathbf{q},\mathbf{Q}=0),\\ {\hat{H}}_{SB}& =\sum _{i=1}^{3N-M}{\lambda }_i(\mathbf{q})Q_i+\frac{1}{2}\sum _{i,j=1}^{3N-M}{\mathrm{\Lambda }}_{i,j}(\mathbf{q})Q_i{Q}_j,\\ {\hat{H}}_B& =\sum _{i=1}^{3N-M}(-\frac{{\hslash }^2}{2}\frac{{\partial }^2}{\partial Q_i^2}+\frac{1}{2}{\omega }_i^2{Q}_i^2).\end{array}\end{array}$$

The one– and two–phonon coupling functions ${\lambda }_i(\mathbf{q})$ and ${\mathrm{\Lambda }}_{i,j}(\mathbf{q})$ (corresponding to one– and two–phonon relaxations) can be obtained by calculating the derivative of $V(\mathbf{q},\mathbf{Q})$ as

$$\begin{array}{}\text{(4.12)}& \begin{array}{rl}{\lambda }_i(\mathbf{q})& =\frac{\partial V}{\partial Q_i}{|}_{\mathbf{q},\mathbf{Q}=0}=\sum _{k=1}^{3N}{O}_{M+i,k}\frac{\partial V}{\partial r_k}{|}_{\mathbf{q},\mathbf{Q}=0},\\ {\mathrm{\Lambda }}_{i,j}(\mathbf{q})& =\frac{{\partial }^{2}V}{\partial Q_i\partial Q_j}{|}_{\mathbf{q},\mathbf{Q}=0}=\sum _{k,l}^{3N}{O}_{M+i,k}{O}_{M+j,l}\frac{{\partial }^{2}V}{\partial r_k\partial r_l}{|}_{\mathbf{q},\mathbf{Q}=0}.\end{array}\end{array}$$

These couplings are used within Fermi’s Golden Rule (FGR) to compute the phonon–induced transition rates between system eigenstates ${\varphi }_i$ and ${\varphi }_f$ as

$$\begin{array}{}\text{(4.13)}& {\mathrm{\Gamma }}_{i,f}=\frac{2\pi }{\hslash }|⟨{\mathrm{\Phi }}_f|{\hat{H}}_{\mathrm{SB}}|{\mathrm{\Phi }}_i⟩{|}^2\delta (E_i-E_f),\end{array}$$

where one can choose the initial and final states ${\mathrm{\Phi }}_i$ and ${\mathrm{\Phi }}_f$ as product states of system and bath states which can be integrated out. Thermal averaging over the initial states and summing over the final bath states allows to account for the thermal population of bath oscillator at finite temperature ($T\ne 0K$) and the respective upward rates.

Following [270, 272, 273, 275], and in particular the derivation in [274], the one–phonon coupling term in (4.11) and (4.12) yields the following rate expressions

$$\begin{array}{}\text{(4.14)}& \begin{array}{rl}{\mathrm{\Gamma }}_{i,f}^{(1),\downarrow }& =\pi \sum _k|⟨{\varphi }_i|{\lambda }_k(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k⟩+1}{{\omega }_k}\delta ({\epsilon }_i-{\epsilon }_f-\hslash {\omega }_k)\\ {\mathrm{\Gamma }}_{i,f}^{(1),\uparrow }& =\pi \sum _k|⟨{\varphi }_i|{\lambda }_k(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k⟩}{{\omega }_k}\delta ({\epsilon }_i-{\epsilon }_f+\hslash {\omega }_k),\end{array}\end{array}$$

where $⟨n_k⟩=\sum _{n}n{p}_{n,k}(T)$ is the thermally averaged quantum number of the $k$th bath mode with $p_{n,k}(T)$ being the thermal population of its $n$th level. However, if the energy of the system mode is above ${\omega }_{\max }$, the highest phonon frequency, two–phonon terms become important. The corresponding rates for the simultaneous (de–) excitation of two different bath modes are given by

$$\begin{array}{}\text{(4.15)}& \begin{array}{rl}{\mathrm{\Gamma }}_{i,f}^{(2a),\downarrow }& =\frac{\pi \hslash }{8}\sum _{k\ne l}|⟨{\varphi }_i|{\mathrm{\Lambda }}_{k,l}(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k⟩+1}{{\omega }_k}\frac{⟨n_l⟩+1}{{\omega }_l}\delta ({\epsilon }_i-{\epsilon }_f-\hslash {\omega }_k-\hslash {\omega }_l),\\ {\mathrm{\Gamma }}_{i,f}^{(2a),\uparrow }& =\frac{\pi \hslash }{8}\sum _{k\ne l}|⟨{\varphi }_i|{\mathrm{\Lambda }}_{k,l}(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k⟩}{{\omega }_k}\frac{⟨n_l⟩}{{\omega }_l}\delta ({\epsilon }_i-{\epsilon }_f+\hslash {\omega }_k+\hslash {\omega }_l)\end{array}\end{array}$$

whereas the expression for two identical bath modes is

$$\begin{array}{}\text{(4.16)}& \begin{array}{rl}{\mathrm{\Gamma }}_{i,f}^{(2b),\downarrow }& =\frac{\pi \hslash }{8}\sum _k|⟨{\varphi }_i|{\mathrm{\Lambda }}_{k,k}(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k^2⟩+3⟨n_k^2⟩+2}{{\omega }_k^2}\delta ({\epsilon }_i-{\epsilon }_f-2\hslash {\omega }_k)\\ {\mathrm{\Gamma }}_{i,f}^{(2b),\uparrow }& =\frac{\pi \hslash }{8}\sum _k|⟨{\varphi }_i|{\mathrm{\Lambda }}_{k,k}(\mathbf{q})|{\varphi }_f⟩{|}^2\frac{⟨n_k^2⟩-⟨n_k^2⟩}{{\omega }_k^2}\delta ({\epsilon }_i-{\epsilon }_f+2\hslash {\omega }_k)\end{array}\end{array}$$

For the sake of completeness two–phonon processes also give rise to the simultaneous creation and annihilation of phonons, however, this contribution can be neglected [274]. For all practical applications at finite temperatures, the delta–functions in (4.14), (4.15) and (4.16) need to be replaced by functions with a finite broadening, e.g. Lorentzians

$$\begin{array}{}\text{(4.17)}& \delta (E)\to \frac{1}{\pi }\frac{\gamma }{{\gamma }^2+E^2},\end{array}$$

with the broadening parameter $\gamma$.

The studied system is again the Si(100)/$a$–SiO${}_2$ interface system containing 472 atoms for which the Si–H system mode has been identified in Chap. 3 and in [MJJ3] using DFT. However, using DFT to calculate the coupling terms (4.12) would have made the calculations prohibitively expensive since these couplings need to be calculated at various system coordinates $\mathit{q}$. Thus, all linear and quadratic coupling terms (4.14), (4.15) and (4.16) have been calculated using the classical ReaxFF force–field [197] implemented in the LAMMPS code [198]. A comparison of the phonon spectrum calculated with DFT and the classical force–field for the Si/SiO2 interface system, as well as for bulk crystalline Si and bulk SiO${}_2$ is given in Fig. 4.2. One can see that ReaxFF is able to qualitatively reproduce the spectra calculated with DFT, albeit the individual spectra seem to extend towards higher energies. Surprisingly, the agreement between both methods is best for the most complex Si/SiO${}_2$ model.

The vector ${\mathbf{e}}_q$, which serves as the degree of freedom for the system mode $q$ and defines its fundamental frequency, see Fig. 4.2, has been identified using the results presented in Sec. 3.2. Furthermore, for all subsequent calculations its 3$N$ components representing the system motion have been restricted to atoms with displacements of at least 5% of the maximum displacement, all the other entries were set to zero. This was done to avoid an artificially large coupling due to the finite cluster size. In order to ensure orthogonality between the system mode vector ${\mathbf{e}}_q$ and the environment modes ${\mathbf{e}}_Q$, 3$N$-1 random 3$N$ dimensional vectors ${\mathbf{u}}_i$ have been selected. Subsequently, the system $({\mathbf{e}}_q,{\mathbf{u}}_i)$ was orthogonalized and a $(3N-1)\times 3N$ transformation matrix $\mathbf{U}$ was constructed using the orthogonalized vectors ${\mathbf{u}}_i$. The matrix $\mathbf{U}$ has been further used to transform the $3N\times 3N$ Hessian $\mathbf{H}=({\partial }^{2}V/\partial \mathbf{x}\partial \mathbf{x}){|}_0$ into a constrained (a subspace orthogonal to the system mode) $(3N-1)\times (3N-1)$ Hessian given by

$$\begin{array}{}\text{(4.18)}& {\mathbf{H}}^{\prime }=\mathbf{U}\mathbf{H}{\mathbf{U}}^T,\end{array}$$

which upon diagonalization yields $3N-1$ constrained bath normal modes ${\mathbf{e}}_Q^{\prime }$, each of length $3N-1$. Using the inverse transformation

$$\begin{array}{}\text{(4.19)}& {\mathbf{e}}_Q={\mathbf{U}}^T{\mathbf{e}}_Q^{\prime },\end{array}$$

finally gives $3N-1$ bath mode vectors with a length of $3N$. This approach guarantees one anharmonic system mode, represented by ${\mathbf{e}}_q$, and $3N-1$ harmonic bath modes, ${\mathbf{e}}_Q$, all being orthogonal to each other, ${\mathbf{e}}_{Q_i}{\mathbf{e}}_{Q_j}={\delta }_{i,j}$ and ${\mathbf{e}}_{Q_i}{\mathbf{e}}_q=0$, see [274].

The calculations for all individual one– (${\mathrm{\Gamma }}^{(1)}$) and two–phonon transition (${\mathrm{\Gamma }}^{(2)}$) rates have been conducted along the anharmonic system DOF ${\mathbf{e}}_q$ at 31 displacements $\mathbf{q}=q{\mathbf{e}}_q$ within the interval $[-0.75,+0.75]a_0\sqrt{m}$. Such a variation of ${\mathbf{e}}_q$ is sufficient to cover at least the lowest eigenstates. Ultimately, the described procedure yields the vibrational lifetime of the fundamental transition which is given by:

$$\begin{array}{}\text{(4.20)}& {\tau }_1^{-1}=({\mathrm{\Gamma }}_{1,0}^{(1),\downarrow }+{\mathrm{\Gamma }}_{1,0}^{(2),\downarrow })-({\mathrm{\Gamma }}_{0,1}^{(1),\uparrow }+{\mathrm{\Gamma }}_{0,1}^{(2),\uparrow }).\end{array}$$

The final result together with the individual contributions is summarized in Table 4.2. One can clearly see that the total lifetime ${\tau }_1^{\mathrm{total}}$ is dominated by a two–phonon relaxation process. One–phonon transition rates are around one order of magnitude smaller, however, their origin is questionable. Such relaxations are actually energy–forbidden and only arise as a direct consequence of the finite broadening used to replace the delta–functions in (4.17). Therefore, these contributions should be completely discarded. Furthermore, the dependence of the rates on the broadening parameter $\gamma$ has been investigated as well. While ${\mathrm{\Gamma }}_{1,0}^{(1)}$ exhibits a strong dependence between $\gamma \in [1,10]{\mathrm{cm}}^{-1}$ (${\tau }_1^{(1)}=18.3ns-1.89ns$), ${\mathrm{\Gamma }}_{1,0}^{(2)}$ is stable against changes in $\gamma$ (${\tau }_1^{(2)}=152.9ps-97.5ps$). Since the latter is the dominant contributor, the overall solution is not very sensitive to the precise choice of the fitting parameter $\gamma$, provided $\gamma$ is chosen in a physically reasonable range.

Further insight into the dissipation mechanism can be gained by analyzing the two–phonon process in detail. Fig. 4.3 shows the individual contribution of phonon pairs $\{{\omega }_k,{\omega }_l\}$ within a 2D–histogram with a bin size of $5/\mathrm{c}\mathrm{m}\times 5/\mathrm{c}\mathrm{m}$. One can see a pronounced dissipation pathway via two similar phonons (${\omega }_k\sim {\omega }_l$). However, as an alternative route, two phonons which in total fit $\mathrm{\Delta }{E}_{1,0}$ can be seen as well in Fig. 4.3. Additionally, the lifetime ${\tau }_1$ and its temperature dependence is shown in the right panel of Fig. 4.3. The decreasing lifetime with increasing temperature is a consequence of temperature–dependent weight factors in the rate expressions based on Fermi’s Golden Rule and also temperature–dependent upwards rates ${\mathrm{\Gamma }}_{0,1}^{\uparrow }$. For the sake of completeness, Fig. 4.3 also illustrates the influence of the only free parameter within the above framework, $\gamma$, (4.13), as described above. Based on the results $\gamma =5/\mathrm{c}\mathrm{m}$ has been chosen as the reference here⁵, which yields ${\tau }_{0\mathrm{K}}=0.126ns$ (${\tau }_{300\mathrm{K}}=0.120ns$).

In practice, the model potential possesses overall 32 bound states for both wells, see Section 4.1.2, for which it would be necessary to calculate all individual relaxation rates ${\mathrm{\Gamma }}_{i,f}$, ideally using the approach described here. However, for the problem at hand, the next vibrational state already yields an energy difference to the vibrational ground state $\mathrm{\Delta }{E}_{2,0}>2\hslash {\omega }_{\max }$, which is beyond the two–phonon coupling included⁶. Therefore, for all practical applications an idealized harmonic model [101, 105, 270, 276] can be employed where the selection rule $\mathrm{\Delta }n=1$ applies. Such a harmonic model possesses a linear scaling law ${\mathrm{\Gamma }}_{i,i-1}^{\downarrow }=i{\mathrm{\Gamma }}_{1,0}^{\downarrow }$ and has been demonstrated to show reasonable agreement with accurate calculations even for anharmonic system–modes [105, 270, 276]. In order to estimate the introduced error and further motivate this approximation, the one– and two–phonon rates as well as the lifetimes for the $2\to 1$ transition have been calculated, see Table 4.2. Since the energy difference $\mathrm{\Delta }{E}_{2,1}$ is still outside the Si/SiO2 phonon band, ${\mathrm{\Gamma }}_{2,1}^{(2)}$ again dominates the resulting lifetime. Furthermore, it is very likely that the direct transition $2\to 0$ is less important, due to the fact that at least three bath modes are necessary to dissipate the energy⁷. Finally, one can see in Tab. 4.2 that the resulting lifetime of ${\tau }_2=46ps$ agrees well with the simple model prediction of ${\tau }_2={\tau }_1/2=61.5ps$.

Mapping the harmonic scaling law onto the double well model potential introduced above splits $V(q)$ into two harmonic oscillators where only neighbouring transitions within one well are allowed, whereas transitions between the oscillator are prohibited. The corresponding single–phonon scaling law for the vibrational lifetime of an eigenstate $|{\varphi }_n⟩$ is then given by

$$\begin{array}{}\text{(4.21)}& {\tau }_n=\frac{{\tau }_1}{n},\end{array}$$

with ${\tau }_1$ being the lifetime of the first eigenstate as calculated above.

A detailed comparison of the utilized methods including the ReaxFF potential to well established results associated with the Si–H bending mode is given in Appendix D.

4.1.4 Dipole Scattering

Next, one can consider the rate expressions for inelastic excitation channels which are usually induced due to external perturbations. As already mentioned in the experimental Section 2.1.1, one inelastic process is related to the electric field caused by moving electrons. The resulting field can interact with the dipole created by the vibration of the adsorbate complex and induce transitions between vibrational eigenstates. The dipole scattering rate can be described due to an electrostatic coupling to the transition dipole moment $\mathit{\mu }$ and can be calculated by a formula derived by Persson et al. [277–279]

$$\begin{array}{}\text{(4.22)}& {\mathrm{\Gamma }}_{i,f}^{\mathrm{dip}}=\frac{I}{e}{\left|\frac{⟨{\varphi }_f|\mathit{\mu }|{\varphi }_i⟩}{e{a}_0}\right|}^2,\end{array}$$

where $e$ is the elementary charge, $a_0$ the Bohr radius, $\hat{\mu }$ the dipole moment and $I$ the current. Using the dipole moment vector extracted in Sec. 3.3, in particular the $x$ and $y$ component due to the current density in the channel of a MOSFET being parallel to the (100) surface, one can estimate the importance of this component. The dipole–induced transition rates ${\mathrm{\Gamma }}_{i,f}^{\mathrm{dip}}$ for $i\Rightarrow i+1$ and $i\Rightarrow i+2$ transitions have been calculated for different currents and compared to the vibrational upward rates ${\mathrm{\Gamma }}_{i,f}^{\mathrm{vib}}$⁸ at $T=300K$, see Fig. 4.4. Current densities in scaled MOSFETs already exceed $100nA/n{m}^2$, which is of similar order as tunneling currents in Si–H related STM experiments ($1-10nA$), assuming that the STM electron beam is of atomic dimensions. One can see that only at higher currents the respective upward rates for $\mathrm{\Delta }n=+1$ are mainly determined by dipole induced scattering, whereas transitions to higher excited states ($\mathrm{\Delta }n=+2$) are already orders of magnitude lower. Although Fig. 4.4 actually shows a non–negligible contribution of ${\mathrm{\Gamma }}_{i,f}^{\mathrm{dip}}$, the importance of dipole–induced excitations is limited in hot–carrier related device reliability issues. Only at the source side, where carriers are still close to equilibrium, but a high current density is present, they potentially play a (minor) role. Nevertheless, starting from the channel region, carriers have already gained enough energy to (potentially) scatter into an available resonance which clearly dominates the upward transitions.

4.1.5 Resonance Scattering

As already mentioned in Section 2.1.1 another possible excitation channel is related to adsorbate resonances. A resonance scattering model accounts for carriers tunneling into available molecular orbitals (MOs) of the adsorbate. In the case of the Si–H bond the resonance states are associated with the LUMO (HOMO), forming a temporal negative (positive) ion resonance. The modified internuclear potential thereby induces a nuclear relaxation of the system. When the electron returns to the substrate an inelastic relaxation process transfers energy to the system, leaving the neutral ground state PEC $V(q)$ in a vibrationally excited state. Two formulations of a vibrational heating model for this mechanism have been developed: an incoherent multiple single-step process [280] and a theory of coherent multiple vibrational excitations [281]. Both variants are schematically illustrated in Fig. 4.5.

Whereas the incoherent model requires a rather high current density (scattering electrons) compared to the vibrational lifetime, the coherent formulation accounts for overtone transitions which results in an excitation path with a smaller number of intermediate states [282, 283]. Recent results [119] have shown that indeed fewer electrons are needed to dissociate H than expected from the incoherent formulation and support the coherent formulation proposed by Persson et al. [284] and Salam et al. [281]. The utilized expression within this work for the transition rates is therefore based on the above descriptions and was already applied to similar problems [285] as discussed here. In this formulation it is assumed that an electron with incident energy $ϵ$ can induce an excitation from state $|{\varphi }_i⟩$ to $|{\varphi }_f⟩$ via the vibrational eigenstate $|{\psi }_j⟩$ of the resonance with

$$\begin{array}{}\text{(4.23)}& {\mathrm{\Gamma }}_{i,f}^{\mathrm{res}}=\frac{4{\mathrm{\Delta }}_{\mathrm{res}}^2}{\pi \hslash }\int dϵf(ϵ)[1-f(ϵ)]\cdot |\sum _j\frac{⟨{\varphi }_f|{\psi }_j⟩⟨{\psi }_j|{\varphi }_i⟩}{{ϵ}_{\mathrm{res}}-ϵ+{ϵ}_i-{ϵ}_j+i{\mathrm{\Delta }}_{\mathrm{res}}}{|}^2.\end{array}$$

This expression includes the energetic position of the resonance with respect to the Fermi level, ${ϵ}_{\mathrm{res}}$, its width as the inverse of the resonance lifetime, ${\mathrm{\Delta }}_{\mathrm{res}}=\hslash /{\tau }_{\mathrm{res}}$, as well as the energies of the eigenstates $|{\varphi }_i⟩$ and $|{\psi }_j⟩$. Furthermore, it accounts for the energy distribution function $f(ϵ)$ of charge carriers at the interface, which is a crucial ingredient as pointed out in the Outline, see the Chapters 1 and 2, of this thesis. All parameters entering the above formula can be extracted (or at least reasonably approximated) from DFT results, see Chap. 3.

4.1.6 Model Formulation

Based on the theoretical treatment introduced above, a model for describing hot-carrier induced damage at the Si/SiO2 interface can be developed. The approach includes all the aforementioned processes, namely the vibrational relaxation mechanism, (4.13) and (4.21), as well as the inelastic channels, dipole, (4.22), and resonance scattering, (4.23). The total transition rate from the vibrational eigenstate $|{\varphi }_i⟩$ to $|{\varphi }_f⟩$ of the Si–H system is given by (4.7).

Experimental quantification of HCD is done on a macroscopic timescale. The device is stressed above operating conditions between seconds and hours to extract a gradual change of the device characteristic to ultimately evaluate its lifetime. On the other hand, the Si–H system interacting with the environment will tend towards a steady state on a much faster timescale, within femto– to nanoseconds. Thus, the time scales can be separated and only a quasi-equilibrium solution of (4.3) is required

$$\begin{array}{}\text{(4.24)}& \frac{\partial {\rho }_{\mathrm{S}}}{\partial t}={\mathcal{L}}_{\mathrm{S}}{\rho }_{\mathrm{S}}+{\mathcal{L}}_{\mathrm{D}}{\rho }_{\mathrm{S}}=0,\end{array}$$

to calculate the population $P_i$ of each individual vibrational state of the system. Once the new equilibrium of the ground state potential $V(q)$ is known, one can calculate the bond breaking rate, defined as the transition from the left to the right well in the ground state potential, see Fig. 4.1, as tunneling rates through a potential barrier, employing the WKB–aproximation:

$$\begin{array}{}\text{(4.25)}& {\mathrm{\Gamma }}^{\mathrm{break}}=\sum _i{\mathrm{\Gamma }}_{i,\mathrm{WKB}}{P}_i=\sum _i\exp (-\frac{2}{\hslash }\underset{x_{1,i}}{\overset{x_{2,i}}{\int }}\sqrt{2(V-E_i)}\mathrm{d}x)P_i,\end{array}$$

with $P_i$ being the population of the eigenstate $i$ and ${\mathrm{\Gamma }}_{i,\mathrm{WKB}}$ the tunneling rate which takes the classical turning points $x_1$ and $x_2$ for the respective energy levels $E_i$ into account. Within this approach only eigenstates localized in the left well are considered by implicitly projecting out all $|{\varphi }_{\mathrm{R}}⟩$ in the above calculation. A more detailed description of how to calculate the final rate ${\mathrm{\Gamma }}^{\mathrm{break}}$ is given in Appendix E.

Knowing the bond–breaking rate ${\mathrm{\Gamma }}^{\mathrm{break}}$ allows to formulate a reaction equation for the evolution of the concentration of broken and intact Si–H bonds with stress time. Two states for the forward reaction can be assumed, see Secs. 3.2 and [MJJ3], which is in analogy to previous observations [35, 201, 244], an intact Si–H bond (left well) and a broken state (right well) described as a $P_{\mathrm{b}}$ center and the released H atom. Once the hydrogen is in the next but one bond–center site, it is likely that it is mobile along the interface, see Sec. 3.6 and  [251, 253, 254], facing a rather small barrier for hopping laterally within the subinterfacial Si region. Thus, once the H is released, it potentially triggers an additional reaction, e.g. the formation of defects or H${}_2$. The backward reaction, on the other hand, passivating the defect with one hydrogen, is assumed to proceed via a different reaction pathway, namely by breaking an H${}_2$ modelcule, see 3.6 and [201, 221, 248].

The reaction equation including the aforementioned chemical reactions is therefore

$$\begin{array}{}\text{(4.26)}& \frac{\mathrm{d}\left[\mathrm{SiH}\right]}{\mathrm{d}t}=-{\mathrm{\Gamma }}^{\mathrm{break}}\left[\mathrm{SiH}\right]+{\mathrm{\Gamma }}^{\mathrm{pass}}\left[{\mathrm{H}}_2\right]\left[P_{\mathrm{b}}\right],\end{array}$$

where ${\mathrm{\Gamma }}^{\mathrm{pass}}$ describes cracking of H${}_2$ and the subsequent passivation of $P_{\mathrm{b}}$ centers. According to  [221, 222, 244, 248] this process is purely thermally activated and the rate constant is given by an Arrhenius equation ${\mathrm{\Gamma }}^{\mathrm{pass}}={\mathrm{\Gamma }}_0^{\mathrm{pass}}\exp (-E_{\mathrm{pass}}/k_{\mathrm{B}}T)$ superimposed with a Gaussian distribution of the passivation energy $E_{\mathrm{pass}}$⁹. The inferred parameters are $E_{\mathrm{p}}\sim 1.51eV$, ${\sigma }_{E_{\mathrm{p}}}\sim 0.12eV$ and ${\mathrm{\Gamma }}_0^{\mathrm{pass}}\sim 5\times {10}^{-6}c{m}^3/\mathrm{s}$¹⁰, which were also used in a recent study to model the recovery of HCD induced damage in an actual MOSFET [248, 249]. In addition, one equation describing the maximum number of pristine Si–H bonds is introduced as

$$\begin{gathered} \begin{array}{}\text{(4.27)}& \left[\mathrm{SiH}\right]+\left[\mathrm{H}\right]={\left[\mathrm{H}\right]}_{\mathrm{tot}}, \\ \end{array} \end{gathered}$$

while the value for [H${}_2$] was set to a constant concentration of $5\times {10}^{17}c{m}^{-3}$, given by the physical solubility of H${}_2$ in vitreous silica [286]. Introducing the quantity ${\mathrm{N}}_0\widehat{=}{\left[\mathrm{SiH}\right]}_{\max }$ and solving for $f_{\mathrm{Pb}}\widehat{=}\left[{\mathrm{P}}_{\mathrm{b}}\right]/{\mathrm{N}}_0=1-\left[\mathrm{SiH}\right]/{\mathrm{N}}_0$, the probability that a bond is broken is then given by

$$\begin{array}{}\text{(4.28)}& \frac{\partial f_{\mathrm{Pb}}}{\partial t}=(1-f_{\mathrm{Pb}}){\mathrm{\Gamma }}^{\mathrm{desorb}}-f_{\mathrm{Pb}}{\mathrm{\Gamma }}^{\mathrm{pass}}[{\mathrm{H}}_2].\end{array}$$

In order to represent the density of bonds at the Si/SiO${}_2$ interface, it is necessary to randomly sample the normally distributed quantities of the Si–H bond and scale them according to the technology dependent pristine density ${\mathrm{N}}_0$, which is, virtually, the only fitting parameter in our model.

4.1.7 Relation & Connection to Previous Models

The framework introduced above is computationally challenging compared to previously used TCAD models describing degradation phenomena, in particular due to the need to randomly sample the rather large parameter space together with the numerical solution for the energy profiles. In order to decrease its complexity and simultaneously increase the usability of the model within TCAD tools, further approximations can be applied. The ground and resonance potential profiles $V(q)$ and $V^{-}(q)/V^{+}(q)$ can be approximated using harmonic potentials, see the right panel of Fig. 4.6. Within this simplification the proposed resonance scattering formalism, see Sec. 4.1.5, is still valid. A charged carrier with an incident energy $ϵ$ triggers a vibrational excitation from state $|{\varphi }_i⟩$ to $|{\varphi }_f⟩$ in the ground state $V(q)$ via the eigenstate $|{\psi }_j⟩$ of the resonance. Comparing the resonance based description using harmonic potentials to previous modeling attempts by Hess [47, 49, 50, 52], Bravaix [34, 53, 56] and Tyaginov [MJJ5, 35, MJJ12, MJC13, MJC16, MJC4] reveals an unambiguous connection. Obviously, the coherent formulation of the inelastic resonance scattering mechanism, see Sec. 4.1.5, naturally maps the single particle (SP) and multiple particle (MP) mechanisms discussed in Chapter 2 onto a single fundamental excitation channel. Furthermore, the analogy allows to reinterpret previous formulations in a physically meaningful context. For instance the model developed by Tyaginov et al. uses the so–called Acceleration Integral to calculate the corresponding SP and MP transition rates

$$\begin{array}{}\text{(4.29)}& {\mathrm{\Gamma }}_{\mathrm{MP}/\mathrm{SP}}={\int }_{{ϵ}_{\mathrm{th}}}f(ϵ)g(ϵ)v(ϵ){\sigma }_{\mathrm{MP}/\mathrm{SP}}(ϵ)\mathrm{d}ϵ,\end{array}$$

where $f(ϵ)$ is the EDF, $g(ϵ)$ the DOS, $v(ϵ)$ corresponds to the carrier velocity and $\sigma (ϵ)$ is a phenomenologically derived capture cross section. The threshold energy, i.e. the minimum carrier energy needed to trigger the respective mechanism, is denoted as ${ϵ}_{\mathrm{th}}$. The energy dependent cross section is assumed to follow an empirical Keldysh–like expression $\sigma (E)={\sigma }_0(ϵ-{ϵ}_{\mathrm{th}}{)}^p$ with $p_{\mathrm{MP}}=1$ and $p_{\mathrm{SP}}=11$. Various exponents have been extracted in the literature, depending on the actual implementation of the utilized model, and have been initially used to capture the multiple dependencies of HCD [33–35, MJJ12, 53, 57].

Using the quantum–chemistry formulation described within this work allows to reformulate (4.29) in terms of a single physical mechanisms

$$\begin{array}{}\text{(4.30)}& {\mathrm{\Gamma }}_{i,f}=\frac{4{\mathrm{\Delta }}_{\mathrm{res}}^2}{\pi }{\int }_{ϵ}f(ϵ)g(ϵ)v(ϵ){\sigma }_{\mathrm{eff}}(ϵ)\mathrm{d}ϵ,\end{array}$$

with the effective cross section

$$\begin{array}{}\text{(4.31)}& {\sigma }_{\mathrm{eff}}(ϵ)={\sigma }_0|\sum _j\frac{⟨{\varphi }_f|{\psi }_j⟩⟨{\psi }_j|{\varphi }_i⟩}{{ϵ}_{\mathrm{res}}-ϵ+{ϵ}_i-{ϵ}_j+i{\mathrm{\Delta }}_{\mathrm{res}}}{|}^2.\end{array}$$

In analogy to (4.23) it is the probability of the vibrational transition from $|{\varphi }_i⟩$ to $|{\varphi }_f⟩$ in the ground state potential via the resonance eigenstate $|{\psi }_j⟩$, given by the overlap of the wavefunctions, weighted by the resonance position represented by a Lorentzian. Thus, this formulation simultaneously accounts for neighbouring transitions ${\mathrm{\Delta }}_{n,n+1}$ (ladder climbing, MP mechanism), overtone transitions ${\mathrm{\Delta }}_{n,m}$ as well as transitions directly causing dissociation ${\mathrm{\Delta }}_{n,f+1}$ where $f$ denotes the last bound state in the potential (SP process). Note that such a resonance based description was already proposed by the group of Hess [52], albeit not rigorously included in their calculations. A detailed comparison of the cross section’s energy dependence is given in Fig. 4.7. While the empirical Keldysh–like expression does not depend on the current eigenstate, i.e. the level of bond excitation, such subtleties are included in (4.31) via ${ϵ}_{\mathrm{res}}+{ϵ}_i-{ϵ}_j$. Note that within the proposed effective cross section overtone transitions can potentially have a higher excitation probability due to the wavefunction overlaps.

Additionally, due to the small effective dipole moment ${\mathit{\mu }}_{\mathrm{eff}}$ of the Si–H bond, see Sec. 3.3, and the the marginal relevance of dipole scattering rates, resulting therefrom Sec. 4.1.4, this component can be neglected. On the other hand, the harmonic scaling law for the vibrational lifetime (4.21), see Sec. 4.1.3, is now strictly fulfilled and can be applied. Therefore, a total transition rate can again be calculated, using the individual contributions of the resonance scattering mechanism (4.30) and (4.31) and the vibrational dissipation (4.21) and (4.13)

$$\begin{array}{}\text{(4.32)}& {\mathrm{\Gamma }}_{i,f}^{\mathrm{tot}}={\mathrm{\Gamma }}_{i,f}^{\mathrm{res}}+{\mathrm{\Gamma }}_{i,f}^{\mathrm{vib}}.\end{array}$$

Subsequently, again the new quasi–equilibrium solution of the Pauli Master equation, see (4.24), can be calculated which yields the individual state populations $P_i$ of the ground state harmonic oscillator. Including the first continuum state (the first state above the Si–H bond breakage energy), the final bond breakage rate ${\mathrm{\Gamma }}^{\mathrm{break}}$ can be classically defined as the vibrational eigenstate population $P_i$ times the total transition rate from $i$ to the continuum state $f+1$. The evolution of broken and intact Si–H bonds with stress time is again given by the reaction equation (4.26) and (4.28), respectively.

The major advantage of harmonic approximations is the availability of analytic solutions for the eigenvalues and eigenvectors. Furthermore, due to the simpler potential energy profile it effectively reduces the parameter set and allows to efficiently sample the distribution on a 3D grid.