2  Theory of States and Reactions

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.1  Problem Statement: The Molecular Hamiltonian

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, ⃗R and ⃗r denote the position vector of all nuclei and all electrons, respectively. For a system containing n electrons and N nuclei, ⃗r will be 3n dimensional and ⃗R will have 3N dimensions. The positions of the individual particles will be denoted ⃗ra and ⃗RA and the associated masses are me for the electrons (which is a fundamental constant) and MA 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 [889091]

ˆH Ψ(⃗r, ⃗R ) = E Ψ(⃗r,R⃗),

where Ψ(⃗r,⃗
R) is called the vibronic wave function as it expands over the electronic and the nuclear (vibrational) degrees of freedom, and

Hˆ = ˆTe + ˆTN + ˆVee + ˆVNN + VˆeN

is the Hamiltonian of the system, which is sometimes also termed the molecular Hamiltonian. It consists of

Tˆe = - a=1n ℏ2
2me⃗∇a kinetic energy operator for electrons (2.3)
TˆN = - A=1N  2
2MA⃗∇A kinetic energy operator for nuclei (2.4)
Vee = -q0--
4πϵ0 a<b---1----
|⃗ra - ⃗rb| electron-electron repulsion (2.5)
VˆNN =  q
4πϵ0 A<B  Z  Z
|R⃗A  - ⃗RB | nucleus-nucleus repulsion (2.6)
VeN = --q0--
4π ϵ0 a=1N A=1M---ZA----
|⃗r - ⃗R  |
  a    A electron-nucleus attraction. (2.7)
Despite the striking elegance of this equation, it is unsolvable in its pure form for all but the simplest molecules. However, a profound simplification can be drawn from a separation of the total Schrödinger equation into separate problems for the electronic and the nuclear degrees of freedom. The first theoretically well-founded approximation of this type was put forward by Max Born and Robert Oppenheimer in 1927.

2.2  The Born-Oppenheimer Approximation

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

ˆH = Hˆ0 +  κ4 ˆH1


κ = ∘ ------
4 me∕MÆ (2.9)
Ĥ0 = ˆTe + Vˆee + VˆeN + VˆNN (2.10)
Ĥ1 = - ℏ2
2me A=1N⃗∇ A, (2.11)
where M is the average ionic mass. The solution of the molecular Schrödinger equation is then expanded in powers of κ  [8890]. The zero-order term of this expansion became known as the ‘clamped nuclei’ Hamiltonian. In this case, the problem reduces to the Schrödinger equation for the electrons with the nuclei fixed in their respective positions
Hˆ0 φi(⃗r;R ⃗) = Ei(⃗R)φi(⃗r; ⃗R).

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 Ei(⃗
R) associated with the electronic state φi(⃗r;⃗
R) will act as a potential energy on the nuclei. Thus, the molecular wave function is obtained as

ΨiI(⃗r, ⃗R) = φi(⃗r;R ⃗)ηiI(R⃗),

where φi(⃗r;⃗R) is one of the solutions of (2.12) and ηiI(⃗R) is determined from

(ˆTN + Ei(⃗R))ηiI(R⃗) = EiIηiI(⃗R ).

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 adiabatically, i.e. without any transfer of momentum from or to the nuclei. This approximation is consequently called the adiabatic Born-Oppenheimer approximation. The potential Ei(R⃗) is the potential energy surface of the molecule which is a central concept of chemistry and plays a significant role in the theory of chemical states and reactions. It is important to realize that different electronic states will give rise to different potential energy surfaces, as illustrated in Fig. 2.1. This difference between the potentials is at the heart of the electrochemical reactions discussed in this document. Especial importance is assigned to the potential energy surface E0(⃗
R) of the lowest electronic state φ0(⃗r; ⃗
R), which is the dominant state in thermal equilibrium.

The original derivation of Born and Oppenheimer  [8890] 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 [9290]. 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.

Figure 2.1: Illustration of the dependence of different electronic states on the nuclear configuration. The adiabatically adjusting energies of the electronic states act as potential for the nuclei. Different electronic states may give rise to drastically different potentials. For example, while the potentials E0 and E1 in the figure show only a shift of the optimum position and a change of the curvature, E2 does not show any optimum position in the shown configuration range. In a molecule, E0 and E1 correspond to bonding states, while E2 leads to an antibonding (repulsive) behavior.

SVG-Viewer needed.

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.

2.3  Electronic Structure Methods

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 electronic structure methods or quantum chemistry methods. The term ‘quantum chemistry method’ may sometimes be a source of confusion, as especially in the materials science community this term is used for wave function based methods to distinguish them from density functional methods. However, in this document the term quantum chemistry method and electronic structure method are both used interchangeably for methods that calculate the potential energy surface based on a solution of the electronic Schrödinger equation.

2.3.1  Self-Consistent Field Methods: Hartree-Fock

SVG-Viewer needed.

Figure 2.2: An intuitive picture of the self-consistent field approach. Instead of considering all the interactions between all particles (left), the problem is reformulated as a system of isolated particles moving in a mean field (right).

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  [939495]. 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 V 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).

HˆH  = ˆTe + ˆVeN + VˆH (φ)

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.

φ(⃗r ,⃗r ,...) = φ (⃗r )φ  (⃗r) ...
   1  2        1  1  2 2

A central problem of this ansatz is that it does not have the required fermionic symmetries upon particle exchange [969798]. 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.

             |                    |
             || φ1(⃗r1)  φ2(⃗r1) ... ||
φ(⃗r1,⃗r2,...) = || φ1(⃗r2) φ2(⃗r2) ... ||
             ||    ..      ..    ... ||
                  .      .

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 ˆVF.

ˆHHF  = ˆTe + VˆeN + (ˆVH(φ )- ˆVF(φ))

