The degradation mechanisms discussed in this text proceed via chemical or electrochemical reactions at defects in the oxide or at the semiconductor-oxide interface. This chapter reviews the theories of theoretical and physical chemistry that are relevant for the applications presented in the later chapters, starting from the highest level of detail and moving to the less detailed but computationally more efficient macroscopic descriptions. Most of the concepts presented in this chapter are much more generally applicable than just to describe the dynamics of oxide defects. Examples are given to link those concepts to the context of the present application.

2.1 Problem Statement: The Molecular Hamiltonian

2.2 The Born-Oppenheimer Approximation

2.3 Electronic Structure Methods

2.3.1 Self-Consistent Field Methods: Hartree-Fock

2.3.2 Self-Consistent Field Methods: Density Functional Theory

2.3.3 Explicit Many-Body Methods

2.3.4 Empirical Methods

2.4 The Nuclear Problem: Chemical Microstates

2.5 Chemical States and Reactions

2.6 The Chemical Master Equation

2.7 The Calculation of Rates

2.8 Barrier Hopping Transitions

2.9 Multi-Phonon Transitions

2.9.1 Vibronic Coupling

2.9.2 Quantum Mechanical Theory of Vibronic Transitions

2.9.3 Model Matrix Elements

2.9.4 The Line Shape Function

2.9.5 Line Shapes for Non-Radiative Transitions

2.10 Vibronic Transitions with Classical Nuclei

2.11 Solution of the Master Equation

2.11.1 Subsystems, Well-Stirredness, and Diffusion

2.12 Reaction Rate Equations

2.2 The Born-Oppenheimer Approximation

2.3 Electronic Structure Methods

2.3.1 Self-Consistent Field Methods: Hartree-Fock

2.3.2 Self-Consistent Field Methods: Density Functional Theory

2.3.3 Explicit Many-Body Methods

2.3.4 Empirical Methods

2.4 The Nuclear Problem: Chemical Microstates

2.5 Chemical States and Reactions

2.6 The Chemical Master Equation

2.7 The Calculation of Rates

2.8 Barrier Hopping Transitions

2.9 Multi-Phonon Transitions

2.9.1 Vibronic Coupling

2.9.2 Quantum Mechanical Theory of Vibronic Transitions

2.9.3 Model Matrix Elements

2.9.4 The Line Shape Function

2.9.5 Line Shapes for Non-Radiative Transitions

2.10 Vibronic Transitions with Classical Nuclei

2.11 Solution of the Master Equation

2.11.1 Subsystems, Well-Stirredness, and Diffusion

2.12 Reaction Rate Equations

The theory employed in the present work rests upon the fundament of physical chemistry,
sometimes also referred to as “low energy physics”, which is applicable to the processes
dominating our everyday experience like the dynamics of fluids and gases, electrical and
thermal conduction, and chemical reactions. In contrast to high energy physics, where
subatomic particles and at least three of the four fundamental interactions have to be
considered, the central actors in physical chemistry are electrons and nuclei which only interact
through electrostatic forces. The theoretical foundation of quantum chemistry thus lies in the
Schrödinger equation of the system of electrons and nuclei. In the following, and denote
the position vector of all nuclei and all electrons, respectively. For a system containing _{a} and _{A} and the associated masses are
_{A} for the nuclei (which is
different for each species). In the following, the term “vibrational” will be used for the
nuclear degrees of freedom. This refers to the nature of the nuclear motion, which is usually
oscillatory and confined to a potential well of some sort, as will be explained later. Especially in
older works [88] and generally in works on molecules the nuclear degrees of freedom are
split into vibrational and rotational parts, where the latter mean a total rotation of the
molecule. As this thesis is concerned with phenomena in a solid-state context, the internal
degrees of freedom (center-of-mass position and rotation) present in theories on molecules are
omitted [89].

The Schrödinger equation for the system of electrons and nuclei is [88, 90, 91]

| (2.1) |

