\(\newcommand{\footnotename}{footnote}\)
\(\def \LWRfootnote {1}\)
\(\newcommand {\footnote }[2][\LWRfootnote ]{{}^{\mathrm {#1}}}\)
\(\newcommand {\footnotemark }[1][\LWRfootnote ]{{}^{\mathrm {#1}}}\)
\(\let \LWRorighspace \hspace \)
\(\renewcommand {\hspace }{\ifstar \LWRorighspace \LWRorighspace }\)
\(\newcommand {\mathnormal }[1]{{#1}}\)
\(\newcommand \ensuremath [1]{#1}\)
\(\newcommand {\LWRframebox }[2][]{\fbox {#2}} \newcommand {\framebox }[1][]{\LWRframebox } \)
\(\newcommand {\setlength }[2]{}\)
\(\newcommand {\addtolength }[2]{}\)
\(\newcommand {\setcounter }[2]{}\)
\(\newcommand {\addtocounter }[2]{}\)
\(\newcommand {\arabic }[1]{}\)
\(\newcommand {\number }[1]{}\)
\(\newcommand {\noalign }[1]{\text {#1}\notag \\}\)
\(\newcommand {\cline }[1]{}\)
\(\newcommand {\directlua }[1]{\text {(directlua)}}\)
\(\newcommand {\luatexdirectlua }[1]{\text {(directlua)}}\)
\(\newcommand {\protect }{}\)
\(\def \LWRabsorbnumber #1 {}\)
\(\def \LWRabsorbquotenumber "#1 {}\)
\(\newcommand {\LWRabsorboption }[1][]{}\)
\(\newcommand {\LWRabsorbtwooptions }[1][]{\LWRabsorboption }\)
\(\def \mathchar {\ifnextchar "\LWRabsorbquotenumber \LWRabsorbnumber }\)
\(\def \mathcode #1={\mathchar }\)
\(\let \delcode \mathcode \)
\(\let \delimiter \mathchar \)
\(\let \LWRref \ref \)
\(\renewcommand {\ref }{\ifstar \LWRref \LWRref }\)
\(\newcommand {\toprule }[1][]{\hline }\)
\(\let \midrule \toprule \)
\(\let \bottomrule \toprule \)
\(\def \LWRbooktabscmidruleparen (#1)#2{}\)
\(\newcommand {\LWRbooktabscmidrulenoparen }[1]{}\)
\(\newcommand {\cmidrule }[1][]{\ifnextchar (\LWRbooktabscmidruleparen \LWRbooktabscmidrulenoparen }\)
\(\newcommand {\morecmidrules }{}\)
\(\newcommand {\specialrule }[3]{\hline }\)
\(\newcommand {\addlinespace }[1][]{}\)
\(\def \LWRpagenote {1}\)
\(\newcommand {\pagenote }[2][\LWRpagenote ]{{}^{\mathrm {#1}}}\)
\(\newcommand {\Vg }{V_\mathrm {G}}\)
\(\newcommand {\Vd }{V_\mathrm {D}}\)
\(\newcommand {\Vdd }{V_\mathrm {DD}}\)
\(\newcommand {\interface }{Si/SiO$_2$}\)
\(\newcommand {\sio }{SiO$_2$}\)
\(\newcommand {\VgVd }{\{$V_\mathrm {G},V_\mathrm {D}$\}}\)
\(\newcommand {\taue }{$\tau _\mathrm {e}$}\)
\(\newcommand {\tauc }{$\tau _\mathrm {c}$}\)
\(\newcommand {\eox }{$F_\mathrm {ox}$}\)
\(\newcommand {\idx }{_{i,j}}\)
\(\newcommand {\idxrm }{_\mathrm {i,j}}\)
\(\newcommand {\kb }{k_\mathrm {B}}\)
\(\newcommand {\nmpeq }{NMP$_\mathrm {eq.}$}\)
\(\newcommand {\nmpeqext }{NMP$_\mathrm {eq.+II}$}\)
\(\newcommand {\nmpneq }{NMP$_\mathrm {neq.}$}\)
\(\newcommand {\units }[2]{$\SI {#1}{#2}$}\)
\(\newcommand {\fptOne }{$\tau _\mathrm {e,FPT}^{1^\prime ,1}$}\)
\(\newcommand {\fptTwo }{$\tau _\mathrm {e,FPT}^{2^\prime ,1}$}\)
\(\newcommand {\ts }{$t_\mathrm {s}$}\)
\(\newcommand {\dvth }{\Delta V_\mathrm {th}}\)
\(\newcommand {\Pb }{$P_\mathrm {b}$}\)
\(\newcommand {\Pbzero }{$P_\mathrm {b0}$}\)
\(\newcommand {\Pbone }{$P_\mathrm {b1}$}\)
\(\newcommand {\Nit }{$N_\mathrm {it}(x)$}\)
\( \newcommand {\multicolumn }[3]{#3}\)
\(\newcommand {\tothe }[1]{^{#1}}\)
\(\newcommand {\raiseto }[2]{{#2}^{#1}}\)
\(\newcommand {\LWRsiunitxEND }{}\)
\(\def \LWRsiunitxang #1;#2;#3;#4\LWRsiunitxEND {\ifblank {#1}{}{\num {#1}\degree }\ifblank {#2}{}{\num {#2}^{\unicode {x2032}}}\ifblank {#3}{}{\num {#3}^{\unicode {x2033}}}}\)
\(\newcommand {\ang }[2][]{\LWRsiunitxang #2;;;\LWRsiunitxEND }\)
\(\newcommand {\LWRsiunitxnumscientific }[2]{\ifblank {#1}{}{\ifstrequal {#1}{-}{-}{\LWRsiunitxprintdecimal {#1}\times }}10^{\LWRsiunitxprintdecimal {#2}} }\)
\(\def \LWRsiunitxnumplus #1+#2+#3\LWRsiunitxEND {\ifblank {#2} {\LWRsiunitxprintdecimal {#1}}{\ifblank {#1}{\LWRsiunitxprintdecimal {#2}}{\LWRsiunitxprintdecimal {#1}\unicode
{x02B}\LWRsiunitxprintdecimal {#2}}}}\)
\(\def \LWRsiunitxnumminus #1-#2-#3\LWRsiunitxEND {\ifblank {#2} {\LWRsiunitxnumplus #1+++\LWRsiunitxEND }{\LWRsiunitxprintdecimal {#1}\unicode {x02212}\LWRsiunitxprintdecimal {#2}}}\)
\(\def \LWRsiunitxnumpm #1+-#2+-#3\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumminus #1---\LWRsiunitxEND }{\LWRsiunitxprintdecimal {#1}\unicode {x0B1}\LWRsiunitxprintdecimal {#2}}}\)
\(\def \LWRsiunitxnumx #1x#2x#3x#4\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumpm #1+-+-\LWRsiunitxEND }{\ifblank {#3}{\LWRsiunitxprintdecimal {#1}\times \LWRsiunitxprintdecimal
{#2}}{\LWRsiunitxprintdecimal {#1}\times \LWRsiunitxprintdecimal {#2}\times \LWRsiunitxprintdecimal {#3}}}}\)
\(\def \LWRsiunitxnumD #1D#2D#3\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumx #1xxxxx\LWRsiunitxEND }{\mathrm {\LWRsiunitxnumscientific {#1}{#2}}}}\)
\(\def \LWRsiunitxnumd #1d#2d#3\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumD #1DDD\LWRsiunitxEND }{\mathrm {\LWRsiunitxnumscientific {#1}{#2}}}}\)
\(\def \LWRsiunitxnumE #1E#2E#3\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumd #1ddd\LWRsiunitxEND }{\mathrm {\LWRsiunitxnumscientific {#1}{#2}}}}\)
\(\def \LWRsiunitxnume #1e#2e#3\LWRsiunitxEND {\ifblank {#2}{\LWRsiunitxnumE #1EEE\LWRsiunitxEND }{\mathrm {\LWRsiunitxnumscientific {#1}{#2}}}}\)
\(\def \LWRsiunitxnumcomma #1,#2,#3\LWRsiunitxEND {\ifblank {#2} {\LWRsiunitxnume #1eee\LWRsiunitxEND } {\LWRsiunitxnume #1.#2eee\LWRsiunitxEND } }\)
\(\newcommand {\num }[2][]{\LWRsiunitxnumcomma #2,,,\LWRsiunitxEND }\)
\(\newcommand {\si }[2][]{\mathrm {#2}}\)
\(\def \LWRsiunitxSIopt #1[#2]#3{{#2}\num {#1}{#3}}\)
\(\newcommand {\LWRsiunitxSI }[2]{\num {#1}\,{#2}}\)
\(\newcommand {\SI }[2][]{\ifnextchar [{\LWRsiunitxSIopt {#2}}{\LWRsiunitxSI {#2}}}\)
\(\newcommand {\numlist }[2][]{\mathrm {#2}}\)
\(\newcommand {\numrange }[3][]{\num {#2}\,\unicode {x2013}\,\num {#3}}\)
\(\newcommand {\SIlist }[3][]{\mathrm {#2\,#3}}\)
\(\newcommand {\SIrange }[4][]{\num {#2}\,#4\,\unicode {x2013}\,\num {#3}\,#4}\)
\(\newcommand {\tablenum }[2][]{\mathrm {#2}}\)
\(\newcommand {\ampere }{\mathrm {A}}\)
\(\newcommand {\candela }{\mathrm {cd}}\)
\(\newcommand {\kelvin }{\mathrm {K}}\)
\(\newcommand {\kilogram }{\mathrm {kg}}\)
\(\newcommand {\metre }{\mathrm {m}}\)
\(\newcommand {\mole }{\mathrm {mol}}\)
\(\newcommand {\second }{\mathrm {s}}\)
\(\newcommand {\becquerel }{\mathrm {Bq}}\)
\(\newcommand {\degreeCelsius }{\unicode {x2103}}\)
\(\newcommand {\coulomb }{\mathrm {C}}\)
\(\newcommand {\farad }{\mathrm {F}}\)
\(\newcommand {\gray }{\mathrm {Gy}}\)
\(\newcommand {\hertz }{\mathrm {Hz}}\)
\(\newcommand {\henry }{\mathrm {H}}\)
\(\newcommand {\joule }{\mathrm {J}}\)
\(\newcommand {\katal }{\mathrm {kat}}\)
\(\newcommand {\lumen }{\mathrm {lm}}\)
\(\newcommand {\lux }{\mathrm {lx}}\)
\(\newcommand {\newton }{\mathrm {N}}\)
\(\newcommand {\ohm }{\mathrm {\Omega }}\)
\(\newcommand {\pascal }{\mathrm {Pa}}\)
\(\newcommand {\radian }{\mathrm {rad}}\)
\(\newcommand {\siemens }{\mathrm {S}}\)
\(\newcommand {\sievert }{\mathrm {Sv}}\)
\(\newcommand {\steradian }{\mathrm {sr}}\)
\(\newcommand {\tesla }{\mathrm {T}}\)
\(\newcommand {\volt }{\mathrm {V}}\)
\(\newcommand {\watt }{\mathrm {W}}\)
\(\newcommand {\weber }{\mathrm {Wb}}\)
\(\newcommand {\day }{\mathrm {d}}\)
\(\newcommand {\degree }{\mathrm {^\circ }}\)
\(\newcommand {\hectare }{\mathrm {ha}}\)
\(\newcommand {\hour }{\mathrm {h}}\)
\(\newcommand {\litre }{\mathrm {l}}\)
\(\newcommand {\liter }{\mathrm {L}}\)
\(\newcommand {\arcminute }{^\prime }\)
\(\newcommand {\minute }{\mathrm {min}}\)
\(\newcommand {\arcsecond }{^{\prime \prime }}\)
\(\newcommand {\tonne }{\mathrm {t}}\)
\(\newcommand {\astronomicalunit }{au}\)
\(\newcommand {\atomicmassunit }{u}\)
\(\newcommand {\bohr }{\mathit {a}_0}\)
\(\newcommand {\clight }{\mathit {c}_0}\)
\(\newcommand {\dalton }{\mathrm {D}_\mathrm {a}}\)
\(\newcommand {\electronmass }{\mathit {m}_{\mathrm {e}}}\)
\(\newcommand {\electronvolt }{\mathrm {eV}}\)
\(\newcommand {\elementarycharge }{\mathit {e}}\)
\(\newcommand {\hartree }{\mathit {E}_{\mathrm {h}}}\)
\(\newcommand {\planckbar }{\mathit {\unicode {x210F}}}\)
\(\newcommand {\angstrom }{\mathrm {\unicode {x212B}}}\)
\(\let \LWRorigbar \bar \)
\(\newcommand {\bar }{\mathrm {bar}}\)
\(\newcommand {\barn }{\mathrm {b}}\)
\(\newcommand {\bel }{\mathrm {B}}\)
\(\newcommand {\decibel }{\mathrm {dB}}\)
\(\newcommand {\knot }{\mathrm {kn}}\)
\(\newcommand {\mmHg }{\mathrm {mmHg}}\)
\(\newcommand {\nauticalmile }{\mathrm {M}}\)
\(\newcommand {\neper }{\mathrm {Np}}\)
\(\newcommand {\yocto }{\mathrm {y}}\)
\(\newcommand {\zepto }{\mathrm {z}}\)
\(\newcommand {\atto }{\mathrm {a}}\)
\(\newcommand {\femto }{\mathrm {f}}\)
\(\newcommand {\pico }{\mathrm {p}}\)
\(\newcommand {\nano }{\mathrm {n}}\)
\(\newcommand {\micro }{\mathrm {\unicode {x00B5}}}\)
\(\newcommand {\milli }{\mathrm {m}}\)
\(\newcommand {\centi }{\mathrm {c}}\)
\(\newcommand {\deci }{\mathrm {d}}\)
\(\newcommand {\deca }{\mathrm {da}}\)
\(\newcommand {\hecto }{\mathrm {h}}\)
\(\newcommand {\kilo }{\mathrm {k}}\)
\(\newcommand {\mega }{\mathrm {M}}\)
\(\newcommand {\giga }{\mathrm {G}}\)
\(\newcommand {\tera }{\mathrm {T}}\)
\(\newcommand {\peta }{\mathrm {P}}\)
\(\newcommand {\exa }{\mathrm {E}}\)
\(\newcommand {\zetta }{\mathrm {Z}}\)
\(\newcommand {\yotta }{\mathrm {Y}}\)
\(\newcommand {\percent }{\mathrm {\%}}\)
\(\newcommand {\meter }{\mathrm {m}}\)
\(\newcommand {\metre }{\mathrm {m}}\)
\(\newcommand {\gram }{\mathrm {g}}\)
\(\newcommand {\kg }{\kilo \gram }\)
\(\newcommand {\of }[1]{_{\mathrm {#1}}}\)
\(\newcommand {\squared }{^2}\)
\(\newcommand {\square }[1]{\mathrm {#1}^2}\)
\(\newcommand {\cubed }{^3}\)
\(\newcommand {\cubic }[1]{\mathrm {#1}^3}\)
\(\newcommand {\per }{/}\)
\(\newcommand {\celsius }{\unicode {x2103}}\)
\(\newcommand {\fg }{\femto \gram }\)
\(\newcommand {\pg }{\pico \gram }\)
\(\newcommand {\ng }{\nano \gram }\)
\(\newcommand {\ug }{\micro \gram }\)
\(\newcommand {\mg }{\milli \gram }\)
\(\newcommand {\g }{\gram }\)
\(\newcommand {\kg }{\kilo \gram }\)
\(\newcommand {\amu }{\mathrm {u}}\)
\(\newcommand {\pm }{\pico \metre }\)
\(\newcommand {\nm }{\nano \metre }\)
\(\newcommand {\um }{\micro \metre }\)
\(\newcommand {\mm }{\milli \metre }\)
\(\newcommand {\cm }{\centi \metre }\)
\(\newcommand {\dm }{\deci \metre }\)
\(\newcommand {\m }{\metre }\)
\(\newcommand {\km }{\kilo \metre }\)
\(\newcommand {\as }{\atto \second }\)
\(\newcommand {\fs }{\femto \second }\)
\(\newcommand {\ps }{\pico \second }\)
\(\newcommand {\ns }{\nano \second }\)
\(\newcommand {\us }{\micro \second }\)
\(\newcommand {\ms }{\milli \second }\)
\(\newcommand {\s }{\second }\)
\(\newcommand {\fmol }{\femto \mol }\)
\(\newcommand {\pmol }{\pico \mol }\)
\(\newcommand {\nmol }{\nano \mol }\)
\(\newcommand {\umol }{\micro \mol }\)
\(\newcommand {\mmol }{\milli \mol }\)
\(\newcommand {\mol }{\mol }\)
\(\newcommand {\kmol }{\kilo \mol }\)
\(\newcommand {\pA }{\pico \ampere }\)
\(\newcommand {\nA }{\nano \ampere }\)
\(\newcommand {\uA }{\micro \ampere }\)
\(\newcommand {\mA }{\milli \ampere }\)
\(\newcommand {\A }{\ampere }\)
\(\newcommand {\kA }{\kilo \ampere }\)
\(\newcommand {\ul }{\micro \litre }\)
\(\newcommand {\ml }{\milli \litre }\)
\(\newcommand {\l }{\litre }\)
\(\newcommand {\hl }{\hecto \litre }\)
\(\newcommand {\uL }{\micro \liter }\)
\(\newcommand {\mL }{\milli \liter }\)
\(\newcommand {\L }{\liter }\)
\(\newcommand {\hL }{\hecto \liter }\)
\(\newcommand {\mHz }{\milli \hertz }\)
\(\newcommand {\Hz }{\hertz }\)
\(\newcommand {\kHz }{\kilo \hertz }\)
\(\newcommand {\MHz }{\mega \hertz }\)
\(\newcommand {\GHz }{\giga \hertz }\)
\(\newcommand {\THz }{\tera \hertz }\)
\(\newcommand {\mN }{\milli \newton }\)
\(\newcommand {\N }{\newton }\)
\(\newcommand {\kN }{\kilo \newton }\)
\(\newcommand {\MN }{\mega \newton }\)
\(\newcommand {\Pa }{\pascal }\)
\(\newcommand {\kPa }{\kilo \pascal }\)
\(\newcommand {\MPa }{\mega \pascal }\)
\(\newcommand {\GPa }{\giga \pascal }\)
\(\newcommand {\mohm }{\milli \ohm }\)
\(\newcommand {\kohm }{\kilo \ohm }\)
\(\newcommand {\Mohm }{\mega \ohm }\)
\(\newcommand {\pV }{\pico \volt }\)
\(\newcommand {\nV }{\nano \volt }\)
\(\newcommand {\uV }{\micro \volt }\)
\(\newcommand {\mV }{\milli \volt }\)
\(\newcommand {\V }{\volt }\)
\(\newcommand {\kV }{\kilo \volt }\)
\(\newcommand {\W }{\watt }\)
\(\newcommand {\uW }{\micro \watt }\)
\(\newcommand {\mW }{\milli \watt }\)
\(\newcommand {\kW }{\kilo \watt }\)
\(\newcommand {\MW }{\mega \watt }\)
\(\newcommand {\GW }{\giga \watt }\)
\(\newcommand {\J }{\joule }\)
\(\newcommand {\uJ }{\micro \joule }\)
\(\newcommand {\mJ }{\milli \joule }\)
\(\newcommand {\kJ }{\kilo \joule }\)
\(\newcommand {\eV }{\electronvolt }\)
\(\newcommand {\meV }{\milli \electronvolt }\)
\(\newcommand {\keV }{\kilo \electronvolt }\)
\(\newcommand {\MeV }{\mega \electronvolt }\)
\(\newcommand {\GeV }{\giga \electronvolt }\)
\(\newcommand {\TeV }{\tera \electronvolt }\)
\(\newcommand {\kWh }{\kilo \watt \hour }\)
\(\newcommand {\F }{\farad }\)
\(\newcommand {\fF }{\femto \farad }\)
\(\newcommand {\pF }{\pico \farad }\)
\(\newcommand {\K }{\mathrm {K}}\)
\(\newcommand {\dB }{\mathrm {dB}}\)
\(\newcommand {\kibi }{\mathrm {Ki}}\)
\(\newcommand {\mebi }{\mathrm {Mi}}\)
\(\newcommand {\gibi }{\mathrm {Gi}}\)
\(\newcommand {\tebi }{\mathrm {Ti}}\)
\(\newcommand {\pebi }{\mathrm {Pi}}\)
\(\newcommand {\exbi }{\mathrm {Ei}}\)
\(\newcommand {\zebi }{\mathrm {Zi}}\)
\(\newcommand {\yobi }{\mathrm {Yi}}\)
\(\newcommand {\intertext }[1]{\text {#1}\notag \\}\)
\(\let \Hat \hat \)
\(\let \Check \check \)
\(\let \Tilde \tilde \)
\(\let \Acute \acute \)
\(\let \Grave \grave \)
\(\let \Dot \dot \)
\(\let \Ddot \ddot \)
\(\let \Breve \breve \)
\(\let \Bar \bar \)
\(\let \Vec \vec \)
\(\require {physics}\)
\(\newcommand {\LWRsfrac }[2][/]{{}^\LWRsfracnumerator \!#1{}_{#2}}\)
\(\newcommand {\sfrac }[2][]{\def \LWRsfracnumerator {#2}\LWRsfrac }\)
\(\newcommand {\nicefrac }[3][]{\mathinner {{}^{#2}\!/\!_{#3}}}\)
\(\newcommand {\bm }[1]{\boldsymbol {#1}}\)
\(\def \LWRsiunitxprintdecimalsub #1,#2,#3\LWRsiunitxEND {\mathrm {#1}\ifblank {#2}{}{.\mathrm {#2}}}\)
\(\newcommand {\LWRsiunitxprintdecimal }[1]{\LWRsiunitxprintdecimalsub #1,,,\LWRsiunitxEND }\)

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.

This part of the thesis focuses on the theoretical formulation of transition rates \(\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.

Ideally, all dissipative and inelastic couplings leading to changes of the states in both *subsystems*^{1} 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{equation} \hat {H}=\hat {H}_\mathrm {S}+\hat {H}_\mathrm {B}+\hat {V}_\mathrm {SB}=\hat {H}_0+\hat {V}_\mathrm {SB}, \end{equation}

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{equation} i\hbar \frac {\partial }{\partial t}\hat {\rho }_\mathrm {S}=i\hbar \mathcal {L}\hat {\rho }_\mathrm {S}(t) =\left [\hat {H}_\mathrm {0},\hat {\rho }_\mathrm {S}(t)\right ]+\int _0^td\tau \mathcal {L}_\mathrm {D}(t,\tau ){\hat {\rho }_\mathrm {S}(\tau )}. \label {eq:chapter1_LvN} \end{equation}

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^{2}. 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{equation} \frac {\partial }{\partial t}\hat {\rho _\mathrm {S}}(t)=\mathcal {L}_\mathrm {S}\hat {\rho }_\mathrm {S}+\mathcal {L}_\mathrm {D}\hat {\rho }_\mathrm {S}=-\frac {i}{\hbar }\left [\hat {H}_\mathrm {S},\hat {\rho }_\mathrm {S}\right ]+\mathcal {L}_\mathrm {D}\hat {\rho }_\mathrm {S}, \label {eq:chapter1_Lb-LvN} \end{equation}

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 \(\ket {\phi _n}\) of the *system* Hamiltonian, \(\rho _{nm}=\mel {\phi _n}{\hat {\rho }}{\phi _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

\(\seteqnumber{0}{4.}{3}\)\begin{equation} \mathcal {L}_\mathrm {D}\hat {\rho }_\mathrm {S}=\sum _k\big (\hat {C}_k\hat {\rho }_\mathrm {S}\hat {C}_k^\dagger -\frac {1}{2}\left [\hat {C}_k^\dagger \hat {C}_k,\hat {\rho }_\mathrm {S}\right ]_+\big ), \end{equation}

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{equation} \hat {C}_k=\sqrt {\Gamma _{i,f}}\ket {\phi _{f}}\bra {\phi _{i}}, \end{equation}

which is the projection of an initial state onto a final state weighted with a rate \(\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

\(\seteqnumber{0}{4.}{5}\)\begin{equation} \Gamma _{i,f}=\frac {2\pi }{\hbar }\sum \limits _{I,F}f_{I}(T)\left [1-f_{F}(T)\right ]\cdot |\!\mel {\phi _{f}\Phi _F}{\hat {V}}{\phi _{i}\Phi _F}\!|^2\delta (E_{i,I}-E_{f,F}). \label {eq:chapter1_FGR} \end{equation}

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 \(\ket {\phi _i}\) and \(\ket {\phi _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 \(\Gamma _{i,f}\) is the sum of all individual contributions:

\(\seteqnumber{0}{4.}{6}\)\begin{equation} \Gamma _{i,f}^\mathrm {tot}=\Gamma _{i,f}^\mathrm {vib}+\Gamma _{i,f}^\mathrm {dip,inel}+\Gamma _{i,f}^\mathrm {res,inel} \label {eq:chapter1_Rtot} \end{equation}

^{1} The Si–H complex and its environment

^{2} i.e. it is trace preserving and completely positive by mathematical construction

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^{3}, however, an analytic fit of the following form has been used:

\begin{equation} V_\mathrm {g}(q)= v_0q^4-v_2q^2+\sigma \sqrt {\frac {v_0}{2v_2}}q, \label {eq:chapter1_analytic_pot} \end{equation}

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=\SI {-3.2}{a_0\sqrt {m}}\) and \(q=\SI {3.0}{a_0\sqrt {m}}\) with a barrier of \(\SI {2.7}{eV}\) 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 \SI {1.7}{eV}\) (\(\sim \SI {2.1}{eV}\)) and a shifted equilibrium position by \(\sim \SI {0.9}{a_0\sqrt {m}}\) (\(\sim \SI {0.4}{a_0\sqrt {m}}\)), see Fig. 4.1 and Tab. 4.1a. Note that the positive potential \(V^+(q)\) additionally exhibits a shifted transition point.

\(v_0\) | \(v_2\) | \(\sigma \) | |

\(V^0(q)\) | 0.023 | 0.442 | 1.062 |

\(V^-(q)\) | 0.046 | 0.459 | 0.955 |

\(V^+(q)\) | 0.072 | 0.640 | 1.241 |

\(n\) | \(E_{n}\) (meV) | \(E_{n,\mathrm {min}}\) (meV) | \(E_{n,\mathrm {max}}\) (meV) |

\(\phi _\mathrm {0,L}\) | 84.44 | 73.22 | 113.54 |

\(\phi _\mathrm {1,L}\) | 251.98 | 219.34 | 309.56 |

\(\phi _\mathrm {2,L}\) | 417.25 | 386.72 | 491.47 |

⋮ | ⋮ | ⋮ | ⋮ |

\(\phi _\mathrm {6,L}\) | 1054.09 | 597.78 (\(\phi _\mathrm {4,L}\)) | 1724.06 (\(\phi _\mathrm {10,L}\)) |

\(\phi _\mathrm {7,R}\) | 1138.60 | 693.31 (\(\phi _\mathrm {5,R}\)) | 1788.93 (\(\phi _\mathrm {11,R}\)) |

⋮ | ⋮ | ⋮ | ⋮ |

\(\phi _\mathrm {30,L}\) | 2627.90 | 2131.35 (\(\phi _\mathrm {26,L}\)) | 3154.57 (\(\phi _\mathrm {34,L}\)) |

\(\phi _\mathrm {31,R}\) | 2675.69 | 2127.86 (\(\phi _\mathrm {25,R}\)) | 3266.36 (\(\phi _\mathrm {35,L}\)) |

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.

^{3} In order to account for a distribution of paramaters, see Sec. 3.2.

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 \(\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^{4}. 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{equation} R_i=\sum _{j=1}^{3N}O_{i,j}r_j,\quad \mathrm { with}\quad \sum _{i=1}^{3N}O_{i,k}O_{i,l}=\delta _{k,l}, \end{equation}

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{equation} \hat {H} = \hat {H}_S+\hat {H}_{SB}+\hat {H}_{B}, \end{equation}

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

\begin{equation} \begin{aligned} \hat {H}_S &= -\frac {\hbar ^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}\Lambda _{i,j}(\mathbf {q})Q_iQ_j,\\ \hat {H}_B &= \sum _{i=1}^{3N-M}\big (-\frac {\hbar ^2}{2}\frac {\partial ^2}{\partial Q_i^2}+\frac {1}{2}\omega _i^2 Q_i^2\big ). \end {aligned} \label {eq:chapter1_dpt_hamiltonians} \end{equation}

The one– and two–phonon coupling functions \(\lambda _i(\mathbf {q})\) and \(\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

\(\seteqnumber{0}{4.}{11}\)\begin{equation} \begin{aligned} \lambda _i(\mathbf {q})&=\frac {\partial V}{\partial Q_i}\Bigr |_{\mathbf {q},\mathbf {Q}=0}=\sum _{k=1}^{3N}O_{M+i,k}\frac {\partial V}{\partial r_k}\Bigr |_{\mathbf {q},\mathbf {Q}=0},\\ \Lambda _{i,j}(\mathbf {q})&=\frac {\partial ^2 V}{\partial Q_i\partial Q_j}\Bigr |_{\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}\Bigr |_{\mathbf {q},\mathbf {Q}=0}. \end {aligned} \label {eq:chapter1_dpt_derivatives} \end{equation}

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

\begin{equation} \Gamma _{i,f} = \frac {2\pi }{\hbar }|\!\mel {\Phi _{f}}{\hat {H}_\mathrm {SB}}{\Phi _{i}}\!|^2\delta (E_i-E_f), \label {eq:chapter1_dpt_fgr} \end{equation}

where one can choose the initial and final states \(\Phi _i\) and \(\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\neq \SI {0}{K}\)) 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{equation} \begin{aligned} \Gamma _{i,f}^{(1),\downarrow }&=\pi \sum _k|\!\mel {\phi _{i}}{\lambda _k(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k\right \rangle +1}{\omega _k}\delta (\varepsilon _i-\varepsilon _f-\hbar \omega _k)\\ \Gamma _{i,f}^{(1),\uparrow }&=\pi \sum _k|\!\mel {\phi _{i}}{\lambda _k(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k\right \rangle }{\omega _k}\delta (\varepsilon _i-\varepsilon _f+\hbar \omega _k), \end {aligned} \label {eq:chapter1_dpt_onePhonon} \end{equation}

where \(\left \langle n_k\right \rangle =\sum _n np_{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 _\mathrm {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{equation} \begin{aligned} \Gamma _{i,f}^{(2a),\downarrow }&=\frac {\pi \hbar }{8}\sum _{k\ne l}|\!\mel {\phi _{i}}{\Lambda _{k,l}(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k\right \rangle +1}{\omega _k}\frac {\left \langle n_l\right \rangle +1}{\omega _l} \delta (\varepsilon _i-\varepsilon _f-\hbar \omega _k-\hbar \omega _l),\\ \Gamma _{i,f}^{(2a),\uparrow }&=\frac {\pi \hbar }{8}\sum _{k\ne l}|\!\mel {\phi _{i}}{\Lambda _{k,l}(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k\right \rangle }{\omega _k}\frac {\left \langle n_l\right \rangle }{\omega _l} \delta (\varepsilon _i-\varepsilon _f+\hbar \omega _k+\hbar \omega _l) \end {aligned} \label {eq:chapter1_dpt_twoPhonon_1} \end{equation}

whereas the expression for two identical bath modes is

\(\seteqnumber{0}{4.}{15}\)\begin{equation} \begin{aligned} \Gamma _{i,f}^{(2b),\downarrow }&=\frac {\pi \hbar }{8}\sum _{k}|\!\mel {\phi _{i}}{\Lambda _{k,k}(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k^2\right \rangle +3\left \langle n_k^2\right \rangle +2}{\omega _k^2} \delta (\varepsilon _i-\varepsilon _f-2\hbar \omega _k)\\ \Gamma _{i,f}^{(2b),\uparrow }&=\frac {\pi \hbar }{8}\sum _{k}|\!\mel {\phi _{i}}{\Lambda _{k,k}(\mathbf {q})}{\phi _{f}}\!|^2\frac {\left \langle n_k^2\right \rangle -\left \langle n_k^2\right \rangle }{\omega _k^2} \delta (\varepsilon _i-\varepsilon _f+2\hbar \omega _k) \end {aligned} \label {eq:chapter1_dpt_twoPhonon_2} \end{equation}

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

\(\seteqnumber{0}{4.}{16}\)\begin{equation} \delta (E)\rightarrow \frac {1}{\pi }\frac {\gamma }{\gamma ^2+E^2}, \label {eq:chapter1_lorentzians} \end{equation}

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 \(\bm {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

\(\seteqnumber{0}{4.}{17}\)\begin{equation} \mathbf {H}^\prime =\mathbf {U}\mathbf {H}\mathbf {U}^T, \end{equation}

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

\(\seteqnumber{0}{4.}{18}\)\begin{equation} \mathbf {e}_Q=\mathbf {U}^T\mathbf {e}_Q^\prime , \end{equation}

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– (\(\Gamma ^{(1)}\)) and two–phonon transition (\(\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:

\(\seteqnumber{0}{4.}{19}\)\begin{equation} \tau _1^{-1} = \big (\Gamma _{1,0}^{(1),\downarrow }+\Gamma _{1,0}^{(2),\downarrow }\big )-\big (\Gamma _{0,1}^{(1),\uparrow }+\Gamma _{0,1}^{(2),\uparrow }\big ). \end{equation}

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 \(\Gamma _{1,0}^{(1)}\) exhibits a strong
dependence between \(\gamma \in [1,10]\mathrm {cm}^{-1}\) (\(\tau _1^{(1)}=\SI {18.3}{ns}-\SI {1.89}{ns}\)), \(\Gamma _{1,0}^{(2)}\) is stable against changes in \(\gamma \) (\(\tau _1^{(2)}=\SI
{152.9}{ps}-\SI {97.5}{ps}\)). 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.

Mode | T [K] | \(\tau _1^\mathrm {total}\) | \(\Gamma _{1,0}^{(1)}\) | \(\tau _1^{(1)}\) | \(\Gamma _{1,0}^{(2\mathrm {a})}\) | \(\Gamma _{1,0}^{(2\mathrm {b})}\) | \(\tau _1^{(2)}\) |

0 | 0.1237 | 0.1930 | 5.1801 | 7.3745 | 0.5101 | 0.1268 | |

Si–H\(_\mathrm {break}^{1\rightarrow 0}\) | 300 (\(\downarrow \)) | 0.1180 | 0.2116 | 4.9531 | 7.7452 | 0.5244 | 0.1209 |

300 (\(\uparrow \)) | 0.0098 | 0.0009 | 0.0001 | ||||

Si–H\(_\mathrm {break}^{2\rightarrow 1}\) | 0 | 46.344 | 0.4238 | 2.359 | 19.9757 | 1.1779 | 47.273 |

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 \(\SI {5}{\per \cm }\times \SI {5}{\per \cm }\). 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 \(\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 \(\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 =\SI {5}{\per \cm }\) has been chosen as the reference here^{5}, which
yields \(\tau _\mathrm {0K}=\SI {0.126}{ns}\) (\(\tau _\mathrm {300K}=\SI {0.120}{ns}\)).

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 \(\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 \(\Delta E_{2,0}>2\hbar \omega _\mathrm {max}\), which is beyond the two–phonon coupling included^{6}. Therefore, for all practical applications an idealized harmonic model [101, 105, 270, 276] can be employed where the selection rule \(\Delta n=1\) applies. Such
a harmonic model possesses a linear scaling law \(\Gamma _{i,i-1}^{\downarrow } = i\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\rightarrow 1\) transition have been calculated, see Table 4.2. Since the energy difference \(\Delta E_{2,1}\) is still outside
the Si/SiO2 phonon band, \(\Gamma _{2,1}^{(2)}\) again dominates the resulting lifetime. Furthermore, it is very likely that the direct transition \(2\rightarrow 0\) is less important, due to the fact that at least three bath
modes are necessary to dissipate the energy^{7}. Finally, one can see in Tab. 4.2 that the resulting lifetime of \(\tau
_2=\SI {46}{ps}\) agrees well with the simple model prediction of \(\tau _2=\tau _1/2=\SI {61.5}{ps}\).

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 \(\ket {\phi _n}\) is then given by

\begin{equation} \tau _n=\frac {\tau _1}{n}, \label {eq:chapter1_vibrational_rates} \end{equation}

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} Note that in a fully rigorous approach the 2D potential energy surface (PES \(V(\bm {q})\)) of the Si–H bond would be considered, which effectively enables coupling between the low–energy bending and the
high–energy stretching modes, see [270]. However, this requires a well–parametrized 2D PES which is outside of the scope of this work.

^{5} The calculations show an approximated deviation of \(25\%\) compared to \(\gamma =\SI {1}{\per \cm }/\SI {10}{\per \cm }\).

^{6} Or by including excited vibrational states of the 2D potential surface, see 3.2.

^{7} The inclusion of a three phonon term in the system–bath Hamiltonian is proportional to \(Q_iQ_jQ_k\) which would render the calculation of the coupling element \(\kappa _{i,j,k}\) computationally
extremely challenging.

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 \(\bm {\mu }\) and can be
calculated by a formula derived by Persson *et al.* [277–279]

\begin{equation} \Gamma _{i,f}^\mathrm {dip}=\frac {I}{e}\abs {\frac {\mel {\phi _f}{\bm {\mu }}{\phi _i}}{e a_0}}^2, \label {eq:chapter1_dipole_rates} \end{equation}

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 \(\Gamma ^\mathrm {dip}_{i,f}\) for \(i\Rightarrow i+1\) and \(i\Rightarrow i+2\) transitions have been calculated for different currents and
compared to the vibrational upward rates \(\Gamma ^\mathrm {vib}_{i,f}\)^{8} at \(T=\SI {300}{K}\), see Fig. 4.4. Current densities in scaled MOSFETs already exceed \(\SI {100}{nA/nm^2}\), which is of similar order as tunneling currents in
Si–H related STM experiments (\(\SI {1}{}-\SI {10}{nA}\)), assuming that the STM electron beam is of atomic dimensions. One can see that only at higher currents the respective upward rates for \(\Delta n=+1\) are mainly
determined by dipole induced scattering, whereas transitions to higher excited states (\(\Delta n=+2\)) are already orders of magnitude lower. Although Fig. 4.4 actually shows a non–negligible contribution of \(\Gamma ^\mathrm {dip}_{i,f}\), 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.

^{8} \(\Gamma ^\mathrm {vib}_{i,f}\) are only available for \(\Delta n+1/\Delta n-1\) transitions, see Sec. 4.1.3.

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 \(\epsilon \)
can induce an excitation from state \(\ket {\phi _{i}}\) to \(\ket {\phi _{f}}\) via the vibrational eigenstate \(\ket {\psi _{j}}\) of the resonance with

\begin{equation} \Gamma _{i,f}^\mathrm {res}=\frac {4\Delta _\mathrm {res}^2}{\pi \hbar }\int d\epsilon f(\epsilon )\left [1-f(\epsilon )\right ]\cdot |\sum \limits _j\frac {\braket {\phi _{f}}{\psi _{j}}\braket {\psi _{j}}{\phi _{i}}}{\epsilon _\mathrm {res}-\epsilon +\epsilon _i-\epsilon _j+i\Delta _\mathrm {res}}|^2. \label {eq:chapter1_resonance_rates} \end{equation}

This expression includes the energetic position of the resonance with respect to the Fermi level, \(\epsilon _\mathrm {res}\), its width as the inverse of the resonance lifetime, \(\Delta _\mathrm {res}=\hbar /\tau _\mathrm {res}\), as well as the energies of the eigenstates \(\ket {\phi _i}\) and \(\ket {\psi _j}\). Furthermore, it accounts for the energy distribution function \(f(\epsilon )\) 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.

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 \(\ket {\phi _i}\) to \(\ket {\phi _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{equation} \frac {\partial \rho _\mathrm {S}}{\partial t}=\mathcal {L}_\mathrm {S}\rho _\mathrm {S}+\mathcal {L}_\mathrm {D}\rho _\mathrm {S}=0, \label {eq:chapter1_quasi_equilibrium} \end{equation}

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{equation} \Gamma ^\mathrm {break}=\sum \limits _{i}\Gamma _{i,\mathrm {WKB}}\,P_{i}=\sum \limits _{i}\exp \big (-\frac {2}{\hbar }\int \limits _{x_{1,i}}^{x_{2,i}}\sqrt {2(V-E_i)} \mathrm {d}x\big )P_i, \label {eq:FinalDesorbRate} \end{equation}

with \(P_{i}\) being the population of the eigenstate \(i\) and \(\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 \(\ket {\phi _\mathrm {R}}\) in the above calculation. A more detailed description of how to calculate the final rate \(\Gamma ^\mathrm {break}\) is given in Appendix E.

Knowing the bond–breaking rate \(\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

\(\seteqnumber{0}{4.}{25}\)\begin{equation} \frac {\mathrm {d}\left [\mathrm {SiH}\right ]}{\mathrm {d} t}=-\Gamma ^\mathrm {break}\left [\mathrm {SiH}\right ]+\Gamma ^\mathrm {pass}\left [\mathrm {H}_2\right ]\left [P_\mathrm {b}\right ], \label {eq:ReacEqu} \end{equation}

where \(\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 \(\Gamma ^\mathrm {pass}=\Gamma _0^\mathrm
{pass}\mathrm {exp}(-E_\mathrm {pass}/k_\mathrm {B}T)\) superimposed with a Gaussian distribution of the passivation energy \(E_\mathrm {pass}\)^{9}. The inferred parameters are \(E_\mathrm {p}\sim \SI
{1.51}{eV}\), \(\sigma _{E_\mathrm {p}}~\sim \SI {0.12}{eV}\) and \(\Gamma _0^\mathrm {pass}\sim \SI {5 e-6}{cm^3\per \second }\)^{10}, 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{equation} \left [\mathrm {SiH}\right ]+\left [\mathrm {H}\right ]=\left [\mathrm {H}\right ]_\mathrm {tot},\\ \label {eq:MaxConcEqu} \end{equation}

while the value for [H\(_2\)] was set to a constant concentration of \(\SI {5 e17}{cm^{-3}}\), given by the physical solubility of H\(_2\) in vitreous silica [286]. Introducing the quantity \(\mathrm {N}_0\mathrel {\widehat {=}}\left [\mathrm {SiH}\right ]_\mathrm {max}\) and solving for \(f_\mathrm {{Pb}}\mathrel {\widehat {=}}\left [\mathrm {P_b}\right ]/\mathrm {N_0}=1-\left [\mathrm {SiH}\right ]/\mathrm {N_0}\), the probability that a bond is broken is then given by

\(\seteqnumber{0}{4.}{27}\)\begin{equation} \frac {\partial f_\mathrm {{Pb}}}{\partial t}=\left (1-f_\mathrm {{Pb}}\right )\Gamma ^\mathrm {desorb}-f_\mathrm {{Pb}}{\Gamma }^\mathrm {pass}[\mathrm {H_2}]. \label {eq:ReacEquSol} \end{equation}

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.

^{9} Recent studies show that also an applied gate bias substantially influences the recovery dynamics [249, 250]. It is assumed that the Fermi level affects the defects charge state and hence the passivation barrier.

^{10} Note that while the reported passivation barrier \(E_\mathrm {pass}\) is consistent throughout the literature, the reaction rate constant \(\Gamma _0^\mathrm {pass}\) varies between \(\SI {e-4}{}\) and
\(\SI {e-9}{cm^3\per \second }\).

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 \(\epsilon \) triggers a vibrational excitation from state \(\ket
{\phi _{i}}\) to \(\ket {\phi _{f}}\) in the ground state \(V(q)\) via the eigenstate \(\ket {\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{equation} \Gamma _\mathrm {MP/SP}=\int _{\epsilon _\mathrm {th}}f(\epsilon )g(\epsilon )v(\epsilon )\sigma _\mathrm {MP/SP}(\epsilon )\mathrm {d}\epsilon , \label {eq:chapter1_hcd_tyaginov} \end{equation}

where \(f(\epsilon )\) is the EDF, \(g(\epsilon )\) the DOS, \(v(\epsilon )\) corresponds to the carrier velocity and \(\sigma (\epsilon )\) 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 \(\epsilon _\mathrm {th}\). The energy dependent cross section is assumed to follow an empirical Keldysh–like expression \(\sigma (E)=\sigma _0(\epsilon -\epsilon _\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

\(\seteqnumber{0}{4.}{29}\)\begin{equation} \Gamma _{i,f}=\frac {4\Delta _\mathrm {res}^2}{\pi }\int _\epsilon f(\epsilon )g(\epsilon )v(\epsilon )\sigma _\mathrm {eff}(\epsilon )\mathrm {d}\epsilon , \label {eq:chapter1_resonance_rates_approximated} \end{equation}

with the effective cross section

\(\seteqnumber{0}{4.}{30}\)\begin{equation} \sigma _\mathrm {eff}(\epsilon )=\sigma _0\big |\sum \limits _j\frac {\braket {\phi _{f}}{\psi _{j}}\braket {\psi _{j}}{\phi _{i}}}{\epsilon _\mathrm {res}-\epsilon +\epsilon _i-\epsilon _j+i\Delta _\mathrm {res}}\big |^2. \label {eq:chapter1_effective_cross_section} \end{equation}

In analogy to (4.23) it is the probability of the vibrational transition from \(\ket {\phi _i}\) to \(\ket {\phi _f}\) in the ground state potential via the resonance eigenstate \(\ket {\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 \(\Delta _{n,n+1}\) (ladder climbing, MP mechanism), overtone transitions \(\Delta _{n,m}\) as well as transitions directly causing dissociation \(\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 \(\epsilon _\mathrm {res}+\epsilon _i-\epsilon _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 \(\bm {\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)

\(\seteqnumber{0}{4.}{31}\)\begin{equation} \Gamma _{i,f}^\mathrm {tot}=\Gamma _{i,f}^\mathrm {res}+\Gamma _{i,f}^\mathrm {vib}. \end{equation}

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 \(\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.