The resulting method was named Hartree-Fock self-consistent field method, the new term in the equation is usually called exchange- or Fock-operator1 . 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 [9697]. 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  [9997].

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 (O(N4), where N is the number of basis functions). The sharp rise in computational demand for increasing number of basis functions results in the need for highly efficient basis sets for practical Hartree-Fock-type calculations. The employed basis sets are usually linear combinations of atomic orbitals (LCAO), most typically constructed using Gaussian functions for computational reasons  [9996]. These basis sets themselves require an appreciable amount of fine-tuning and are not to be considered general-purpose tools. For the study of defects, Hartree-Fock methods are sometimes used in embedded  [7766] or isolated  [5457] cluster models.

2.3.2  Self-Consistent Field Methods: Density Functional Theory

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

  1. for a given external potential Vext which acts equally on all of the electrons, there is a ground state electron density n0. This density is unique to the external potential apart from a trivial constant energy.
  2. there is a functional E[n] of the electron density which, when minimized in n for a given Vext, gives the ground state energy E0 and the ground state electron density n0 of the many-body system.

From these two findings it follows that the exact ground state energy of the (3N-dimensional) electronic system, i.e. the ground state potential energy surface, could be found by minimizing the Hohenberg-Kohn functional of the (three-dimensional) electron density, if this functional was known exactly [10096]. However, as the exact functional itself depends on the wave functions of the many-body system [96], additional approximations are necessary. To make the problem solvable, Kohn and Sham made the assumption that a non-interacting auxiliary system can be defined

HˆKS  = Tˆe + ˆVeN + ˆVH + ˆVxc

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 ˆVxc that itself depends on the ground state electron density. Due to its favorable scaling properties (O(N3)) 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 Hybrid-functionals. These functionals combine the original (semi-) local density functional of standard DFT approaches with a fraction of Hartree-Fock exchange. Hybrid functionals significantly improve the description of semiconductors and insulators, at the expense of increased computational complexity.

2.3.3  Explicit Many-Body Methods

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 [9697]. 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.

2.3.4  Empirical Methods

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  [103104]. A comparison of semi-empirical methods for the calculation of the oxygen vacancy defect in α-quartz is given in  [105], which shows a strong dependence of the predicted defect levels on the employed method.

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  [106107108], but also for defects in SiO2  [55]. Compared to semiempirical methods, the tight-binding Hamiltonians rely on even stronger simplifications which makes the parametrizations of this method even less transferable. Thus, when employing tight-binding methods one cannot rely on standard parametrizations as is usually done in semiempirical methods.

A review of different empirical molecular dynamics potentials applied to silicon or SiO2 systems is found in  [109]. Although empirical potentials are computationally very efficient and make calculations with millions of atoms possible  [78], they have too little predictive power for defect calculations as those involve unpaired spins and charged states which are difficult to include in this formalism in a physically meaningful way. However, empirical potentials are used extensively for the generation of amorphous structures that are subsequently used for defect calculations  [11078111].

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 SiO2, the calibration of all discussed empirical methods requires an extensive training set. Unfortunately there is no reliable experimental data available that can be used as a reference for parametrization, so for the present work we are bound to ab-inito methods.

2.4  The Nuclear Problem: Chemical Microstates

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 1and 2 2. For situations where the electronic state changes non-adiabatically, as in the charge transition reactions 1 2and 2 1, or at very low temperatures the quantum-mechanical nature of the nuclei must be considered and the harmonic approximation, or a similar simplification, is necessary [91].

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 ηiI(⃗R) in a quantum-mechanical description and the instantaneous position in phase-space (⃗R,⃗P) — ⃗P denoting the momenta of all particles — in a classical description. In the following, these combinations of electronic and nuclear states will be called microstates in accord with the theory of statistical physics [112]. As will be explained in more detail in Sec. 2.9.1, the adiabatic Born-Oppenheimer approximation can be assumed to be valid for all regions of the nuclear configuration space in which the potential energy surfaces are sufficiently separated. The systems of interest for the present document spend most of their time inside these regions and the electronic transitions happen either through radiative transitions or non-radiatively in the vicinity of potential energy surface crossings. Both of these transitions are sufficiently described as occurring instantaneously. Consequently, the transitions between the microstates are governed by the dynamics of the nuclei and the transitions between different potential energy surfaces. Although in principle the dynamics of any molecule or solid can be described in this way, in practice it is usually tedious to describe the long-term evolution of even the simplest chemical system from the trajectories of its nuclei and transitions between electronic states. Fortunately, the time-evolution of chemical systems usually consists of longer time-periods of thermal quasi-equilibrium which are interrupted by transitional motions of very short duration. This behavior is exploited and formalized in the theory of chemical states and reactions.

2.5  Chemical States and Reactions

Figure 2.3: The potential energy surface of an atomic system typically has several minima which are separated by energetic barriers. The nuclei will spend a significant amount of time oscillating randomly around those minima and will eventually pass the barriers at time scales that are too short to be observable. Along the barriers, the configuration space can be separated into chemical states, illustrated as grey/white region. The potential energy surface could for example arise from two different bonding states of a defect, as schematically shown in the boxes.

SVG-Viewer needed.

In order to understand chemical kinetics, it is first necessary to define in a physically meaningful way the chemical states, which are the central actors of chemistry, and the chemical reactions, which change the system under consideration from one chemical state to another. Quite generally, chemical states are groups of microstates. The microstates that belong to a certain chemical state could in principle be selected completely arbitrarily, but usually a chemical state is defined and receives its meaning from a certain observable2 . 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 elementary, i.e. they are not further separable into long-term stable chemical states.