where Ψ(

| (2.2) |

is the Hamiltonian of the system, which is sometimes also termed the

The Born-Oppenheimer approximation relies on the fact that the mass of even the lightest nuclei exceeds the electron mass by orders of magnitude. In order to make use of this property, Born and Oppenheimer rewrote the molecular Hamiltonian as

| (2.8) |

with

= | (2.9) | |

_{0} | = | (2.10) |

_{1} | = _{A=1}^{N}_{
A} | (2.11) |

| (2.12) |

Although this is an eigenproblem in the electronic degrees of freedom, the results depend
parametrically on the position of the nuclei. The higher order elements of the series expansion show
that the expectation value of the energy

| (2.13) |

where

| (2.14) |

This approximation implicitly assumes that the electronic quantum number is not affected by the
momentum of the nuclei, which means that the nuclei are assumed to move slowly enough that the
electrons can adjust

The original derivation of Born and Oppenheimer [88, 90] was restricted to molecules in which the potential energy surface was essentially a quadratic function with the nuclei moving around the minimum. Interestingly, it also gives very accurate results for situations where these conditions are not even approximately fulfilled [92, 90]. This behavior and also the breakdown of the approximation can be understood using an improved approach that is discussed later in Sec. 2.9.1.

In summary, the essence of the Born-Oppenheimer approximation is that due to the large difference in masses, the heavy nuclei are standing still from the perspective of the electrons. The electrons will thus adjust their orbit instantaneously to the motion of the nuclei and the energy of this orbit acts back as a potential on the nuclei. The resulting separation of the electronic and the nuclear degrees of freedom is, however, only a slight improvement from the perspective of solubility, as a typical solid or molecule still contains too many electrons and ions for an exact solution of (2.12) and (2.14). However, both the electronic as well as the nuclear problem can be greatly simplified by the exploitation of inherent symmetries and some well-defined approximations.

The solution of the electronic many-particle problem arising from the Born-Oppenheimer
approximation is a classic topic of quantum chemistry, and a multitude of approaches exist using
different approximations. Depending on the size of the system and the available computational
resources, different levels of sophistication can be employed to calculate the potential energy surface of
a molecule or a solid. While the simplest approaches employ analytic potentials which try to mimic
the shape of the electronic energy, the more sophisticated approaches actually attempt to solve the
electronic Schrödinger equation (2.12) for its ground (and sometimes also for an excited) state. The
latter methods are referred to as

For all but the smallest systems, an exact solution of the electronic Schrödinger equation is
computationally unfeasible, as the dimensionality of the problem increases linearly with the number of
particles, leading to an exponential increase in computational complexity. The first approach to
the solution of the many-electron problem for isolated atoms was reported by Hartree
in three papers in 1928 [93, 94, 95]. His approach, described in the second part of the
series, was termed ‘self-consistent field method’ as it approximates the electron-electron
interaction through an effective potential _{H}, calculated from the total charge density less the
charge density of the electron under consideration in order to avoid self-interaction (see
Fig. 2.2).

| (2.15) |

The formulation of Hartree’s method initially followed physical intuition rather than a mathematical derivation. Its connection to the many-body Schrödinger equation was clarified later as a variational determination of the ground state using a product-wave function ansatz.

| (2.16) |

A central problem of this ansatz is that it does not have the required fermionic symmetries upon particle exchange [96, 97, 98]. An improved description was found by Fock using a properly symmetrized ansatz, which can be written in a compact form as determinant in which the single-particle wave function is the same along a column and the particle coordinate is the same along a row.

| (2.17) |

Following the same basic steps as in the Hartree-method again leads to an effective single-particle
problem, with an additional self-consistent non-local potential

| (2.18) |

The resulting method was named Hartree-Fock self-consistent field
method, the new term in the equation is usually called exchange- or
Fock-operator^{1} .
The Hartree-Fock self-consistent field method gives good results in many situations and is a standard
method of quantum chemistry. The wave function resulting from the method as well as the total
(ground state) energy are physically relevant within the bounds of the approximation, in contrary
to the density functional approach which is discussed below. Although it is considered
a reasonable starting point, the accuracy of the theory is limited due to the neglect of
electronic correlation [96, 97]. A popular method that improves on this situation is the
Møller-Plesset perturbation series, which brings correlation into the Hartree-Fock-Hamiltonian
using a perturbation approach that dramatically improves the accuracy for some systems
[99, 97].

Despite their sound theoretical foundation and the good accuracy of the perturbation
approach, Hartree-Fock based methods are rarely used for large molecules or solid-state
problems due to the large computational effort arising from the non-local exchange operator
(^{4}), where

Density functional theory is a more mathematical approach to the many-body problem. It is based on the findings by Hohenberg and Kohn, which state that

- for a given external potential
V ext which acts equally on all of the electrons, there is a ground state electron densityn _{0}. This density is unique to the external potential apart from a trivial constant energy. - there is a functional
E [n ] of the electron density which, when minimized inn for a givenV ext , gives the ground state energyE _{0}and the ground state electron densityn _{0}of the many-body system.

From these two findings it follows that the exact ground state energy of the (3

| (2.19) |

that has the same ground state density and energy as the interacting system. All exchange and
correlation effects of the interacting system are represented in the auxiliary system through the
effective local potential ^{3})) compared to the Hartree-Fock approach, and the remarkable
accuracy of the method for a broad set of problems, DFT has become the standard approach in many
fields of physical chemistry. The exchange-correlation potential has been a focus of active research for
decades, yielding a manifold of different local and non-local potentials for different specific
applications. A standard general-purpose functional that is also used extensively in the present work is
the first-principles gradient-corrected functional by Perdew, Burke and Ernzerhof [101]. It is evident
from the explanation above that, in contrary to the HF method, the wave functions in Kohn-Sham
DFT do not have any physical relevance. Rather, the wave functions and eigenlevels are
features of an auxiliary system that is only introduced as part of an approximation to the
many-body problem. Nevertheless, it is common to draw conclusions from the Kohn-Sham
eigenlevels, e.g. when calculating band structures [96] or defect levels [102]. Problems of the
Kohn-Sham eigenspectrum, such as the underestimation of band-gaps, as well as errors in total
energies are sometimes reduced by the application of so-called

The most accurate methods available today are the diffusion and variational quantum Monte Carlo methods, and the configuration-interaction (CI) or coupled-cluster (CC) methods [96, 97]. These methods are able to appropriately handle correlation effects and are often used as a high-accuracy reference to benchmark self-consistent field methods. However, due to their vast computational demand those methods are usually limited to very small molecules or require massively parallel computers with thousands of CPUs. For the present work, these methods are way too expensive and also more accurate than necessary, considering the various approximations that have to be made.

The broad class of empirical electronic structure methods offers inexpensive, yet accurate solutions to various problems, especially in the field of microelectronics research. The level of physical accuracy here ranges from the classic semi-empirical quantum chemistry methods (CNDO, MNDO, PM3, …) [97] over tight-binding methods, which use a Hamiltonian that is built from analytic expressions [96], to the simplest empirical molecular dynamics potentials, which are based on analytical functions of the positions of the nuclei.

Semi-empirical quantum chemistry methods use a Hartree-Fock-based Hamiltonian, which replaces
some of the more expensive integrals with parametrizable expressions. These methods have been under
active development for decades, and are especially popular in the field of organic chemistry. The
parameters for these methods are usually calibrated to small molecules, and the standard
parametrizations have been shown to fail in the solid-state context [103, 104]. A comparison of
semi-empirical methods for the calculation of the oxygen vacancy defect in

Empirical tight-binding methods are often used for the study of the electronic structure of ideal
bulk or surface structures [96]. They have also been employed in defect studies, primarily for defects
in crystalline and amorphous silicon [106, 107, 108], but also for defects in SiO

A review of different empirical molecular dynamics potentials applied to silicon or SiO

The advantage of semiempirical methods is the inclusion of experimental input, to account
efficiently for the physics neglected in the approximations. For the study of defects in SiO

For electronic structure problems it is reasonable to treat the correlated motion of the electrons as a
second-order effect, which does not influence the basic behavior for most practical cases. For the
motion of the nuclei, an approximation like this is not possible. Vibrational states are usually heavily
correlated motions of nuclei like stretching or bending modes in molecules or phonons in a solid. A
fully quantum-mechanical treatment is again impossible due to the large number of degrees of freedom
and even a quantum Monte-Carlo-like treatment is unfeasible due to the computational effort required
to calculate the potential energy surface. It is thus necessary to either approximate the quantum
mechanical nuclear problem as a classical one or to approximate the potential energy surface as an
essentially quadratic function of the position of the nuclei so that the atomic motion can be
separated into uncoupled harmonic oscillations. The former method is broadly applied
throughout theoretical chemistry for molecules and solids and is well suited for situations where
temperatures are high enough so that quantization effects may be neglected and where
the adiabatic Born-Oppenheimer approximation holds. With respect to the BTI model
introduced in Sec. 1.5, this concerns the structural transitions 1

Once the potential energy surface for a given atomic system is known and a sufficiently accurate
description of the nuclear dynamics has been found, it is possible to calculate observable
quantities and predict the evolution of that system over time. The state of the atomic
system under consideration is fully described by the state of the electronic system, which is
usually identified by its potential energy surface, and the nuclear state, which is

In order to understand chemical kinetics, it is first necessary to define in a physically
meaningful way the ^{2} .
In the case of a defect within a semiconductor device the corresponding observable may be for
example its charge state, which is observed through its effects on the device characteristics, or its
bonding state, which will give rise to different capture and emission times. The microstates belonging
to a chemical state are consequently those which are indistinguishable from the viewpoint of the
associated observable. It is the nature of these observables that they need to be stable for a sufficient
amount of time, so that they can be measured by a human-made apparatus. As for the microstates,
the electronic subsystem often defines its own chemical states as it tends to establish a stationary state
on a very short time scale (femtoseconds) and this state is then quite stable in the absence of
potential energy surface crossings or radiative transitions. For the vibrational subsystem, the
requirement on stability over a certain time-range means that states are defined as the
neighborhood of minima in the potential energy surface. This follows from the nature of
the atomic motion, which is essentially random due to the permanent interaction with
the environment. Usually within picoseconds after a change in state the momenta of the
nuclei will assume the thermal equilibrium distribution induced by the heat-bath. Thus the
momentum, or in the quantum mechanical view the vibrational quantum state, of the
nuclei is not relevant in the definition of chemical states and the only relevant information
comes from the shape of the region of the configuration space in which they are moving,
which is defined by the energetic barriers of the potential energy surface, as illustrated in
Fig. 2.3.

Once a chemical state is defined in terms of microstates, the measurable observables are obtained from thermal averages over those microstates belonging to this chemical state. If the nuclei are treated as classical particles, those thermal averages consist of integrals over a certain potential energy surface and over the momentum-space of the nuclei, which is commonly assumed to follow the behavior of the ideal gas [112].