2.6  The Chemical Master Equation

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 [114115]. For this purpose, the state of the chemical system is described as a vector ⃗x. 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 ⃗xα at time t, for every reaction channel γ a reaction rate cγ can be defined, so that cγdt is the probability of the reaction taking place between t and t + dt [115]. Due to the assumed instantaneous equilibration of the microstates, cγ 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 cγ = aγ(⃗xα). These functions are called the propensity functions  [115]. The change induced by the reaction channel γ is described using the state change vector ⃗νγ. 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  [114115], the evolution of this system over time can then be described by a chemical master equation

∂P(⃗x,t)-  ∑Γ
  ∂t    =    [aγ(⃗x - ⃗νγ)P (⃗x- ⃗νγ,t)- aγ(⃗x)P (⃗x,t)],

where P(⃗x,t) = P(⃗X = ⃗x,t|⃗x0,t0) is the probability that the stochastic process ⃗X(t) equals ⃗x at time t, given that X⃗(t0) = ⃗x0.

Figure 2.4: Illustration of the example system discussed in the text. The system exists in one of the two states ⃗x1 and ⃗x2. Transitions between these states are caused by two reaction channels with the state-change vectors ⃗ν1 and ⃗ν2, and the associated rates k12 and k21.

SVG-Viewer needed.

The master equation approach can be illustrated using the simple example of a system with two states ⃗x1 and ⃗x2  [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

⃗x1 + ⃗ν1 = ⃗x2 and   ⃗x2 + ⃗ν2 = ⃗x1

The propensity functions a1 and a2 assume the form

a1(⃗x1) = k12, a1(       ⃗x2) = 0, (2.22)
a2(⃗x1) = 0,and a2(       ⃗x2) = k21. (2.23)
The master equation for this system consequently reads
∂P (⃗x1,t)
---∂t---- = k21P(⃗x2,t) - k12P(⃗x1,t) (2.24)
∂P (⃗x  ,t)
   ∂t = k12P(⃗x1,t) - k21P(⃗x2,t) (2.25)
As the system can only exist in one of the two states at a time, it follows that
P (⃗x1,t) = 1 - P(⃗x2,t) = p(t),

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

∂p(t)=  k (1 - p(t))-  k p(t).
 ∂t      21             12

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 aγ 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.

2.7  The Calculation of Rates

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.

2.8  Barrier Hopping Transitions

Figure 2.5: Illustrative example of a potential energy surface along a reaction path. The reaction path itself is the optimal path on the 3N-dimensional potential energy surface between two minima. The point of maximum energy, which is a saddle-point of the potential energy surface, is sometimes called the transition state. Its energy with respect to the energetic minimum of the initial state is the activation energy Ea. The transition point is also the border which separates the initial and the final state in a classical picture.

SVG-Viewer needed.

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  [119120112] 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  [2637], which makes them rare events from a microstate perspective.

The theory employed for the calculation of rates in this picture is the transition state theory. This theory is the mathematical formulation of the assumptions of Sec. 2.5. As in Fig. 2.3, the chemical states in the transition state theory are defined by regions Rα 3N in the configuration space of the N nuclei. These regions are non-overlapping and defined by the barriers of the potential energy surface. The transition from a region R1 to an adjacent region R2 is assumed to proceed as a flux F(R12) through the boundary surface R12, which is called the transition state [121].

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 R1 and the energetic minimum point of R2, as illustrated in Fig. 2.5. This path is usually termed the reaction path. The maximum energy along this path is called the activation energy. Macroscopic theories usually model the temperature activation of hopping transitions using an Arrhenius law that has the activation energy as parameter. However, one has to be aware that this neglects entropic properties which are only sufficiently described in the full transition state theory.

2.9  Multi-Phonon Transitions

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 [91122123124125] 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.

2.9.1  Vibronic Coupling

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

Ψ (⃗r, ⃗R) =   φi(⃗r; ⃗R)ηi(⃗R),

where the φi(⃗r;⃗R) are the solutions of the electronic Hamiltonian Ĥ0 = ˆTe + ˆVee + ˆVNN + VˆeN in (2.12), which are assumed to be known, and ηi(⃗R) can be seen as a weighting-factor for the electronic wave function i. In a way, this ansatz expands the molecular wave function using the adiabatic Born-Oppenheimer wave functions as a basis set. This expansion considers the correlation of the electronic motion to the motion of the nuclei, but neglects the effects on the nuclei induced by the instantaneous positions of the electrons. Thus, it is not exact although stated otherwise sometimes.

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 Vˆ = ˆVee + VˆNN + VˆeN = V (⃗r,⃗R) and are applied as

ˆ     ⃗    ∑       ⃗      ⃗    ⃗
V Ψ(⃗r,R ) =    V(⃗r,R )φi(⃗r;R )ηi(R ).

The electronic kinetic energy operator (2.3) acts only on the φi(⃗r;⃗R)

 ˆ     ⃗    ∑     ⃗ ˆ      ⃗
TeΨ (⃗r,R) =    ηi(R)Teφi(⃗r;R).

The application of the electronic Hamiltonian to Ψ(⃗r,⃗R) thus gives

R) = ( ˆ
Te +  ˆ
R) (2.31)
= iηi(⃗R)ˆTeφi(⃗r;⃗R) + V (⃗r,⃗R)φi(⃗r;⃗R)ηi(⃗R) (2.32)
= iηi(⃗R)(                          )
 Tˆeφi(⃗r; ⃗R )+ V (⃗r,R⃗)φi(⃗r;R ⃗) (2.33)
= iηi(⃗R)(          )
 Hˆ0 φi(⃗r; ⃗R) = iηi(R⃗)(             )
 Ei (R⃗)φi(⃗r; ⃗R) (2.34)
where Ei(R⃗) is again the eigenvalue for φi(⃗r;⃗R) with respect to Ĥ0. The nuclear kinetic energy operator (2.4) acts on both the φi(⃗r;⃗R) and the ηi(⃗R). Using the definition
     ∑N  PˆA ⋅PˆA
ˆTN =     -2M----
     A=1     A


ˆPA = - iℏ ⃗∇A,

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

ˆ      ∑     ˆ       ˆ      ∑N  -1--(ˆ   )  ( ˆ   )
TNΨ  =    φiTN ηi + ηiTNφi +    MA   PAφi  ⋅ PA ηi .
        i                   A=1

Applying the Hamilton operator

ˆH Ψ(⃗r, ⃗R ) = (Hˆ0 + TˆN )Ψ(⃗r, ⃗R ),

one arrives at

ĤΨ (⃗r,⃗R) =
iηiEiφi + φiˆTNηi + ηiˆTNφi + A=1N 1
MA(     )
  ˆPAφi(     )
  ˆPAηi = E iφiηi. (2.39)
Finally, the equation is expanded in the electronic basis set. Therefore we multiply with φj* from the left and integrate over the electronic degrees of freedom ⃗r, yielding the system of equations
(           )
 TˆN + Ej (⃗R) ηj(R⃗)+
i(                                                           )
  ∫                       ∑N   1  ∫
     φ*j(⃗r; ⃗R)ˆTN φi(⃗r;R⃗)d⃗r +   ----   φ*j(⃗r; ⃗R)ˆPA φi(⃗r;R⃗)d⃗r ⋅PˆA
                          A=1 MAηi(R⃗)
= j(⃗R) (2.40)
for i running over all eigenstates of Ĥ0. This system of equations is conveniently rewritten into matrix form as
H⃗η(R⃗) = E ⃗η(⃗R )


       (        )
       ||  η2(⃗R) ||
⃗η(⃗R ) = |( η3(⃗R) |)


H = HA + HNA (2.43)
HA = (  ˆTN + E1(R⃗)      0            0       ⋅⋅⋅)
|       0       ˆT  + E (⃗R)       0       ⋅⋅⋅|
||                N    2      ˆ       ⃗      ||
(       0           0       TN + E3 (R)  ⋅⋅⋅)
        ...            ...           ...       ... (2.44)
HNA = (                    )
   ˆL11  ˆL12  ˆL13  ⋅⋅⋅
||  ˆL21  ˆL22  ˆL23  ⋅⋅⋅||
|(  ˆL31  ˆL32  ˆL33  ⋅⋅⋅|)
    ..    ..    ..   ..
    .    .    .     . (2.45)
Lij = φ* j(⃗r; ⃗
R)d⃗r + A=1N-1--
MA φ* j(⃗r;⃗
R)d⃗r  ˆ
PA. (2.46)
HA is now again the adiabatic Hamiltonian known from the Born-Oppenheimer approximation, the transfer of kinetic energy from the vibrational to the electronic system is described by the non-adiabatic part of the Hamiltonian HNA. As one would expect, the adiabatic Hamiltonian is diagonal in the matrix representation and the coupling between the potential-energy surfaces Ei(⃗R) is solely due to HNA, more specifically the off-diagonal matrix elements ˆLij.

If HNA is neglected, the original Born-Oppenheimer approximation is retained and the resulting eigenstates will be of the form (2.13)

ΨiI(⃗r, ⃗R) = φi(⃗r;R ⃗)ηiI(R⃗),

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

|ΨiI⟩ = |i(⃗R )⟩|iI⟩⟩.

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 ⃗
R. 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 ⟨⋅|⋅⟩ and over the nuclear degrees of freedom ⟨⟨⋅|⋅⟩ ⟩.

2.9.2  Quantum Mechanical Theory of Vibronic Transitions

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(⃗
R)⟩ in the following. The transition rate to a specific final state |f( ⃗
R)⟩ is to be calculated. The vibrational degrees of freedom assumingly have established an equilibrium with the surrounding lattice, i.e. the probability p of finding the vibrational subsystem in a state |iI⟩ ⟩, given that the electronic system is in state i depends on the energy EiI as

p  =  Z-1e- kiBIT ,

where Z is the canonical partition function

    ∑    - EiI
Z =     e kBT .

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

             |          |
WiI→jJ =  2π-||⟨ΨiI|ˆL|ΨjJ⟩||2δ(EjJ - EiI),

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

The total rate from |i(⃗
R)⟩ to |f(⃗
R)⟩ is obtained by averaging (2.51) over all possible initial states |iI⟩ ⟩ and summing over all possible final states |fF⟩ ⟩

        2π    ∑   ||                   ||2
Wi →f = ---avg    |⟨⟨fF |⟨f(⃗R)|ˆL|i(⃗R )⟩|iI ⟩⟩| δ(EfF - EiI).
         ℏ  I  F

where avg I is just a shorthand for an averaging over all states I, using (2.49) as the weighting factors.

2.9.3  Model Matrix Elements

Now that the basic theory of vibronic transitions is laid out, the next important step is the calculation of the matrix element ⟨ ⟨fF|⟨f(⃗
R)⟩|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  [12712812991130] contain several approximations that simplify the molecular Hamiltonian (2.2) and the corresponding Born-Oppenheimer states. These approximations are described in the following.

Unaffected Electrons

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 ˆVNN will give rise to a potential for the nuclei, the latter interactions are considered in the effective electronic Hamiltonian Ĥe. The resulting Hamiltonian is

 ˆ   ˆ    ˆ     ˆ
H  = He + HN + HeN,

where Ĥe only acts on the (considered) electronic degrees of freedom and

HˆN  = ˆTN + ˆVN