As mentioned above, in practical calculations it is assumed that transitions between chemical states occur instantaneously and the thermal equilibrium of the microstates in the new chemical state is established immediately. This assumption is obviously justified for transitions of the electronic subsystem, as mentioned above. Transitions between the minima of a potential energy surface or between two different electronic states are usually followed by an equilibration phase where the momentum distribution of the nuclei returns to the ideal gas distribution. This equilibration usually proceeds at time-scales of picoseconds, which is also reasonably instantaneous for most observables. However, especially for the case of hydrogen-passivated dangling bonds at the silicon surface, long-term stable vibrational modes have been reported which lead to a considerably prolonged equilibration [113]. Thus, regarding the nuclear subsystem the assumption of instantaneous transitions needs to be applied cautiously.

For the sake of completeness it is mentioned here that, as chemical states may be arbitrarily defined,
it is possible that certain microstates are part of more than one chemical state. Additionally, a
chemical state may consist of an arbitrary number of lower-level chemical states. To avoid confusion,
in the following all chemical states and reactions are

Once the chemical states and reactions that comprise the chemical system under consideration are
defined, their dynamics can be described as a random process that switches between states [114, 115].
For this purpose, the state of the chemical system is described as a vector . Additionally, a set of
reaction channels is established, which cause the transitions between the discrete states of this vector.
Due to the unpredictable nature of the dynamics of the microstates, the time at which a reaction takes
place is not a deterministic quantity. Instead, if the chemical system is in a given state _{α} at time _{γ} can be defined, so that _{γ}d_{γ} is a constant in time. Different chemical states have different
reaction rate constants for their reaction channels. Again due to the assumed instantaneous
equilibration, these reaction rate constants depend only on the current state of the chemical
system irrespective of the previous states of the system. In this case a function can be
defined for every reaction channel that assigns a specific rate to a specific state _{γ} = _{γ}(_{α}).
These functions are called the _{γ}. The thus formulated model
describes a memory-less random process with discrete states, which is usually called a
Markov process [37]. The removal of memory from the system occurs through the thermal
equilibration, which is assumed to happen much faster than the chemical reactions. Examples
for non-markovian behavior are correlated reactions as in the recombination-enhanced
defect reaction (REDR) theory [116] or time-dependent propensity functions as suggested
for reemission after the multi-phonon capture of carriers [117]. In the context of BTI,
non-markovian behavior has been used to explain accelerated recovery observed in charge-pumping
measurements [118]. Interestingly, in standard BTI experiments with constant or AC stress,
even up to 5MHz, no memory-related behavior has been found yet. We thus assume that
the accelerated recovery is a specific effect of the charge-pumping measurements, which
requires deeper experimental investigation, and ignore non-markovian behavior in this
work.

According to the theory of stochastic chemical kinetics [114, 115], the evolution of this system over time can then be described by a chemical master equation

| (2.20) |

where _{0}_{0}.

The master equation approach can be illustrated using the simple example of a system with two
states _{1} and _{2} [37], see Fig. 2.4. The system has two reaction channels 1 and 2, which connect the
two states through the state change vectors _{1} and _{2} as

| (2.21) |

The propensity functions _{1} and _{2} assume the form

_{1}(_{1}) | = _{12} | _{1}(_{2}) | = 0 | (2.22) |

_{2}(_{1}) | = 0 | _{2}(_{2}) | = _{21} | (2.23) |

= _{21}_{2}_{12}_{1} | (2.24) | |

= _{12}_{1}_{21}_{2} | (2.25) |

| (2.26) |

which reduces the master equation of the two-state system to

| (2.27) |

This is the rate-equation of the two-state system, which is equivalent to the master equation for this simple example.

Within the theoretical frame work of the chemical master equation, all the microphysical details
elaborated in the previous section are now contained in the propensity functions _{γ} and the
state-change vectors _{γ} for the Γ reaction channels. These can be either obtained from
experiment, or calculated from microphysical theories, some of which are explained in the
following.

The assumptions introduced above lay the foundation for the development of theoretical methods to calculate rates between chemical states. As mentioned in Sec. 2.4, the description of the nuclear state can either be classical with the full potential energy surface or quantum-mechanical with a highly idealized potential. Both approaches contain specific approximations and thus need to be cautiously applied to real-world situations. Quite generally, for reactions occurring at low temperatures the quantum-mechanical approach will be the most suitable as the effects of quantization and tunnelling may be pronounced, yet the nuclear system will be close to a minimum of the potential energy surface where idealizations are quite accurate. For reactions at higher temperatures, quantum mechanical effects will be less pronounced but the system will occupy states that are further away from the potential minima which are poorly described by the usually harmonic model Hamiltonians.

Barrier hopping transitions such as the structural transitions in the BTI model are reactions which the electronic system follows adiabatically and thus only change the state of the nuclear system. For the temperature ranges relevant to this document these transitions are sufficiently described using classical nuclei. This is a common working-hypothesis in practical calculations [119, 120, 112] and is sometimes supplemented with quantum mechanical corrections [112] where necessary. The microstates consist of the electronic state, which is unchanged except for the adiabatic adjustment of the electronic wave functions, and the position of the nuclear system in phase-space [119], which follows classical mechanics.