acts only on the nuclei. The details of Ĥe and ĤeN depend heavily on the system under consideration and cannot be given in a general form.

Phonon Modes

Another important approximation concerns the potentials ˆVN and ĤeN. It is usually assumed that they are smooth enough so that a Taylor series expansion is suitable and this expansion is usually canceled after the second order for  ˆ
VN. For a parabolic potential

  ⃗         ⃗   ⃗   ⃗T ˜ ⃗
V(R ) = V0 + V1 ⋅R + R V2R

it is always possible to define a coordinate system R⃗which is shifted so that the linear part of the potential is zero. On top of this, a coordinate system ⃗Q = SQS⃗eS is obtained from the diagonalization of the quadratic part as

˜V2⃗eS = ΛS⃗eS.

The eigenvalue ΛS is conveniently written as

ΛS =    2

using the mass

MS  =    (e2SAx + e2SAy + e2SAz)MA,

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, MS is the modal mass, and the coordinates QS are called modal coordinates. The vibrational Hamiltonian thus assumes the form

      ∑     ℏ2   ∂2    M  ω2
ˆHN  =    - --------2+  --S-S-Q2S.
       S   2MS  ∂Q S     2

An additional assumption that is always employed but rarely stated explicitly is that the electron-phonon interaction ĤeN does not change the modal structure, i.e. it is also diagonal in  ⃗
Q [129131]. In this case the interaction Hamiltonian can be written as

 ˆ     ∑  ˆ
HeN =     USFS (QS ),

where ŪS is a potential for the electrons which is coupled to the instantaneous nuclear position by FS(QS).

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

      ⃗         ⃗     ⃗       ⃗
ΨiI(⃗x,Q ) = φi(⃗x;Q )ηiI(Q ) = |i(Q )⟩|iI⟩⟩

and the non-adiabacity operator

ˆLij = SLˆijS (2.62)
LijS = --ℏ2--
MS⟨j( ⃗
∂QS. (2.63)
The |i(⃗Q)⟩ are the solutions of the equation
( ˆHe + ˆHeN )|i(⃗Q)⟩ = Ei(⃗Q)|i(⃗Q )⟩,

Perturbation Theory

In addition to the time dependent perturbation theory employed to calculate the transition rates, nearly all published model matrix elements treat the dependence of the electronic wave functions and energies on ⃗
Q using a perturbation expansion based on the ( ⃗
Q-independent) solutions of the electronic Hamiltonian

ˆHe|i⟩ = Ei|i⟩.

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 Q⃗-dependence of the electronic wave functions is obtained from first-order perturbation theory, i.e.

|i(Q⃗)⟩ = |i⟩ + ki⟨k| ˆH |i⟩
 Ei - Ek|k⟩
= |i⟩ + ki S⟨k|UˆS-|i⟩FS-(QS-)-
   Ei - Ek|k⟩ (2.66)
Q) = Ei + ⟨i|ĤeN|i ⟩
= Ei + S⟨i|ŪS|i⟩FS(QS). (2.67)
This approximation was used initially by Huang and Rhys [127132] and was subsequently employed in many publications. The matrix element for a general operator Ō is
⟨j(⃗Q)|Ō|i(⃗Q)⟩ = ⟨j|Ō|i⟩
+ S( ki⟨j|ŌFS(QS)|k⟩⟨k|--UˆS----
Ei - Ek|i ⟩ + kj⟨j|---ˆUS---
Ej - Ek|k⟩⟨k|FS*(Q S)Ō|i⟩)
+ second order terms. (2.68)
Unfortunately, the matrix elements for non-radiative transition rates calculated using the Condon approximation were too small by orders of magnitude  [132133] and so different approximations were developed. The most popular approaches are the static coupling scheme  [13413589132] and the non-Condon approximations  [130], which were based on an “infinite order” perturbation expansion. Both methods give transition rates closer to experiment and the relative advantages of each of the approximations have been discussed for more than two decades. In the early 1980s it was conclusively shown by several authors that in principle all three approximations are equivalent and that significant differences arise from the employed wave functions, which makes this part of the theory slightly arbitrary [132136137138]. In the remaining part of this document the standard Condon approximation is applied for its simplicity. The electronic matrix elements employed contain parameters that are calibrated to experimental data, which is assumed to properly account for the shortcomings of the approximation.

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   ∂2    M  ω2
( ˆHN + Ei(⃗Q ))|iI⟩⟩ = Ei +    - --------2-+ --S-S-Q2S + ⟨i|UˆS |i⟩FS (QS )  |iI⟩⟩ = EiI|iI ⟩⟩.
                          S   2MS  ∂Q S     2

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

|iI⟩⟩ = |iI ⟩⟩ ⊗ |iI ⟩⟩⊗ |iI ⟩⟩⊗ ⋅⋅⋅|iI ⟩⟩⋅⋅⋅
        1      2      3        S

with the corresponding energy

EiI = Ei +    EiIS

and the constituents |iIS⟩ ⟩ and EiIS are obtained from

(     2   2         2                   )
 - -ℏ----∂---+ MS-ω-SQ2 + ⟨i|UˆS |i⟩FS (QS )  |iIS⟩⟩ = EiI |iIS⟩⟩.
   2MS  ∂Q2S     2     S                            S

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 I that belongs to the electronic state i now consists of the independent states of all phonon modes IS.


In most publications, the coupling function F(⃗
Q) is a linear expression in the modal coordinates [130]

ˆHeN =     ˆUSaSQS,

where aS is a normalization factor that is sometimes absorbed into the coupling potential. This approximation is known as the linear coupling approximation and is widely used throughout literature  [127130132].

Inserting the linear relation for FS(QS) into (2.72) gives

(     2   2         2                 )
  - -ℏ---∂---+ MS--ωS(QS + QSi )2 - E ′ |iIS ⟩⟩ = EiI |iIS⟩⟩
    2MS ∂Q2S      2                  Si           S


       ⟨i|ˆUS|i⟩aS       ′    ⟨i|ˆUS |i⟩2a2S
QSi =  -M---ω2---and ESi = --2M--ω2--.
          S  S                  S S

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

|iIS⟩ ⟩ = hλSIS (QS - QSi) (2.76)
λS = MSωS (2.77)
hλn(Q) =   1
  2 n!Hn(Q)hλ0(Q) (2.78)
hλ0(Q) = (    )
  πℏ14e-Q22λℏ, (2.79)
where Hn(Q) are the Hermite polynomials [139]. Each mode S has an associated vibrational center QSi that is determined by the electronic state i. Upon a transition to a different electronic state, the vibrational wave functions are shifted to the new potential minimum as illustrated in Fig. 2.6.


Figure 2.6: Potential energy surfaces and the corresponding vibrational wave functions for one mode in the linear coupling approximation. (left) Different electronic states give rise to harmonic potential energy surfaces that differ only in the energetic and spatial position of the minimum. The corresponding set of vibrational wave functions is shifted accordingly but remains unchanged otherwise. (right) The vibronic transition theory (2.51) describes the transition between two specific vibrational sub-states as indicated by the arrow. The total rate is computed by thermally averaging over all possible transitions.


Figure 2.7: Potential energy surfaces and vibrational wave functions for a transition in the quadratic coupling (left) and in the combined linear and quadratic coupling approximation (right). In the purely quadratic coupling approximation, there is no shift in the equilibrium position but a change in the oscillation frequency. In the combined version, both the equilibrium position and the oscillation frequency change.

In addition to the linear coupling approximation, some works also investigate systems with purely quadratic coupling which has the form  [128129]

HˆeN =     ˆUSbSQ2S.

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  ∂2    MS ω2    )
  ---------2-+ -----SiQ2S  |iIS⟩⟩ = EiIS|iIS⟩⟩
   2MS  ∂Q S      2


ω2Si = ω2S + 2bS-⟨i|ˆUS|i⟩.

The vibrational wave functions thus assume the form

|iaS⟩ ⟩ = hλSiaS (QS) (2.83)
λSi = MSωSi. (2.84)
Both the linear coupling approximation and the quadratic approximation have been studied in numerous papers. As is shown in Chap. 4 and also described in our publications  [86140], the defects that are considered relevant for MOS reliability need a combination of both couplings to even capture the basic features of the potential energy surfaces calculated from density functional theory. Thus we are dealing with an interaction of the form
ˆHeN =     ˆUS(aSQS  + bSQ2 )
       S                S

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   2         2                  )
 - -ℏ----∂---+ MS-ω-Si(Q   - Q  )2 + E ′ |iI ⟩⟩ = E   |iI ⟩⟩
   2MS  ∂Q2S      2     S    Si      Si    S     iIS   S


QSi =      aS
2bS + ⟨i|ˆU-S|i⟩
         S, (2.87)
ESi = -     2       2
4bS⟨i|ˆUS |i⟩+ 2MS ωS, and (2.88)
ωSi = √ --
  2∘ ----------------
  2bS⟨i|UˆS-|i⟩    2
     MS      + ωS (2.89)

2.9.4  The Line Shape Function

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 Mˆ and not from the non-adiabacity operator ˆ
L in (2.52). It is safe to assume that ˆ
M 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

Q)| ˆ
Q)⟩ = ⟨j|ˆ
+ S(FS(QS) ki⟨j|ˆ
Ei - Ek|i⟩
    + FS*(Q S) kj⟨j|   ˆ
Ej - Ek|k⟩⟨k|M ˆ|i⟩) (2.90)
where the first term represents the direct transition between the two electronic states due to the incident radiation, and the other two terms account for transitions mediated by electron-phonon interactions. It is commonly assumed that the indirect transitions are negligible compared to the direct matrix element ⟨j| ˆM|i⟩. Consequently the sums in (2.90) can be neglected and the electronic matrix element is approximately independent of the position of the nuclei.
⟨j(⃗Q )|M ˆ|i(⃗Q)⟩ ≈ ⟨j|Mˆ|i⟩

This is known as the quantum mechanical Franck-Condon principle [12791]. It is important to stress that the Q⃗-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 ⃗
Q, yielding

⟨ΨjJ|Mˆ|ΨiI⟩ = ⟨⟨jJ |⟨j(⃗Q )|Mˆ|i(Q⃗)⟩|iI⟩⟩ ≈ ⟨j|M ˆ|i⟩⟨⟨jJ|iI⟩⟩.

In essence, this describes the transition as an electronic matrix element with vibrational selection rules. The factors |⟨⟨jJ|iI⟩⟩|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π-||   ˆ  ||2
Wi →f(hν) =  ℏ |⟨f|M |i⟩| f(hν),

where is the energy of the radiation and

           ∑          2
f(hν) = avg    |⟨⟨fF |iI⟩⟩| δ(EfF - EiI - hν)
         I  F

is the line shape function. In optical spectroscopy line shape functions describe the thermal broadening of absorption peaks [12791117], but the concept is also applicable to non-radiative carrier capture at defects in semiconductors, as shown below.

2.9.5  Line Shapes for Non-Radiative Transitions

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 ⟨j(QS)|  ∂
∂QS-|i(QS)⟩ and ⟨j(QS)|∂2
∂Q2S- |i(QS)⟩ that do not induce a direct transition. Using the interaction (2.85), the Condon approximation for these matrix elements gives with (2.66)

       ∂                           ˆUS