As mentioned above, for a given defect model the barrier hopping processes can in principle be described by a propagation of the microstates, i.e. the phase-space position of the nuclei. Calculations of this type are called molecular dynamics simulations [119]. However, the time range that can be treated with molecular dynamics is in the picosecond regime and can only be extended to nanoseconds for classical molecular dynamics potentials, which are not usable in defect calculations. Usual defect reactions in BTI occur in and above the microsecond regime [26, 37], which makes them rare events from a microstate perspective.

The theory employed for the calculation of rates in this picture is the _{α} ^{3N} in the configuration space of the
_{1} to an adjacent region _{2} is assumed to
proceed as a flux _{12}) through the boundary surface _{12}, which is called the

In practice, one is usually only interested in the temperature activation for a certain transition.
Also, the thermal integrals necessary to calculate the flux through the transition state, are quite
tedious to calculate. Due to the assumption of thermal equilibrium, in which the system preferably
occupies states of lower energy, the transition will be dominated by the point of lowest
energy of the transition state. In practical calculations it is thus usually sufficient to find the
energetically optimal path between the energetic minimum point of _{1} and the energetic
minimum point of _{2}, as illustrated in Fig. 2.5. This path is usually termed the

The charge-state transitions in the multi-state multi-phonon transition model for BTI involve the trapping of holes from the silicon valence band. These trapping events are understood as non-radiative multi-phonon (NMP) transitions. Multi-phonon, or vibronic, transitions are transitions involving a change of the electronic and the vibrational state. They play significant roles in many domains of molecular and solid state physics and have been studied extensively from both the experimental and theoretical side [91, 122, 123, 124, 125] especially for transitions at point-defects. The development of the theory of vibronic transitions is tightly linked to the development of the quantum mechanical theory of solids and molecules. As multi-phonon transitions at defects are inherently quantum-mechanical processes that are significantly influenced by a chaotic perturbation (i.e. the heat-bath), their theoretical description is a quite demanding task and usually involves strong idealizations.

The first theory capable of explaining the behavior of molecules on a well-founded quantum-mechanical level was the adiabatic Born-Oppenheimer approximation. This approximation, which is based on a power-series expansion with respect to the ratio of the electronic and nucleic masses, see Sec. 2.2, has some drawbacks that in principle limit its applicability to chemical processes. The first problem is that the potential energy surfaces arising from the power-series expansion are essentially harmonic functions with the equilibrium positions fixed at a certain configuration. In this form, the theory is unable to describe the adiabatic reactions discussed in the previous section. Nevertheless, the comparison with experimental studies showed that “the adiabatic model has a wider application than predicted by this theory” ( [90], Appendix VII). Secondly, within the series expansion formulation there will always be an adiabatic adjustment of the electronic state to the vibrational degrees of freedom, i.e. there will be no transfer of kinetic energy from the nuclei to the electrons and thus no change in the electronic state can be induced from the vibrations of the molecule or solid.

A second approach to the derivation of a molecular wave function was presented by Born more than two decades later [92] (see also [90], Appendix VIII). Instead of a series expansion of the molecular Hamiltonian, this derivation starts with an ansatz-wave function of the form

| (2.28) |

where the _{0} =

Before the ansatz is inserted into the molecular Schrödinger equation (2.1), it is useful to
apply it to the constituent operators of the associated Hamiltonian (2.2) separately. The
potential operators, which are mere functions of the degrees of freedom, are summed up as
=

| (2.29) |

The electronic kinetic energy operator (2.3) acts only on the

| (2.30) |

The application of the electronic Hamiltonian to Ψ(

_{0}Ψ( | = ( | (2.31) |

= _{i} | (2.32) | |

= _{i} | (2.33) | |

= _{i}_{i} | (2.34) |

| (2.35) |

with

| (2.36) |

and the product rules for the gradient and the divergence, it is trivial to show that

| (2.37) |

Applying the Hamilton operator

| (2.38) |

one arrives at

( | |||

_{i}_{A=1}^{N}_{i} | (2.39) |

_{i} | |||

= | (2.40) |

| (2.41) |

with

| (2.42) |

and

= | (2.43) | |

= | (2.44) | |

= | (2.45) | |

_{ij} | = _{A=1}^{N} | (2.46) |

If

| (2.47) |

i.e. the state of the system can be identified by a separate quantum number for the electronic and the vibrational state, respectively. In the more general case where the non-adiabatic part of the Hamiltonian cannot be neglected, the electronic and vibrational quantum numbers are not generally separable and the eigenstates will expand over several electronic states.

There seems to be no established name for the presented approach to the molecular Schrödinger equation. It is sometimes also refered to as Born-Oppenheimer approximation (which is misleading as it is using a completely different derivation), Born-Huang approximation, or Born-ansatz. The non-adiabatic theory not only makes it possible to describe vibronic coupling, but also gives an idea of the conditions under which the adiabatic approximation is reasonable. Following the expressions in [126] and Sec. 3.1.3 of [91], the adiabatic approximation is valid for systems with a large separation between the potential energy surfaces and a weak dependence of the electronic wave functions on the position of the nuclei.

For the derivations in the following sections it is convenient to use a Dirac notation for the Born-Oppenheimer states as

| (2.48) |

Although one may be tempted to view this wave function as a simple tensor product of wave functions
and employ the usual algebra, this is not possible here due to the parametric dependence of the
electronic state on the vibrational coordinate . Nevertheless, this notation is often found in literature
and will also be applied in this document where ever it may serve to improve readability.
However, the reader is advised to handle vectors of the form (2.48) with care. To improve
readability somewhat further this document uses a notation that distinguishes between integrals
over the electronic degrees of freedom

Although in principle a diagonalization of the non-adiabatic Hamiltonian is possible, the resulting
stationary states are not of interest for this work, as we assume transitions between well-defined
electronic states. Indeed, in real-world situations the stationary states corresponding to the
eigenstates of the non-adiabatic Hamiltonian will never be established due to the constant
perturbation from the environment. The influence of the environment is a very complex random
process in nature that cannot be treated, or even formulated, exactly. The theory of vibronic
transitions therefore has to rely on physical intuition and defines the initial and final states of
the transition somewhat heuristically as the initial and final Born-Oppenheimer states
(2.13). It is assumed that initially the electronic subsystem is in a defined state, which is
called i() in the following. The transition rate to a specific final state f() is to be
calculated. The vibrational degrees of freedom assumingly have established an equilibrium
with the surrounding lattice, i.e. the probability _{iI}
as

| (2.49) |

where

| (2.50) |

The Born-Oppenheimer states, however, do not diagonalize the non-adiabatic Hamiltonian and thus
will not give rise to stationary states but will decay over time and thus the system will move between
different electronic states. As the off-diagonals in the non-adiabatic Hamiltonian are usually assumed
to be small compared to the diagonal part, the transitions between the electronic states is treated
using first-order time-dependent perturbation theory [98]. Thus, the rate of transition from Ψ_{iI} to
Ψ_{jJ} reads

| (2.51) |

where the perturbation is the non-adiabacity operator (2.46).

The total rate from i() to f() is obtained by averaging (2.51) over all possible initial states
i

| (2.52) |

where avg _{I} is just a shorthand for an averaging over all states

Now that the basic theory of vibronic transitions is laid out, the next important step is the calculation of the matrix element fFf()i()iI . The integration for this matrix element runs over all electronic and vibrational degrees of freedom and requires the quantum mechanical states of the electronic and the vibrational system to be known. Especially the latter requirement is impossible to fulfill exactly for real-world systems, so additional assumptions have to be made. These assumptions must lead to a model description that captures enough of the basic behavior of the real-world system to give accurate results yet include only so much complexity that a quantum-mechanical treatment is still feasible.

The published models for vibronic transitions [127, 128, 129, 91, 130] contain several approximations that simplify the molecular Hamiltonian (2.2) and the corresponding Born-Oppenheimer states. These approximations are described in the following.

As stated above, the molecular Hamiltonian includes the interactions of all electrons and all nuclei.
However, a vibronic transition will usually leave most of the electrons in their respective state thus it
is not necessary to consider them explicitly in the calculation. What has to be considered is their
interactions with the nuclei and with the explicitly treated electrons (usually only one). The former
interactions together with the Coulomb repulsion

| (2.53) |

where

| (2.54) |

acts only on the nuclei. The details of

Another important approximation concerns the potentials

| (2.55) |

it is always possible to define a coordinate system _{S}_{S}_{S} is obtained from the
diagonalization of the quadratic part as

| (2.56) |

The eigenvalue Λ_{S} is conveniently written as

| (2.57) |

using the mass

| (2.58) |

as in a classical picture this leads to a nuclear motion that is composed of uncoupled harmonic modes
with the frequencies _{S}. The resulting coordinate system is called modal system, _{S} is the modal
mass, and the coordinates _{S} are called modal coordinates. The vibrational Hamiltonian thus assumes
the form

| (2.59) |

An additional assumption that is always employed but rarely stated explicitly is that the
electron-phonon interaction

| (2.60) |

where _{S} is a potential for the electrons which is coupled to the instantaneous nuclear position by
_{S}(_{S}).

The transition rate is now determined from (2.52), with the wave functions in the modal coordinate system

| (2.61) |

and the non-adiabacity operator

_{ij} | = _{S}_{ijS} | (2.62) |

_{ijS} | = | (2.63) |

| (2.64) |

In addition to the

| (2.65) |

An exception is the work of Kubo [128], which uses a highly idealized defect wave function. In the Condon approximation or adiabatic approximation (not to be confused with the adiabatic Born-Oppenheimer approximation), the -dependence of the electronic wave functions is obtained from first-order perturbation theory, i.e.

= _{k≠i} | |||

= _{k≠i} _{S} | (2.66) | ||

_{i}() | = _{i} + | ||

= _{i} + _{S}_{S}_{S}(_{S}) | (2.67) |

= | |||

+ _{S}_{k≠i}_{S}(_{S})_{k′≠j}_{S}^{*}(_{
S}) | |||

+ second order terms | (2.68) |

The perturbation expansion of the energy has an important consequence for the vibrational system, as in the employed Born-Oppenheimer approximation the energy levels act as a potential for the nuclei. The vibrational part of the problem now reads

| (2.69) |

This Hamiltonian is a sum of independent operators for every phonon mode which is solved by the wave function [98]

| (2.70) |

with the corresponding energy

| (2.71) |

and the constituents _{S} and _{iIS} are obtained from

| (2.72) |

In essence, this result makes the true separation of the phonons possible, which is of importance for
the following considerations. The quantum state of the vibrational system _{S}.

In most publications, the coupling function

| (2.73) |

where _{S} is a normalization factor that is sometimes absorbed into the coupling potential. This
approximation is known as the

Inserting the linear relation for _{S}(_{S}) into (2.72) gives

| (2.74) |

with

| (2.75) |

The resulting modal wave functions are harmonic oscillator wave functions [139]

_{S} | = _{λS}^{IS
}(_{S} _{Si}) | (2.76) |

_{S} | = _{S}_{S} | (2.77) |

_{λ}^{n}( | = _{n}(_{λ}^{0}( | (2.78) |

_{λ}^{0}( | = ^{1∕4}^{-} | (2.79) |

In addition to the linear coupling approximation, some works also investigate systems with purely

| (2.80) |

As shown in Fig. 2.7, contrary to the linear coupling, the quadratic coupling does not change the equilibrium coordinate of the mode but the oscillation frequency, leading to a vibrational Schrödinger equation of the form

| (2.81) |

with

| (2.82) |

The vibrational wave functions thus assume the form

_{S} | = _{λSi}^{aS
}(_{S}) | (2.83) |

_{Si} | = _{S}_{Si} | (2.84) |

| (2.85) |

and upon a transition between two electronic states both the equilibrium position and the oscillation frequency of the phonon modes change, yielding the vibrational Schrödinger equation

| (2.86) |

with

_{Si} | = | (2.87) |

_{Si} | = | (2.88) |

_{Si} | = | (2.89) |

Line shapes are a concept originating from the field of optical spectroscopy. As an introduction to the concept, we take a look at the first quantitative theory for the multi-phonon transitions of a point-defect derived by Huang and Rhys [127] in 1952 for the F-Center in potassium bromide. The first part of that work is dedicated to a transition induced by incident radiation. In these transitions, the dominant perturbation arises from the polarization operator and not from the non-adiabacity operator in (2.52). It is safe to assume that acts only on the electrons, again due to the comparably large mass of the nuclei. In the Condon approximation (2.68) the corresponding matrix element reads

+ _{S}_{S}(_{S})_{k≠i} | |||

_{S}^{*}(_{
S})_{k′≠j} | (2.90) |

| (2.91) |

This is known as the quantum mechanical Franck-Condon principle [127, 91]. It is important to stress that the -independence only concerns the matrix element but not its constituents. In this case, the electronic matrix element in (2.51) can be taken out of the integration over , yielding

| (2.92) |

In essence, this describes the transition as an electronic matrix element with vibrational selection
rules. The factors ^{2} arising from the quantum-mechanical Franck-Condon principle are
usually called Franck-Condon factors.

Inserting this matrix element into (2.52) results in the expression [117]

| (2.93) |

where

| (2.94) |

is the

While the Franck-Condon principle is easy to justify for optical matrix elements using the Condon
approximation, the above argument completely fails for the non-adiabacity operator (2.63) which
contains the electronic matrix elements _{S})_{S}) and _{S})_{S}) that do not
induce a direct transition. Using the interaction (2.85), the Condon approximation for these matrix
elements gives with (2.66)

| (2.95) |

and

| (2.96) |

and the transition matrix element reads with (2.63)

| = _{S} _{ijS} | (2.97) |

_{ijS} | = _{ijS}_{IJS} | (2.98) |

_{ijS} | = | (2.99) |

_{IJS} | = _{S} _{S} _{S} _{S} | (2.100) |

Inserting (2.70) for the vibrational wave functions gives

| (2.101) |

This result is interesting as it shows that each mode contributes a Frank-Condon factor to the matrix
element of every other mode. The vibrational part is determined by the change of the potential energy
surface between the states, which is caused by the diagonal part of the electron-phonon
coupling Hamiltonian _{S}_{S}

The target of the present work is to extract relevant parameters for the non-radiative transitions
from atomistic calculations. Although the parameters _{S} and _{S} of the above expression can be
extracted from a suitably designed atomistic model, as will be shown later in this document, this
parameter extraction bears some ambiguity due to the anharmonicity of the potential energy surface
from the electronic structure method. Also, the complicated structure of the investigated defects and
the strong distortion of the host lattice around them makes the assumption of a constant
mode-spectrum quite unlikely. Finally, there is no method available for the estimation of the
cross-terms _{S}_{S′}_{S′} and the conservation of energy. We assume a
perturbation operator

| (2.102) |

This approach is well compatible with published works which either apply the Frank-Condon principle directly to the non-adiabacity operator [127, 117, 131, 142] or where the perturbation does not arise from the non-adiabacity operator but instead from an applied field [143, 144].

As for the barrier hopping transitions in Sec. 2.8, it is also possible to formulate a theory for vibronic transitions that treats the nuclei as classical particles [129, 117]. Several semi-classical approaches have been published over the years, especially for the estimation of reemission probabilities [117] or correlated hopping [116]. In the present work, non-Markovian behavior is not assumed to be relevant, so the theory is again derived using the assumption of thermal equilibrium.

As above, we derive the non-radiative transition from an optical transition by letting the optical energy approach zero. In this approach, the transition between the electronic Born-Oppenheimer wave functions i() and f() proceeds via the optical matrix element . For any position of the nuclei, one can now again employ time-dependent perturbation theory to calculate the transition rate as

| (2.103) |

Again following just the same considerations as before, we assume that the nuclear degrees of freedom are in thermal equilibrium and can be described by classical statistical physics. The probability to find the system in electronic state i at a certain configuration is hence given by

= ^{-1}^{-
} | (2.104) | |

= _{}^{-
} | (2.105) |

= _{ℝN}_{ℝN}^{-1}^{-
} | (2.106) | |

= _{ℝN}^{-1}^{-
}^{2} | (2.107) |

= ^{2} | (2.108) | |

= _{ℝN}^{-1}^{-
} | (2.109) |

Now that the chemical states and reactions are defined and we have obtained rates for the reactions
from the explained microscopic theories, we can calculate the time evolution of the chemical system
from the chemical master equation (2.20). As explained above, this equation is a stochastic differential
equation which assigns a probability to any state vector , given that = _{0} at _{0}. For a system
with a small set of states _{1}_{Ω}

Most chemical systems are constructed of a large number Θ of similar sub-systems with a small number of states Ω, for example a number Θ of well-separated doping atoms within a semiconductor device which can exist in a number Ω of different charge states. The chemical state of the total system can in principle be constructed from all the states of the subsystems, i.e.

| (2.110) |

For every subsystem _{δ1}_{δΓ} and propensity functions _{δ1}(_{δ})_{δΓ}(_{δ}).

In a situation like this a given observable will usually be unable to separate these sub-systems. For
instance for the transport characteristics of a semiconductor device it is irrelevant

| (2.111) |

the chemical state of the total system can be reformulated, considering only the numbers of
subsystems _{ω} in a given state _{ω}. The new total state vector will be

| (2.112) |

the state change vectors will only operate on the numbers _{1}_{Ω}, and the corresponding propensity
functions will take the form

| (2.113) |

corresponding to the probability per unit time that _{2} + O_{2} _{2}O. In those reactions, the corresponding reaction channel changes two or
more numbers of

| (2.114) |

which is the probability per unit time that any two of the _{ω1} specimens of species _{1} and _{ω2}
specimens of species _{2} will undergo the reaction

| (2.115) |

which corresponds to the probability that

Bimolecular reactions usually emerge in situations when at least one of the specimens is able to
move, giving it access to different reaction partners. A mathematical description of the form (2.114) or
(2.115) implies, that the reaction probability is the same for every pair of reactants, i.e. a
homogeneous distribution of the reacting specimens within the volume associated with the chemical
system. This condition is commonly referred to as the system being

As stated above, it is quite common that observables only determine how many constituents of a
chemical system have assumed a certain state, or equivalently how many specimens of a certain species
exist in that system. Usually the most interesting quantities for a chemical system in number
representation is thus the development of the

| (2.116) |

which corresponds to the first moment of the probability distribution

To obtain these averages from the SSA one has to take the average of several simulation runs, which is usually a quite time-consuming task. Especially if there are very fast reactions or a large number of particles in the system, the simulation of longer time-ranges quickly becomes unfeasible using this approach. A computationally more efficient approach is to directly calculate the dynamics of the average values. For convenience we introduce the average number vector = . The time evolution of is calculated from

| (2.117) |

The substitution of the right side of (2.20) for the derivative of

| (2.118) |

For the case of unimolecular reactions (2.113), whose propensity functions are linear in , it is trivial
to show that = _{γ}

| (2.119) |

Due to their efficiency and simplicity, reaction rate equations are a popular method to describe reactive systems. However, one has to be aware that this description is only exact if all the reactions in the system are unimolecular. Nevertheless, reaction rate equations are also commonly used for systems with bimolecular reactions and even reactive-diffusive systems, as in the reaction-diffusion model for NBTI. The consequences and inadequacies of this approach are discussed in Chap. 3.

As a final remark, it is common not to use the number-representation of . Instead, the number of
particles is usually normalized to the volume of the system, leading to a formulation based on