⟨j(Q⃗)|∂Q---|i(⃗Q )⟩ = (aS + 2bSQS )⟨j|E----E-|i⟩
        S                         i    j


        2                  ˆ
⟨j(⃗Q)|-∂-2|i(Q⃗)⟩ = 2bS⟨j|-US---|i⟩,
      ∂Q S              Ei - Ej

and the transition matrix element reads with (2.63)

⟨ ⟨jJ|⟨j(⃗Q)|ˆL|i(⃗Q)⟩|iI⟩ ⟩ = S⟨ ⟨jJ|ˆLijS|iI⟩ ⟩ (2.97)
⟨ ⟨jJ|ˆLijS|iI⟩ ⟩ = -UijSSIJS (2.98)
UijS =  ℏ2
MS⟨j|  Uˆ
Ei - Ef|i⟩ (2.99)
SIJS = bS⟨ ⟨jJ|iI⟩ ⟩ + aS⟨ ⟨jJ|--∂--
∂QS|iI⟩ ⟩ + 2bS⟨ ⟨jJ|QS-∂---
∂QS|iI⟩ ⟩. (2.100)
The matrix elements of the separate modal non-adiabacity operators again resemble the radiative matrix elements obtained from the Franck-Condon principle above as the expression consists of an electronic part and a vibrational part. Unlike the Frank-Condon expression above, however, this vibrational part contains both a dependence on the vibrational coordinates themselves as well as a dependence on their derivative.

Inserting (2.70) for the vibrational wave functions gives

         ∏                                   --∂--                 --∂--
SIJS = (  ′  ⟨⟨jJS′|iIS′⟩⟩)(bS⟨⟨jJS |iIS⟩⟩ + aS⟨⟨jJS|∂QS |iIS⟩⟩+ 2bS⟨⟨jJS|QS ∂QS  |iIS⟩⟩).

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 ⟨i|ŪS|i⟩, see (2.67). The electronic matrix element, on the other hand, arises from the off-diagonal elements ⟨j|ŪS|i⟩. As the diagonals and the off-diagonals of the coupling operator need not be related to each other, there is no direct link between the contribution of a mode to the electronic matrix element and its contribution to the vibrational part. In the literature, modes contributing to the electronic matrix element are called promoting modes and those contributing to the vibrational part are called accepting modes  [131]. Usually, every mode of a defect system will fall into both categories to some degree.

The target of the present work is to extract relevant parameters for the non-radiative transitions from atomistic calculations. Although the parameters aS and bS 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 ⟨j|ŪS|i⟩ required for the promoting contribution due to the ill-defined or inaccurate wave functions in the available electronic structure methods. Under this viewpoint, an approach to calculate the non-radiative transition rate using the accurate operator (2.98) appears unreasonable due to the inaccuracy of the underlying model system. We therefore assume that all modes in our system are purely accepting modes, which enter the transition rate only through the Franck-Condon factors ⟨ ⟨jJS|iIS⟩ ⟩ and the conservation of energy. We assume a perturbation operator ˆLthat contains the promoting contributions of the modal spectrum but acts like an optical matrix element Mˆ. Just as ˆL, Lˆis assumed to couple the initial and final state elastically, i.e. it does not contribute to the conservation of energy. As Lˆis not known exactly, it will be approximated by a parametrizable expression. The resulting non-radiative transition rate can thus be based on the radiative transition rates (2.93) as the limit [14111712386]

                           |      |2
Wi→f =  lim  Wi →f(hν) = 2π-||⟨f|ˆL′|i⟩||f (0).
       hν→0             ℏ

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

2.10  Vibronic Transitions with Classical Nuclei


Figure 2.8: A schematic configuration-coordinate diagram illustrating vibronic transitions in a classical and a quantum mechanical picture. In the classical picture (left), radiative as well as non-radiative electronic transitions occur with the nuclei fixed in their respective position (Franck-Condon principle). Thus, optical transitions occur vertically in this diagram and non-radiative transitions happen at the crossing-points of the potential energy surfaces. As the system most dominantly resides at the minima of the potential energy surfaces, different wavelengths are observed for absorption and emission of radiation. In the quantum-mechanical picture (right), the transitions happen between the vibrational sub-states arising from the different potential energy surfaces.

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 [129117]. 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(⃗R)⟩ and |f(⃗R)⟩ proceeds via the optical matrix element Mˆ. For any position R⃗ of the nuclei, one can now again employ time-dependent perturbation theory to calculate the transition rate as

              2π ||             ||2  (                  )
Wi →f(hν, ⃗R ) =--|⟨f(R⃗)|M ˆ |i(R⃗)⟩| δ hν - Ef(R⃗)- Ei(⃗R ) .

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  ⃗
R is hence given by

R) = Z-1e-  ⃗
EikB(RT) (2.104)
Z = R⃗e-   ′
Eik(⃗RBT) d⃗
R (2.105)
and thus the total transition rate for a given optical energy can be computed from
Wif() = NP(⃗R)Wif(hν,R⃗)d⃗R= NZ-1e-Ei(⃗R′)-
kBT Wif(hν,⃗R)d⃗R (2.106)
= NZ-1e-   ′
Eki(B⃗RT)- 2π-
 ℏ||  ⃗ ′ ˆ   ⃗′  ||
|⟨f(R )|M |i(R )⟩|2δ(        ⃗ ′      ⃗′)
 hν - Ef(R ) - Ei(R )d⃗
R (2.107)
Again assuming the electronic matrix element to be independent of the nuclear configuration one obtains
Wif() = 2π-
ℏ|      |
|  ˆ   |
|⟨f|M |i⟩|2f() (2.108)
f() = NZ-1e-Ei(⃗R′)-
kBT δ(          ′       ′)
 hν - Ef(⃗R )- Ei(R⃗)d⃗R. (2.109)
This time the line shape function is only determined from the probability to find the system in a configuration that fulfills the conservation of energy for the optical energy . For the non-radiative transitions we again replace the optical matrix element by the quasi-optical non-adiabacity operator ˆLand let the optical energy go to zero as illustrated in Fig. 2.8. In summary, when calculating vibronic transitions with classical nuclei, we again employ (2.102), but use (2.109) as the line shape function. This classical line shape is determined from a transition-state theory like expression, where the transition state is defined by the conservation of energy. In comparison to (2.94), the classical line shape has the advantage that it can in principle use the potential energy surfaces from electronic structure methods in their full glory, i.e. without the necessity of a harmonic approximation.

2.11  Solution of the Master Equation

SVG-Viewer needed.

Figure 2.9: Sketch of the stochastic simulation algorithm (SSA)  [114]. The algorithm generates a realization of the stochastic process described by the chemical master equation (2.20).

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 ⃗x, given that ⃗X = ⃗x0 at t = t0. For a system with a small set of states {⃗x1,,⃗xΩ}, a direct solution can be attempted, which results in a coupled system of Ω differential equations [37]. However, in many situations the number of states will be quite big or even infinite, rendering a direct solution of the master equation unfeasible or even impossible. A feasible alternative is the stochastic simulation algorithm (SSA)  [114] explained in Fig. 2.9, which is also known as the kinetic Monte Carlo method. Instead of solving the differential equation (2.20), a realization of the stochastic process X⃗ is generated using pseudo-random numbers. The SSA does not have any algorithmic parameters and is a mathematically exact description of the system defined by the states and reaction channels  [114]. The moments of the probability distribution P(⃗x,t) are trivially calculated by averaging over several simulation runs, although care has to be taken to ensure the randomness of the pseudo-random numbers between two runs to avoid correlation effects.

2.11.1  Subsystems, Well-Stirredness, and Diffusion

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.

    (  ⃗x1  )
    |  ⃗x   |
    ||    2 ||
⃗x = ||  ⃗x3.  || , ⃗x δ ∈ {⃗s1,⃗s2,...,⃗sΩ }, Θ ≫ Ω.
    (   ..  )

For every subsystem δ there is now a number Γ of reaction channels with associated state change vectors ⃗νδ1⃗νδΓ and propensity functions aδ1(⃗sδ)aδΓ(⃗sδ).

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 which dopant has contributed a carrier but only how many of them have. As the sub-systems are similar enough that their propensity functions are the same

a   =  a   = a ,
 δiγ    δjγ    γ

the chemical state of the total system can be reformulated, considering only the numbers of subsystems nω in a given state ⃗sω. The new total state vector will be

     (     )
     |  n2 |
 ′   ||  n  ||
⃗x  = ||  .3 || ,
     (  ..  )

the state change vectors will only operate on the numbers n1,,nΩ, and the corresponding propensity functions will take the form

a ωγ = n ωaγ,

corresponding to the probability per unit time that any of the subsystems in state ω will take reaction channel γ. In theoretical chemistry these reactions in which an isolated sub-system undergoes a transition are termed unimolecular reactions [115]. Especially if the chemical system under consideration consists of diffusing species, such as atoms or molecules in a solvent, two of these subsystems may also interact and form new species, for example 2H2 + O2 2H2O. In those reactions, the corresponding reaction channel changes two or more numbers of ⃗xat the same time and the associated propensity function takes the form

a′ω ω γ = nω1nω2aγ,
  1 2

which is the probability per unit time that any two of the nω1 specimens of species ω1 and nω2 specimens of species ω2 will undergo the reaction γ, where different species in the above context may refer to different types of sub-systems (e.g. different molecules) or different sub-systems of the same type in a different state. Reactions of this type are called bimolecular reactions. For bimolecular reactions which involve two specimens of the same species, the propensity function is

 ′     nω(nω---1)-
aω γ =     2     aγ,

which corresponds to the probability that any pair of specimens of type ω will undergo the reaction γ. Any higher order reactions can be reduced to a series of bimolecular and unimolecular reactions [115].

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 well-stirred. Different approaches exist for chemical systems where the intermixing is not fast enough, which are called reactive-diffusive systems. These approaches and their relevance for the reaction-diffusion model for NBTI will be discussed in Chap. 3.

2.12  Reaction Rate Equations

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 average numbers of (2.112) over time. These averages are obtained as

⟨⃗n⟩ (t) =    ⃗nP (⃗n,t|⃗n0,t0),

which corresponds to the first moment of the probability distribution P.

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 ⃗N = ⟨⃗n ⟩. The time evolution of ⃗N is calculated from

∂ ⃗N-  ∂- ∑                 ∑    -∂
∂t  = ∂t    ⃗nP (⃗n,t|⃗n0,t0) =    ⃗n∂tP (⃗n,t|⃗n0, t0).
          ⃗n                 ⃗n

The substitution of the right side of (2.20) for the derivative of P and a few trivial transformations yield [115]

∂ ⃗N-  ∑
∂t  =    ⃗νγ ⟨aγ⟩.

For the case of unimolecular reactions (2.113), whose propensity functions are linear in ⃗n, it is trivial to show that ⟨aγ⟩ = ⃗aγ N⃗, which leads to the popular description using reaction rate equations

∂-⃗N-=  ∑  ⃗ν (⃗a ⋅N⃗).
 ∂t        γ  γ

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 ⃗
N. Instead, the number of particles is usually normalized to the volume of the system, leading to a formulation based on concentrations instead. In this formulation, the particles within the system are understood as a continuum of a certain density. A formulation like this is usually termed a macroscale formulation.