4 Atomistic Modeling and the BTI Defect
The defect in the multistate multiphonon model for BTI has (at least) two structural and (at least)
two electronic configurations, as explained in Sec. 1.5. The structural reorganization of the
defect goes largely unnoticed by the rest of the system and it is only the charge state of the
defect that acts on the rest of the device via coulombic interaction. The barrier hopping
transitions 1 ⇌ 1′ and 2 ⇌ 2′ appear to proceed independently of the current state of the
device. Only the NMP transitions 1 ⇌ 2′ and 2 ⇌ 1′ need to consider the availability of
electrons and holes. All the states and transitions of the BTI defect can be understood
within the framework laid out in Chap. 2. The present chapter shows how these transitions
can be described using an atomistic model of the defect and a macroscopic model of the
device.
4.1 Atomistic Defect Models
An overview of atomistic defect models employed for the calculation of E′ defects in the literature is
given in Tab. 4.1, including the electronic structure method and the atomic representation employed.
The design of the atomic arrangement used to represent the defect directly influences the predictive
power of the calculations. Three approaches are widely used in the literature [102], which are
illustrated in Fig. 4.1:
 In isolated cluster calculations, the defect and its surrounding atoms are represented as
an isolated molecule. This method is often used in combination with highly sophisticated
quantum chemistry programs which offer highly accurate electronic structure methods.
 Embedded cluster methods embed the pseudomolecules into an atomic host structure
which is treated using empirical potentials.
 In supercell calculations, the defect is introduced into a periodic cell.
All three approaches have their own benefits and drawbacks. Isolated cluster models are relatively easy
to set up and one can choose from the rich pool of quantum chemistry programs for the electronic
structure calculation. Isolated clusters offer a simple way to treat defects that are hard to integrate
into a full atomic host structure. However, they suffer from the inherent neglect of electrostatic and
mechanical long range interactions.
Embedded cluster methods are an approach to improve on these drawbacks. In an embedded cluster
calculation the atomistic defect model is split into three nested regions. The inner region
surrounds the defect and is treated using a quantumchemistry method. This quantum region is
surrounded by the socalled classical region, which contains a large number of atoms whose
interactions are described using empirical potentials [66]. Finally, the quantummechanical
and classical regions are embedded in an infinite continuum model of the host material.
Electrostatic interactions between the quantum region and the classical region are considered
using the shell model for polarization on the classical atoms. In structural optimizations
an inner sphere of atoms, containing the quantum region and a fraction of the classical
region, is moved while the outer shell is held fixed. The atoms at the interface between
the quantum region and the classical region need to be included both in the electronic
structure part and in the empirical potential part of the calculation. The representation
of the interface atoms in the electronic structure method is especially problematic as at
least one valence electron of these atoms needs to be replaced by an interaction with the
empirical potential. To account for the missing valence electron, the interface atoms are
equipped with a parametrizable pseudopotential that acts on the electrons of the quantum
region. Embedded cluster methods improve upon the shortcomings of the isolated cluster
calculations at the expense of a dramatically increased complexity. The setup of an atomistic
defect model using the embedded cluster approach requires a lot of finetuning of both the
empirical potential and the pseudopotential for the interface atom. For this reason, the
embedded cluster approach is not suitable for the original study of defect parameters, but is
more appropriate for refining investigations on defects for which reference calculations
exist.
Supercell defect calculations have become quite popular in the solid state community. This
popularity comes in part from the existence of conduction and valence levels in the calculation, and
the increasing availability of plane wave based general purpose codes [102]. In the present work, we
employ a supercell approach for exactly these reasons. Issues of the supercell approach include the
treatment of charged cells, which will be discussed later, and interactions between the periodic
images.
Reference  Electronic Structure  Boundary Condition  Host Structure 
[54]  LCLOMO  isolated cluster  αquartz 
[55]  TightBinding  isolated cluster  αquartz 
[56, 57]  MINDO/3, MOPN  isolated cluster  αquartz 
[58]  LSDA  supercell  αquartz, amorphous 
[159]  HF, MP2  isolated cluster  αquartz 
[59]  PBE  supercell  αquartz 
[160, 161]  B3LYP  isolated cluster  αquartz 
[162]  HF, MP2  isolated cluster  Free relaxation 
[77]  HF  embedded cluster  αquartz 
[60]  B3LYP  isolated cluster  amorphous 
[63, 163]  HF  embedded cluster  αquartz 
[61]  PBE  supercell  amorphous 
[64]  LSDA  supercell  αquartz, amorphous 
[65]  PBE  supercell  amorphous 
[62]  DFT  supercell  αquartz 
[78, 68]  HF  embedded cluster  amorphous 
[164, 69]  B3LYP  embedded cluster  amorphous 
[165]  PBE+HF  supercell  αquartz 

Table 4.1:  An overview of electronic structure methods and atomistic representations
employed in published defect studies. The electronic structure methods are MINDO
(modified intermediate neglect of differential overlap), HF (HartreeFock), MP2 (second order
MøllerPlesset perturbation theory), PBE (DFT with gradientcorrected functional due to
Perdew and coworkers), B3LYP (hybrid DFT with gradientcorrected functional due to Becke
and coworkers with HF exchange), and PBE+HF (hybrid DFT with PBE functional and HF
exchange). 
The host materials used in the studies of E′centers are mostly αquartz and amorphous SiO2. Due
to its similar density, αquartz is usually considered a reference system for amorphous silica
[59, 57, 163]. Atomistic amorphous silica models are continuous random network structures
generated by molecular dynamics melting and quenching [166, 167]. Although an amorphous host
structure is closer to the physical reality of the oxide in an MOS transistor, the result of the molecular
dynamics generation depends strongly on the parameters of the processing, especially the
cooling rate [110]. As the defect parameters strongly depend on the surrounding atoms, the
accuracy of a defect study depends on a host structure that accurately represents the material
under study. The generation of an amorphous host structure thus requires careful testing of
the dependence of the structural properties on the creation process and their evaluation
against experimental data. The present work, however, calculates quantities that have not
been investigated before. Therefore, it is reasonable to create our defect models upon an
αquartz host structure, in order to build a reference for later studies of defects in amorphous
host structures. Due to the vast amount of literature available for this system, αquartz
seems to be best suited for these reference calculations. Our work follows Blöchl [84, 59]
and uses an orthogonal 72 atom αquartz supercell [168]. The employed functional is the
gradientcorrected functional of Perdew, Burke, and Ernzerhof [101], the core electrons are
represented using the projector augmented wave method as implemented in the Vienna
abinitio simulation program (VASP) [169, 170]. The electronic structure is expanded
using a planewave basis set, employing periodic boundary conditions for the simulation
cell.
4.1.1 Parameters of the DFT Calculations
To be compatible with the work of Blöchl [59], we have constructed an orthogonalized
αquartz supercell for our defect calculations. We use the coordinates given in [168], which
reflect the hexagonal symmetry of the αquartz as a nonorthogonal unit cell, see Fig. 4.2.
As two angles of the cell are already 90^{∘}, the orthogonalization of the cell concerns only the plane of
the γ angle. The procedure is described in Fig. 4.3.
After the construction the atoms in the cell were relaxed keeping the cell shape and volume
constant. The resulting structure is the basis for all following defect calculations and is shown in
Fig. 4.4.
The energy cutoff for the plane wave expansion is 800eV for all calculations in this work. VASP’s
real space projection is disabled to improve the accuracy. The k space is only sampled at the Γ point.
Structures are optimized using the quasiNewton algorithm [171, 96] until the forces are below
10^{2}eV∕Å. Optimizations starting far from the optimum were preoptimized using the more robust
conjugate gradient algorithm.
4.1.2 Defect Structures
As explained in Sec. 1.6, the defect structures investigated in this work are the oxygen vacancy and
the hydrogen bridge. The oxygen vacancy defect was created by removing one oxygen atom
from the cell, the hydrogen bridge by replacing an oxygen atom with a hydrogen atom.
Subsequent optimizations led to the dimer configuration for the oxygen vacancy and the closed
configuration for the hydrogen bridge, see Fig. 4.6 and Fig. 4.7. The puckered oxygen
vacancy and the broken hydrogen bridge were constructed by manual rearrangement of the
defects.
Both defects give rise to one occupied state in the SiO2 band gap. This state is the highest occupied
state of the KohnSham spectrum. As the eigenstates of the auxiliary system are occupied from the
bottom up, the positively charged defects are easily created by removing an electron from the
calculation. Charged defects, however, are problematic in supercell calculations, as it is not possible to
define a meaningful energy for a periodic cell with a net charge. Thus, when an electron is removed
from the DFT calculation a compensating background charge is automatically added by VASP
[102, 171]. The interaction of the charged defect with this artificial charge background
has to be accounted for. In the present work we employ the correction used by Blöchl
[59], which simply increases the energy for charged cells by 0.48eV, see the reference for
details.
In the neutral state, the dimer configuration of the oxygen vacancy is more stable than the
puckered configuration by 2.966eV. In the positive state, the situation is reversed with the
puckered configuration being slightly lower in energy than the dimer configuration by 88meV.
For the neutral hydrogen bridge, the DFT predicts the closed configuration to be 21meV
more stable than the broken configuration. Similar to the oxygen vacancy, the positive
hydrogen bridge is 185meV more stable in the broken configuration than in the closed
configuration.
The optimized structures for the oxygen vacancy and the hydrogen bridge are shown in Fig. 4.6 and
Fig. 4.7, respectively. The geometries obtained from our calculations are in excellent agreement with
the literature as shown in Tab. 4.2 and Tab. 4.3.
 dimer  puckered

 SiSi distance  SiSi distance  SiO(3) distance

 neutral  positive  neutral  positive  neutral  positive 
      
present work  2.437 Å  2.964 Å  4.026 Å  4.358 Å  1.887 Å  1.841 Å 
Blöchl [59]  2.445 Å  3.011 Å  4.052 Å  4.358 Å  1.907 Å  1.852 Å 
Boero [58]  2.520 Å  3.050 Å  ≈ 4.37^{†} Å  4.46 Å   1.82 Å 
Pacchioni [159]  2.530 Å  3.024 Å  

^{†}…sampled from fig. 1
Table 4.2:  Comparison of selected structural parameters of the oxygen vacancy defect in our
calculations with values from the literature. 
  neutral  positive

  present work  Blöchl [59]  present work  Blöchl [59] 
     
closed  SiSi distance  3.351 Å  3.368 Å  3.206 Å  3.225 Å 
 Si _{
1}H distance  1.522 Å  1.534 Å  1.641 Å  1.646 Å 
 Si_{2}H distance  1.985 Å  1.953 Å  1.689 Å  1.710 Å 






broken  SiSi distance  3.948 Å  4.011 Å  4.389 Å  4.394 Å 
 Si_{1}H distance  1.460 Å  1.460 Å  1.471 Å  1.477 Å 
 Si_{2}H distance  2.698 Å  2.742 Å  3.171 Å  3.088 Å 
Table 4.3:  Comparison of selected structural parameters of the hydrogen bridge defect in our
calculations with values from the literature. 
It is common in defect studies to calculate the “formation energy” of a defect. This energy is an
indicator for the abundance of the defect in thermal equilibrium. The structure of the oxide of an
MOS transistor arises from the growth kinetics of the oxidation process. Although the vibrational
state of the atoms is in thermal equilibrium with the environment, the bonding structure is not a
thermal equilibrium arrangement. Thus, the significance of the formation energies calculated here
for the abundance of oxygen vacancies and hydrogen bridges in MOS oxides is limited.
We calculate the formation energies here for the purpose of comparison with published
work. For the calculation of the formation energies it is necessary to define a “reservoir
energy” for each atomic species. The reservoir energy of the atoms is usually defined as
the gas phase energy, i.e. an isolated atom calculation [159, 66]. The formation energy
of the oxygen vacancy E_{fOV} and the hydrogen bridge E_{fHB} have been calculated with
respect to the energies of the isolated oxygen atom E_{O} and the isolated hydrogen atom E_{H}
as
E_{fOV}  = E_{OV}  E_{αQ} + E_{O}  (4.1)

E_{fHB}  = E_{HB}  E_{αQ} + E_{O}  E_{H}  (4.2) 
where E_{OV}  E_{αQ} and E_{HB}  E_{αQ} are the optimized total energies differences between the DFT
calculations with the defect and the ideal crystalline host. The formation energy for the oxygen
vacancy in our calculations is E_{fOV} = 8.468eV, in agreement with earlier works [59, 159].
Considerably lower formation energies of about 4.6eV have been obtained using the embedded cluster
method [66]. This difference partially arises from the neglect of correlation in the HartreeFock
method used for the quantum region in [66], but also comes from the larger number of relaxing atoms
in the embedded cluster approach. This indicates that some strain may remain in the employed
supercell and the results of this work should be compared against calculations involving larger
supercells at some point. However, this check is not part of the present work. The formation energy
of the hydrogen bridge defect in our calculations is E_{fHB} = 7.802eV, in agreement with
[59].
4.2 Barrier Hopping Transitions
In terms of Chap. 2, the structural reorganization of the defects, i.e. the transitions 1 ⇌ 1′ and
2 ⇌ 2′ are barrier hopping transitions. A detailed study of these transitions requires the solution
of the thermal integrals of the transition state theory, which gives the barrier hopping
transition rates. However, as discussed in Sec. 2.8, one can already gain some insight by
just calculating the activation energy for the process, which we consider sufficient for the
evaluation of a candidate BTI defect. The calculation of the activation energy requires to
find the path of minimum energy between the two potential energy surface minima that
correspond to the start and endpoint of the transition. A popular method to find the reaction
path is the nudged elastic band (NEB) method [171, 172], which is also applied in the
present work. In the NEB approach, configurations along the reaction path are connected via
virtual spring potentials. Then the system of all configurations is optimized with respect to
the total energy consisting of the energies of the individual configurations plus the spring
potential.
Our NEB calculations use ten points to sample the reaction path. The path is preoptimized using
the conjugate gradient algorithm until forces are below 1eV∕Å. The final optimization to reduce the
forces below 0.1eV∕Å is done using the quasiNewton algorithm. The resulting energies on the
optimized reaction paths for both defects are shown in Fig. 4.8.
4.3 Macroscopic Device Simulation
The charge state transitions 1 ⇌ 2′ and 2 ⇌ 1′ arise from the interaction of the defect with the
semiconductor device. Their theory consequently needs to encompass both the atomistic model of the
defect and the macroscopic device model. The atomistic level is the main focus of this
thesis and has been described extensively in this document. An indepth discussion of the
fine art of macroscopic device simulation goes beyond the scope of this document and so
this topic is only discussed to the extent that is absolutely necessary for the following
calculations.
The methods of semiconductor device simulation are used to model the occupation dynamics of the
extended electronic states k inside a semiconductor [173]. In the ground state some of these states
are occupied, and some are unoccupied. The occupied extended states kv form the valence band of
the semiconductor, the unoccupied states kc form the conduction band. Their occupation is
described by the occupation probability fe(k) as
 (4.3) 
The time evolution of the electronic system of the semiconductor is described by the time evolution of
the occupancy function fe(k). On the energy scale, the valence band and the conduction band of a
semiconductor are separated by an energy range that has no stationary states, which is called the
band gap. In its ground state, the electronic system of the semiconductor is unable to carry electric
currents. Only when an electron is removed from the valence band (fe(kv) < 1), or added to the
conduction band (fe(kc) > 0), a transport of charge becomes possible. These electronic
configurations are inherently many body excitations, which include not only the nongroundstate
occupancy but also the response of the electronic system to this occupancy. In the theory of
semiconductors, these many body excitations are described as quasi particles [174], where additional
electrons in the conduction band are called “excess electrons”, or electrons for short, and the
unoccupied states are called “holes”, or defect electrons. For the excess electrons, the quasi
particle states {w_{n}} are defined, which correspond to occupied conduction band states
kc. For the holes, the quasi particle states {w_{p}} correspond to unoccupied conduction
band states kv. The occupancies for the excess electron and hole states are then fn and
fp, respectively. The physics of the states {w_{n}} and {w_{p}} is determined by the band
structure. This band structure defines the total energy of each quasi particle state as a
function of its crystal momentum . At material interfaces such as the semiconductoroxide
interface in MOS devices, the band structure abruptly changes, leading to energetic barriers
for the carrier gas. It is common to draw the local band edges Ec0() (conduction band
edge) and Ev0() (valence band edge) as a function of position. These graphs are called
band diagrams and are a popular tool to illustrate processes in semiconductor devices, see
Fig. 4.9.
In semiconductor device modeling, the materials that constitute the semiconductor device are
described using their macroscopic properties, e.g. the permittivity [175]. The coulombic interactions
in the device are accounted for in a mean field fashion by the potential Φ(), which is obtained from
the Poisson equation
 (4.4) 
that accounts for the charge arising from the mean density of holes
 (4.5) 
the mean density of electrons
 (4.6) 
and the ionized dopants ND() and NA() [175]. This potential contributes to the energy of the
particles, which can be accounted for in the band diagram by replacing the band energies
by
 (4.7) 
The dependence of the local band structure on the electric potential is usually termed “band
bending”. The band bending as obtained from a device simulation is shown in Fig. 4.10
The quasi particle states w_{n} and w_{p} as well as their occupancies are determined by the carrier
model. Depending on the complexity of the device structure, the dynamics of the excess
electrons and holes can be described at different levels of physical accuracy. The most
common carrier models, which are also employed in stateoftheart TCAD simulators, are
based on semiclassical electrons and holes. These semiclassical particles are pictured as
wave packets of Bloch wave functions moving on classical trajectories determined by the
local band structure. The time evolution of the carriers is described using the Boltzmann
transport equation [176] in various approximations [177, 178]. The state of the carrier gas in
this description is fully characterized by the distribution f_{n}(,) of the electrons and the
distribution f_{p}(,) of the holes in the phase space (,). The current transport in most
semiconductor devices can be described very accurately using classical carrier models. For our
purposes, however, their applicability is somewhat limited by the inherent neglect of quantum
mechanical effects such as the quantization in the inversion channel and the penetration of
carriers into the oxide through tunneling. While this quantization has only little effect on the
prediction of transport properties, it has a profound effect on the energetics of the carrier
system, which strongly influences the carrier trapping, as will be discussed in the next
section.
Quantum mechanical carrier models can be loosely categorized into closed and open
boundary models. In the former, the Schrödinger equation for the {w_{n}} and {w_{p}} is
solved with Dirichlet boundary conditions. The occupancies of the excess electron and hole
states are derived from thermal equilibrium. As before, the carrier model has to be solved
selfconsistently with the Poisson equation [179], hence these carrier models are sometimes called
SchrödingerPoisson models. Due to the assumed equilibrium occupancy of the particle states
these models give good results for nonconducting situations. In bias temperature stress
experiments, apart from the switches between stress and recovery voltage, the only currents
flowing are the gate current and the monitoring drain current. As these currents have
negligible impact on the states and their occupancies the carriers are typically in thermal
equilibrium with the contacts, so a closed boundary description is justified in principle.
SchrödingerPoisson models have been employed in calculations of defectdevice interactions to
model the quantization in the inversion region of the MOS structure and the penetration
of the quasi particle states into the oxide [180, 181, 182]. The boundary conditions of
the closed boundary quantum mechanical carrier models are of special importance for
the calculation of carrier trapping at defects in MOS oxides. As is more deeply discussed
below, carrier trapping depends on the ability of the carriers to penetrate into the oxide via
quantum mechanical tunneling. The boundaries in SchrödingerPoisson calculations have to
be placed either at the oxidegate or the oxidesubstrate interface and thus potentially
influence the tunneling behavior of the carriers and consequently the predicted trapping
rates.
Open boundary quantum mechanical carrier models are the most accurate description of the carriers
in the device. These models treat the device essentially as a scattering problem. The most popular
open boundary transport models are based on the framework of nonequilibrium Green’s functions
(NEGF) [176]. The biggest advantage of open boundary over closed boundary models is the absence
of artificial boundaries inside the device.
Independent of the physical details of the carrier model, the quantities in a semiconductor device
simulation are the average potential Φ(), the excess electron and hole states {w_{n}} and {w_{p}} as
well as their occupancies f_{n}(w_{n}) and f_{p}(w_{p}). As for many applications the most interesting property
of an electronic state is its energy, the excess electron and hole states are usually replaced by the
corresponding densities of states [176]
D_{n}(x,E)  = ∫
_{wn}^{2}δ(E  E_{
wn′})dw_{n}′  (4.8)

D_{p}(x,E)  = ∫
_{wp}^{2}δ(E  E_{
wp′})dw_{n}′  (4.9) 
and the corresponding occupation densities f_{n}(,E) and f_{p}(,E). These quantities are also of central
importance for the calculations of the following section.
4.4 NonRadiative Transitions in the DeviceDefect System
For the calculation of nonradiative transitions between the defect and the device we start from the
considerations in Sec. 2.9.2, where the general theory of NMP transitions has been laid out. The
transition rates are calculated using (2.102), which requires the determination of the electronic matrix
element and the line shape function, where the latter can either be described in the quantum
mechanical (2.94) or the classical form (2.109). To describe the phonon mediated capture and emission
of electrons and holes, it is necessary to specify the initial electronic state i and the final electronic
state f for the devicedefect system. The electronic structure of this system is essentially a
quantum mechanical many body problem, which we will formulate upon a basis of single
particle states. In the following, d is the localized orbit at the defect site, and k are the
electron states in the semiconductor device, which act as a reservoir for the transition. The
defect state as well as the reservoir states can be either unoccupied or occupied by an
electron.
The mathematical framework for this is the second quantization formalism [183]. In the following,
the creation operators ĉ†, the destruction operators ĉ†, and the number operators = ĉ†ĉ† are
used .
In accord with the derivations in Sec. 2.9.2, the electronic basis states for the vibronic transition are
BornOppenheimer states, which are mixed by the effective nonadiabacity operator ′. The adiabatic
Hamiltonian of the system Ĥ consists of the Hamiltonian of the defect Ĥd and of the device
ĤD
 (4.10) 
The energy of the device is determined by the occupancy of the free states k in the valence
band and the conduction band. The energy E_{k} of these states comes from the employed
carrier model and can be understood as a functional of the mean field of the carrier gas,
yielding
 (4.11) 
The energy contribution from the electronic structure of the defect is the potential energy surface
that we obtain from our atomistic model. In the neutral (occupied) state, the potential
energy surface of the defect is E0d() and in the positive (unoccupied) state it is E+
d (). In
addition, we have to consider the interaction of the defect with the mean potential Φ(). As
we approximate the defect as a point charge at the position d, this interaction is just
q0Φ(d) in the positive state and zero in the neutral state. The resulting defect Hamiltonian
reads
 (4.12) 
The effective nonadiabacity operator mixes the electronic states. In the devicedefect system, it
annihilates an electron in a reservoir state and creates an electron at the defect state. The
corresponding single electron operators are denoted by ′′. The nonadiabacity operator thus
reads
 (4.13) 
Now that all the operators are defined, one can calculate the rates for the capture and
emission of electrons from the reservoir. At first we consider the capture transition from
a reservoir state k′, which corresponds to the removal of an excess electron from the
conduction band or emission of a hole into the valence band of the device. Initially, d is
unoccupied, k′ is occupied, and the occupancy of the other reservoir states is given by f_{k}. In
the second quantization formalism, the wave functions are constructed by applying the
creation operators to the vacuum state ∅. The electronic system changes from the initial
state
 (4.14) 
to a state
 (4.15) 
The energies of the initial and final state read
 (4.16) 
and
 (4.17) 
respectively. Interestingly, the energy of the captured particle as well as the mean potential at the
defect enter as an energetic shift between the potential energy surfaces before and after the capture
transition. Thus, for every reservoir state, a different configuration of the potentials is obtained leading
to a different NMP transition rate, see Fig. 4.11. The dependence of the energy shift on the mean
potential is responsible for the large bias dependence of the BTI defect [32, 85] as mentioned
briefly in Sec. 1.5. The dependence of the energetic shift between the potential energy
surfaces on the potential is a general feature of NMP transitions that is also present at
defects within semiconductors [144, 184]. However, for defects in the oxide this effect
dominates the field dependence of the transition as the oxide due to the large tunneling
distance between the free state and the defect. The shift induced by the oxide field in
the limit of a nearly uncharged oxide can be estimated by a simple model [185]. In this
case, the potential grows or falls linearly with the distance d from the substrateoxide
interface and with the field F. Neglecting the dependence of the free states E_{k} on the
potential, the shift induced by the oxide field is just q0dF. This estimate is of course only
accurate to firstorder. When considering real devices, the potential at the defect site and the
dependence of the E_{k} on the potential have to be taken from a selfconsistent calculation.
The capture rate is now calculated from (2.102) as
 (4.18) 
with
 (4.19) 
As mentioned in Sec. 2.9.5, ′′ is not accurately known, and so this term has to be estimated.
Different approaches for the estimation of the matrix element exist in the literature, ranging from
simple tunneling expressions to model defect wave functions [144, 180, 26]. Due to the large extent of
the free states of the reservoir compared to the extent of the localized state, we approximate the
electronic matrix element as [186]
 (4.20) 
where λ is a parameter that needs to be calibrated to experimental data. In essence, this compresses
the electronic state of the defect into a Dirac peak and reduces the electronic matrix element to a
tunneling expression, which is assumed to be a reasonable approximation. The use of a more
natural wave function as is sometimes employed in the literature enables the modeling of
crystal momentum dependent coupling to the free states, which is not considered at the
moment.
The quantum mechanical line shape is defined as (2.94)
 (4.21) 
The energies EfF and EiI are obtained from the solutions of the vibrational Schrödinger equation
using the electronic energies E_{f}() and E_{i}() as potentials,
(N + ∑
_{k≠k′}f_{k}E_{k}[Φ] + E+
d () + E_{k′}[Φ] + q0Φ(d))iI  = EiIiI ,  (4.22)

(N + ∑
_{k≠k′}f_{k}E_{k}[Φ] + E0
d())fF  = EfF fF .  (4.23) 
In these potentials, the device level quantities only enter as a constant shift of the energies. This shift
has no influence on the eigenvectors iI and fF , and trivially adds to the eigenvalues. We can
thus define the eigensystems
 (4.24) 
and
 (4.25) 
which only depend on the quantities of the atomistic model. For electron capture EiI and EfF can
then be written as
EiI  = E_{+I} + ∑
_{k≠k′}f_{k}E_{k}[Φ] + E_{k′}[Φ] + q0Φ(d)  (4.26)

EfF  = E_{0F } + ∑
_{k≠k′}f_{k}E_{k}[Φ],  (4.27) 
and the vibrational wave functions are
 (4.28) 
Inserting these expressions into (4.21) gives
 (4.29) 
with
 (4.30) 
The line shape as a function of energy is determined by the FrankCondon factors 0F+I as well
as E_{0F }  E_{+I}. All those quantities are local to the defect and can be extracted from the selected
electronic structure method, as shown in Sec. 4.6. The energy conservation expression provides that
the energy of the captured particle is counterbalanced by an equivalent increase or decrease of
the vibrational energy. The energy of the captured or emitted carrier in nonradiative
transitions thus plays the same role as the optical energy in (2.94). For the use of the line
shape in a device simulation, it is favorable to reference the energy scale of (4.30) to the
carrier energy scale present in the macroscopic device model. We choose the local valence
band edge Ev0(d) of the host material of the defect, which is SiO2 in our case, for this
purpose
 (4.31) 
In this case (4.29) becomes
 (4.32) 
where Ev(d) accounts for the band bending as described above. In essence this means that the line
shape follows its local energetic environment and can be referenced to the local band edges at the
defect site, as illustrated in Fig. 4.12.
An analogous derivation can be done for the emission of an electron. In this case the electronic
states are
i  = (ĉ†_{
d} + ∑
_{k≠k′}f_{k}ĉ†_{
k})∅  (4.33)

f  = (ĉ†_{
k′} + ∑
_{k≠k′}f_{k}ĉ†_{
k})∅.  (4.34) 
The relations of the initial and final energies and vibrational wave functions to the defect eigensystem
become
EiI  = E_{0I} + ∑
_{k≠k′}f_{k}E_{k}[Φ] and  (4.35)

EfF  = E_{+F } + ∑
_{k≠k′}f_{k}E_{k}[Φ] + E_{k′}[Φ] + q0Φ(d),  (4.36) 
and
 (4.37) 
The resulting expressions for the electron trapping and detrapping rates with respect to a state k′
are
W_{k′→d}  = λ^{2}f(+∕0)
v (E_{k′}[Φ]  Ev(d)),  (4.38)

W_{d→k′}  = λ^{2}f(0∕+)
v (E_{k′}[Φ]  Ev(d)),  (4.39) 
with
 (4.40) 
The prerequisite for an electron capture from the state k′ is that this state is occupied. On
the other hand The prerequisite for electron emission into a state k′ is that this state is
empty. In a semiconductor device at finite temperature, the occupancy of the quasiparticle
states fluctuates randomly due to the motion of the particles and the various scattering
processes. Therefore, for a real device one has to consider the probability fe(k′) of k′ being
occupied for electron capture and the probability 1  fe(k′) of k′ being unoccupied for
emission. Following the explanations of the previous section, these probabilities, along with the
probability ^{2} to find the extended state k′ at the defect site, is obtained from the device
simulation. The rates for the capture and emission of holes and electrons are then given
as
k_{n→d}  = λ∫
_{∞}^{∞}n(
d,E)f(+∕0)
v (E  Ev(d))dE,  (4.41)

k_{p→d}  = λ∫
_{∞}^{∞}p(
d,E)f(0∕+)
v (E  Ev(d))dE,  (4.42)

k_{d→n}  = λ∫
_{∞}^{∞}(D_{
n}(d,E)  n(d,E))f(0∕+)
v (E  Ev(d))dE,  (4.43)

k_{d→p}  = λ∫
_{∞}^{∞}(D_{
p}(d,E)  p(d,E))f(+∕0)
v (E  Ev(d))dE,  (4.44) 
where the particle densities relate to the density of states as
 (4.45) 
Our approach to the calculation of NMP capture or emission rates using both a macroscopic model of
the device and an atomic electronic structure model of the defect thus contains the following
steps:
 Calculation of the line shape functions f(+∕0) and f(0∕+) from the atomistic model.
 Calculation of Ev(), n(,E), D_{n}(,E), p(,E), and D_{p}(,E) form the macroscopic device
model.
 Combining the two using the expressions above.
4.5 Energy Levels
Before we can calculate the line shapes (4.67) and (4.68), it is necessary to establish a relation between
the potential energy surfaces of the atomistic model and the energy scale of the device model. This is
related to the task of energy level calculation from an electronic structure method. As it turns out,
there are several ways to define an energy level for a given atomistic defect model, which can be a
magnificent source of confusion for the communication between theorists working on atomistic
modeling and those working on device modeling. Therefore, we make an attempt here to clarify the
relations of the defect levels usually given in publications to each other and to our NMPtransition
based approach.
The most commonly published type of energy level calculated from atomistic defect models is the
thermodynamic trap level. It is also known as thermodynamic transition level in solid state theory
and is equivalent to the electron affinity (for reduction) and the ionization potential (for oxidation) in
theoretical chemistry [99]. The thermodynamic transition level is an equilibrium property of a defect
and comes from viewing the point defect as a thermodynamic system coupled to a reservoir
of electrons with a chemical potential μe. Neglecting the effects of pressure and entropy
and identifying the thermodynamic internal energy with the BornOppenheimer energy
Ed^{q}(_{q}^{0}) at the optimum configuration of the defect _{q}^{0} for the charge state q of the
system ,
one can define the formation energy [102]
 (4.46) 
This formation energy accounts for the average work necessary to exchange an electron with the
reservoir. At a certain μe, the charge state with the lowest Ed^{q} is the dominant charge state in thermal
equilibrium.
As can be seen from (4.46), the charge state formation energy depends linearly on the electronic
chemical potential and different charge states have a different slope q and a different shift. The
situation for a hypothetical defect is illustrated in Fig. 4.13. For different intervals of the chemical
potential a different charge state has the lowest formation energy, i.e. it dominates in thermal
equilibrium. The crossing points of the formation energy lines are the values of the chemical potential
where the dominant charge state changes. These are the thermodynamic transition levels.
The thermodynamic transition level from charge state q_{1} to charge state q_{2} is calculated
as
 (4.47) 
In thermal equilibrium, the thermodynamic transition level determines the dominant charge state of a
defect, see Fig. 4.14.
In semiconductor defect modeling, the chemical potential is commonly replaced by the electronic
Fermi level. For defects with vibrational internal states, as considered in the present work, different
potential energy surfaces for the various charge states lead to different vibrational spectra. In this
case, a change from one charge state to another leads to a difference in the vibrational entropy of the
defect [187] and thus the thermodynamic transition level of the defect becomes a Gibbs free
energy [188]
 (4.48) 
The change in the vibrational entropy of the defect thus induces a temperature dependence of the
transition level. As usually the entropy change between two charge states is small, the transition level
calculated from (4.47) is commonly considered a reasonable approximation. In the present work, we
follow the standard approach to neglect entropy for the calculation of transition levels. A special
property of the thermodynamic transition levels amongst the defect levels discussed here is that they
can describe transitions that involve the exchange of more than one electron with the reservoir. This
behavior is found in negativeU defects as illustrated in Fig. 4.15. In SiO2, a negativeU behavior has
been found for the interstitial hydrogen atom [59], which is also of significance for the present work, as
discussed below. Chadi [62] has predicted a negativeU behavior with a thermodynamically unstable
positive state for the oxygen vacancy. In BTI experiments the thermal equilibrium is only
relevant for the initial condition and this equilibrium is broken by the stress. Thus, the
thermal instability of the positive state of the oxygen vacancy is not part of the following
considerations, as these are concerned with the nonequilibrium behavior of the defect.
In addition to the thermodynamic transition level, which describes the thermal equilibrium of the
defect, a so called switching trap level or optical trap level can be defined. This trap level
describes the elastic exchange of an electron between the orbit of the defect and a different state. In
accord with the classical FranckCondon principle these transitions are assumed to take place with the
nuclear coordinates fixed [185, 59, 189, 190, 191]. Before the transition, the nuclear system
is assumed to be in thermal equilibrium, and thus the most dominant configuration is
the optimal configuration for the initial charge state. Contrary to the thermodynamic
transition level, the switching trap level can only be meaningfully defined for transitions
involving the exchange of exactly one electron as it refers to an elementary reaction. Again,
the energy of the defect plus the target state have to be considered. However, contrary
to the derivation of the thermodynamic transition level the energies here are the actual
electronic BornOppenheimer energies instead of thermodynamic potentials. For an elastic
transition the total energy before and after the electron exchange are the same. Energy balance
gives
 (4.49) 
for q∕q + 1 transitions and
 (4.50) 
for q∕q  1 transitions involving a state with energy E_{k′}. Defining the switching trap level E^{q1∕q2} from
the condition E_{sw}^{q1∕q2} = E_{k′}, one obtains
 (4.51) 
for q∕q + 1 transitions and
 (4.52) 
for q∕q  1 transitions. The switching trap level is a special case of a vibronic transition for negligible
phonon energy, i.e. T → 0, as can be easily seen from the considerations in Sec. 2.10. As the
temperature decreases, the probability density P(_{q}^{0}) in (2.104) for the nuclei to reside at the
minimum of the potential energy surface Ed^{q} grows without bounds while integral of the probability
density is normalized due to the partition function (2.105). Thus, in the zero temperature limit we
find
 (4.53) 
In this case, the classical line shape function (2.109) becomes
 (4.54) 
and via the capture and emission line shapes (4.67) and (4.68), the energy conservation expressions
(4.52) and (4.51) are obtained.
The relation of the discussed transition levels to the potential energy surfaces calculated from an
atomistic defect model are illustrated in Fig. 4.17.
Just as the line shapes, the thermodynamic and the switching levels require the definition of a
reference energy in order to set them into relation with the properties of the host material.
The density functional calculations in different charge states implicitly contain their own
reference level for the electrons. In the calculation of isolated molecules, the removal of
an electron from the calculation corresponds to taking the electron infinitely far away
from the molecule [99]. Thus, in isolated molecule calculations the reference level is the
vacuum level. While this referencing is reasonable for isolated molecules, in an infinite
periodic lattice the vacuum level is not a welldefined quantity and it is usually favorable to
select as reference energy one of the band edges or the midgap energy. All these energies,
however, are absent in the atomic defect models. There are different approaches to obtain
reasonable relations between the total energies of the electronic structure method and the band
edges. In principle one can use the valence or conduction level from the density functional
calculation [189, 192] as the reference. Unfortunately these levels come from the auxiliary
system and suffer from the band gap problem. A very popular method for the alignment of
the defect levels with the band edges is the socalled marker method [102, 59]. For this
approach it is necessary to have a defect with a thermodynamic transition state E_{th} that is
known from experiment relative to a reference level, for instance the valence band edge
Ev
 (4.55) 
This defect is called the marker. The same thermodynamic transition level of the marker defect is
calculated from the density functional calculation using (4.47). This level will be relative to the DFT
zero level
 (4.56) 
Assuming that the transition level was measured and calculated accurately, the experimental
valence level in the DFT energy scale is then obtained by subtracting the two expressions
as
 (4.57) 
In our work we follow Blöchl, who uses the thermodynamic (+∕) transition level of the hydrogen
interstitial in SiO2 as the marker. This level has been experimentally determined to be 0.2eV above
the silicon midgap [59]. Using the Si band gap of 1.14eV and the SiSiO2 valence band offset of
4.4eV, the SiO2 valence level is placed 0.2eV + 1.14∕2eV + 4.4eV = 5.17eV below the hydrogen
interstitial (+∕) transition level.
 Oxygen vacancy  Hydrogen bridge

 Dimer  Puckered  Closed  Broken 
(+/0) thermodynamic level  1.761eV  4.814eV  5.161eV  5.366eV 
(0/+) switching trap level  1.109eV  3.399eV  3.980eV  3.354eV 
Table 4.4:  Thermodynamic and switching defect levels for the oxygen vacancy and the hydrogen
bridge in both structural configurations, referenced to the SiO2 valence level. 
Tab. 4.4 shows the thermodynamic and switching defect levels in our defect calculations using
Blöchl’s energy scale alignment scheme. The calculated values are in agreement with results from the
literature [59, 65, 190, 193].
As a final remark, it is worthwhile to discuss one further defect level that is sometimes given in the
literature. This defect level is the eigenenergy, which results from the solution of the Schrödinger
equation in the electronic structure model of the defect. For the HartreeFock mean field
theory, Koopman’s theorem [96] states that the orbital energies are an approximation to
the change in the total electronic energy upon addition or removal of an electron to an
unoccupied state or from an occupied state, respectively. The values are approximate in
the sense that they do not consider orbital relaxation. Although Koopman’s theorem is
specific to the HartreeFock theory, similar theorems exist for density functional theory, so
the orbital energies can be taken as an approximation to the switching levels described
above. Fig. 4.18 shows the KohnSham spectra for the defect structures studied in the
present work. Indeed these spectra show the same trends as the calculated switching (0/+)
defect levels. The KohnSham eigenlevel of the oxygen vacancy in its dimer state is very
close to the valence band edge, all other occupied defect levels are deep in the band gap.
4.6 Extraction of the Quantum Mechanical Line Shape Functions
Once an atomistic defect model and an energy alignment method is selected it is possible to calculate
line shape functions f(+∕0)
v and f(0∕+)
v for NMP transitions between the charge states of
the defect. In order to calculate the NMP line shape of a defect using (4.31), we need
to find the vibrational wave functions 0I and +J and the associated energies. As
discussed in Sec. 2.9.2, for a quantum mechanical treatment of the vibronic transitions it is
necessary to approximate the potential energy surfaces of the different electronic states
as parabolas, see (2.85), which leads to a system of normal modes {Q_{1},Q_{2},Q_{3},…}. To
extract quantum mechanical line shapes from the density functional calculation we have to
approximate the DFT potential energy surfaces E0d() and E+
d () as quadratic functions.
Due to the inherent anharmonicity of these potential energy surfaces, different harmonic
parametrizations can be defined for the same defect, depending on which properties are
to be reproduced as accurately as possible [131]. In the present work, we have chosen a
parametrization that gives the same thermodynamic and switching defect levels as in the full density
functional calculation [86, 85]. Additionally, we assume here that only one mode couples to the
transition. While (2.102) is formulated upon an arbitrarily large set of coupling modes,
the exact modal spectrum obtained from the electronic structure method will generally
include some degree of mode mixing [131, 194], which requires a generalization of the
theories this work is based on. The single mode approach must be considered only a first
order approximation [144, 195]. However, the application of single effective modes is quite
common in the literature [131, 196], and it was shown quite recently that the approximation
obtained from our extraction scheme compares well to line shapes including multiple modes
[194].
In the following, we denote the nuclear optimum configuration of the neutral charge state with 0
and the nuclear optimum configuration of the positive charge state with +. Our single mode
approximation defines the direction vector from 0 to + as mode vector of the coupling mode. Based
on (2.56), we define
 (4.58) 
Consequently, the mass M associated with this mode is given by (2.58). For reasons of convenience, we
set the modal shift of the neutral state to zero, Q0 = 0, and define
 (4.59) 
The approximate potentials for the coupling mode read
E0
a(Q)  = E0
d(0) + Q^{2} and  (4.60)

E+
a (Q)  = E+
d (+) + (Q  Q′)^{2}.  (4.61) 
The switching trap levels are defined by the difference between potential energies at the optimum
positions. Thus, the parameters ω0 and ω+ are obtained from the conditions
 (4.62) 
which gives the oscillation frequencies as
ω0  = and  (4.63)

ω+  = .  (4.64) 
The resulting vibrational wave functions are harmonic oscillator wave functions (2.78)
 (4.65) 
where I,λ,Q is the Ith eigenvector of the harmonic oscillator of normalized frequency λ, displaced
by Q, see (2.78). The associated energies are
 (4.66) 
The line shapes of the defect are now calculated by inserting (4.65) and (4.66) into (4.31)
as
f(+∕0)
v (E)  = avg _{I} ∑
_{F }^{2}  

 δℏω0(F + 1∕2)  ℏω+(I + 1∕2) + E_{th}^{(+∕0)}  E,  (4.67)

f(0∕+)
v (E)  = avg _{I} ∑
_{F }^{2}  

 δℏω0(I + 1∕2)  ℏω+(F + 1∕2) + E_{th}^{(+∕0)}  E,  (4.68) 
with the thermodynamic transition level
 (4.69) 
as defined in Sec. 4.5. The FranckCondon factors F,Mω0,0I,Mω+,Q′ and
F,Mω+,Q′I,Mω0,0 are overlaps of harmonic oscillator wave functions [127]
F,Mω0,0I,Mω+,Q′  = ∫
_{∞}^{∞}h_{
Mω0}^{F }(Q)h_{
Mω+}^{I}(Q  Q′)dQ  (4.70)

F,Mω+,Q′I,Mω0,0  = ∫
_{∞}^{∞}h_{
Mω+}^{F }(Q  Q′)h_{
Mω0}^{I}(Q)dQ.  (4.71) 
Several analytic solutions for this integral can be found in the literature [197, 198, 199, 200]. In our
defect studies, we have compared several different methods. Our initial calculations employed the
expressions of Zapol [198], which give a direct solution for any combination of oscillation frequencies
and for any quantum numbers. However, for quantum numbers I and F larger than about 80, the
factorials in the overlap expressions are numerically problematic as they lead to overflow and roundoff
errors. Numerical integration of the harmonic oscillator wave functions allows the calculation of
overlaps at higher quantum numbers. However for quantum numbers larger than about 120 this
method suffers from the large roundoff errors in the summation. The best results were obtained using
the recurrence relations published recently by Schmidt [200], which are computationally
efficient and numerically stable up to very high quantum numbers [201], see Fig. 4.19.
 M∕amu  ω0∕s^{1}  ω+∕s^{1}   Q′∕Å 
     
Oxygen vacancy dimer  24.662  7.018 × 10^{13}  4.476 × 10^{13}  (36%)  0.504 
Oxygen vacancy puckered  25.644  5.974 × 10^{13}  7.931 × 10^{13}  (+33%)  0.411 
Hydrogen bridge closed  16.329  3.884 × 10^{13}  4.043 × 10^{13}  (+4%)  0.923 
Hydrogen bridge broken  17.928  2.227 × 10^{13}  2.736 × 10^{13}  (+23%)  1.700 
Table 4.5:  Parameters of the extracted potential energy surfaces. 
Fig. 4.20 illustrates the approximate potential energy surfaces for our atomistic models of the
oxygen vacancy and the hydrogen bridge. The parameters with respect to (4.60) and (4.61) are given
in Tab. 4.5. As explained in Sec. 4.4, the relative energetic position of these potentials
depends on the reservoir state E_{k} which is involved in the transition. The reservoir state
for the potentials in Fig. 4.20 is the silicon valence band edge, which corresponds to the
energetic situation for hole capture and emission at zero oxide field. A general feature of the
extracted potentials is the difference in oscillator frequency between the neutral and the
positive state. This is an important finding, as the expressions that are usually used to fit
experimental data are derived assuming linear coupling modes, i.e. ω_{i} = ω_{f} [180, 181, 202].
To check the validity of our line shape calculation method, we compare our numerically calculated
line shapes to the popular formula of Huang and Rhys [127]. The HuangRhys line shape is derived
for linear coupling, which means that the oscillation frequency before and after the transition are the
same, ω+ = ω0 = ω, see Sec. 2.9.2. The electron emission line shape for this case reads
[117]
f(0∕+)
v (E)  = ∑
_{p=∞}^{∞}W_{
p}δ(E_{th}^{(+∕0)}  ℏωp  E),  (4.72) 
with
W_{p}  = e^{S(2n+1)}^{p∕2}I_{
p}(2S,  

n  = , S = .   
Two results of this comparison are shown in Fig. 4.21 for the parameters extracted for
the oxygen vacancy dimer and the closed hydrogen bridge. Due to the limitation of the
HuangRhys formula, we have set ω+ = ω0 for both defects, to make a direct comparison possible.
Both methods give identical results, which we take as an indication of the validity of our
approach.
For the calculation of the line shapes we have considered 200 neutral and 200 positive states.
The line shape functions as calculated by (4.67) and (4.68) are a series of weighted Dirac
impulses, see Fig. 4.22. These Dirac peaks, however, are artifacts of the time dependent
perturbation theory that is used to calculate the rates. This approach assumes that at t = 0
the perturbation operator is “switched on” and then the system evolves freely for very
long time. This clearly is not the case here, as one can expect that the time range during
which the defect evolves without perturbation from the environment is rather short. For
shorter times the Dirac peaks are widened [98] leading to smoother resonances. Additional
spreading of the peaks arises from the nonorthogonality of the selected coupling mode,
the neglect of other coupling modes and the neglect of the energetic contribution of the
perturbation operator. Unfortunately, all these effects are impossible to model rigorously. In
the literature, it is common to account for this “life time broadening” by moving from
integral phonon numbers to fractional ones [127, 117]. Our results are corrected to give
continuous line shapes by smearing with a normal distribution of standard deviation kBT, see
Fig. 4.22. As long as the distributions drop sharper than the weights of the Dirac peaks, the
results are not very sensitive on the smearing parameter. The application of a temperature
dependent smearing parameter is also disputable, Alkauskas et al. have used a fixed smearing
value of 10meV [194]. A more rigorous description of the smearing requires experimental
input, possibly from transitions at low temperatures, that is not available at the moment.
The line shape functions for all transitions are shown in Fig. 4.23, which shows that stronger
bonds in the initial state lead to sharper line shape peaks. The bond of the oxygen vacancy
dimer is weakened by hole capture, so the f(+∕0)
v line shape is wider than the f(0∕+)
v line
shape. The electron in the occupied defect state of the hydrogen bridge has a fraction
of antibinding [59], so the bond that forms over the hydrogen atom is stronger when
this electron is removed and the f(+∕0)
v line shape is slightly sharper than the f(0∕+)
v line
shape.
Fig. 4.24 shows the temperature behavior of the (0∕+) and the (+∕0) line shape of the oxygen
vacancy dimer, both before and after the smearing. An increase in the temperature increases the
weight of the peaks that correspond to higher energy initial states. Especially in the unsmeared line
shape of the oxygen vacancy dimer, one can see a main comb of Dirac peaks that is only weakly
temperature dependent, which corresponds to the initial vibrational states of lower energy. The gaps
between the peaks of the main comb are filled by smaller peaks at higher temperatures. The
oscillations in the smeared line shapes for 100K are partly due to the gaps of the main comb not being
filled up and partly due to the temperature dependent smearing parameter used in this
work.
4.7 Extraction of the Classical Line Shapes
As discussed in Sec. 2.10, line shapes can also be calculated based on classical statistical
physics. In the defectdevice system, the shift induced by the mean field and the free state
enters directly into the energy conservation expression. Again we can define the defect line
shapes
f(+∕0)
v (E)  = ∫
_{ℝN}Z_{0}^{1}e^{
}δd,  (4.73)

f(0∕+)
v (E)  = ∫
_{ℝN}Z_{+}^{1}e^{
}δd,  (4.74) 
with the partition functions
Z_{0}  = ∫
_{ℝN}e^{
}d,  (4.75)

Z_{+}  = ∫
_{ℝN}e^{
}d.  (4.76) 
In these expressions, the integration runs over all nuclear degrees of freedom. For our 72 atom defect
structures, the direct numerical solution would require the evaluation of an integral in 216 dimensions,
which is clearly beyond our computational capacity. Nevertheless, there are atomistic modeling
methods for the evaluation of thermal integrals and we will come back to this below. At first, it is
interesting to see how strong the nuclear quantum effects are for the defects under investigation. For
this purpose, we calculate the classical line shapes for the approximate onedimensional potential
energy surfaces.
Inserting the one dimensional potentials (4.60) and (4.61) into (4.73) and (4.74) yields
f(+∕0)
v (E)  = ∫
_{∞}^{∞}Z_{
+}^{1}e^{
}δdQ′  (4.77)

f(0∕+)
v (E)  = ∫
_{∞}^{∞}Z_{
0}^{1}e^{
}δdQ,  (4.78) 
with the partition functions
Z_{+}  = ∫
_{∞}^{∞}e^{
}dQ =  (4.79)

Z_{0}  = ∫
_{∞}^{∞}e^{
}dQ = .  (4.80) 
For the solution of the integrals above we use
 (4.81) 
where the Q_{i} are the zeros of G(Q). The Dirac distributions thus reduce the line shape integrals to a
summation over the two points
Q_{1,2}(E) = ,   (4.82) 
which are exactly the crossing points of the two parabolas defined by (4.60) and (4.61). The final
expressions for the line shapes read
f(0∕+)
v (E)  = Z_{0}^{1} + 
(4.83)

f(+∕0)
v (E)  = Z_{+}^{1} + 
(4.84) 
These line shapes have the typical Arrhenius form of a thermally activated process. The activation
energies in this case are the crossing points of the potential energy surfaces [123]. The classical line
shapes are simple analytic expressions, which can be easily implemented into a device simulator. Due
to the Boltzmann terms in the classical line shapes, the crossing point that is lower in energy
dominates the transition. As discussed before, the crossing points themselves depend on the energy of
the free state.
The classical f(+∕0)
v line shapes are compared to their quantum mechanically calculated
counterparts in Fig. 4.25 for the oxygen vacancy dimer and the closed hydrogen bridge. The classical
formula underestimates the transition rate at low temperatures and for energies that are below the
maximum of the line shape. This energy range corresponds to the weak coupling regime of the
defect [37]. These underestimations are due to the neglect of tunneling in the classical model. For
strong coupling, which is especially relevant for the exchange of holes with the silicon valence band,
good agreement between the classical and the quantum mechanical version is already given at room
temperature. We take this result as an indication that the quantum mechanical nature of the nuclei
does not have a strong influence on the line shapes in the temperature ranges that are of interest to
BTI.
When quantum effects are negligible, however, it is also possible to calculate line shapes with the
full potential energy surfaces from the density functional calculations. Quite generally, the integrals of
a function X() of the form
 (4.85) 
are called “thermal averages” of the canconical (N,V,T) ensemble and are written as . As the
potential E() that determines the weighting factor for the averaging is usually a function of low
symmetry in ℝ^{N}, these integrals are impossible to evaluate from a simple numerical integration based
on discretization. The calculation methods that are applied to this type of problems are Metropolis
Monte Carlo and molecular dynamics. The Metropolis Monte Carlo method moves through
the configuration space on a random trajectory. This trajectory is steered such that the
regions with the highest occupation probability are most densely sampled. Denoting the
trajectory of the Metropolis Monte Carlo algorithm as _{1},_{2},_{3},…,_{I}, the average of the
function X() over the set of sample points will approach the true average of the canonical
ensemble
 (4.86) 
The Metropolis Monte Carlo algorithm is broadly applied in atomistic modeling to calculate thermally
influenced properties as bond lengths, free energies, etc. and its general formulation makes
it possible to include a particle reservoir and pressure in addition to a thermal reservoir
[120, 112, 203].
The molecular dynamics method on the other hand tries to describe the dynamics of the
molecular system by propagating the positions of the nuclei according to Newton’s equations of
motion [119, 120]. In its pure form, the method just generates the constant energy trajectory of an
atomic system from a given starting point in phase space [119], which corresponds to the
microcanonical (N,V,E) ensemble of statistical physics [112]. As this ensemble has very limited
practical relevance, there are extensions to the method for the coupling to a heat bath, which
results in a canonical (N,V,T) ensemble, and the inclusion of pressure, which results in an
isothermalisobaric (N,p,T) ensemble [120, 112]. These extensions are sometimes referred to as
thermostats and barostats, respectively. Similar to the Monte Carlo Method, the temperature and
pressure control are designed so as to give the correct thermal averages in the long time
limit,
 (4.87) 
(N,V,T) and (N,p,T) molecular dynamics simulations are routinely used to compute the dynamics of
frequent events such as the hopping of fast diffusors or reactions at high temperatures [153]. However,
it has to be mentioned here that strictly speaking the trajectories generated in those runs are not
‘physical’ trajectories, as the thermostats and barostats are only designed to give the right averages,
rather than nuclear dynamics.
The classical line shapes can also be written as thermal averages
f(+∕0)(E)  = δ(E_{
sw}()  E)_{+}  (4.88)

f(0∕+)(E)  = δ(E_{
sw}()  E)_{0},  (4.89) 
with
 (4.90) 
As mentioned in Sec. 2.10, these averages give the probability that for a given carrier energy E an
elastic capture or emission transition is possible. In the terminology of Sec. 4.5, the line shape can be
seen as a thermally broadened switching trap level E+
d ()  E0
d()  Ev0^{SiO2}. In a numerical
calculation, we can approximate this distribution as a normalized histogram over the sample points
drawn from a Metropolis Monte Carlo or molecular dynamics trajectory.
We have tested this calculation method for the closed hydrogen bridge. For the generation of the
trajectory, we use the molecular dynamics implementation of VASP, which implements a Nosé
thermostat [112]. Due to the large computational demand of the molecular dynamics calculation, the
accuracy of the density functional calculation is reduced here, using a plane wave cutoff of 500eV and
the real space projection feature of VASP. We calculate the f(0∕+)
v line shape from the molecular
dynamics trajectory in the neutral state of the defect as a histogram of the switching trap levels of all
configurations in the trajectory. The structure is equilibrated for 3ps, then the line shape is extracted
from a 10ps simulation.
The result is shown in Fig. 4.26. Excellent agreement is found between the line shapes calculated
from molecular dynamics and from the approximate potentials around the maximum and above. We
take this result as an indication for the validity of the harmonic approximation. The calculation of line
shapes from molecular dynamics or Metropolis Monte Carlo is undoubtedly an appealing
approach as it includes all features of the density functional potential energy surfaces.
Its applicability to the calculation of line shapes for BTI, however, is limited by the long
trajectories that have to be calculated in order to get smooth line shapes further away
from the maximum. The line shape in Fig. 4.26 is only smooth to about 0.7eV above the
maximum. As shown below, the hole capture in BTI happens at about 1.2eV above the
maximum, which due to the squaredexponential drop of the line shape we estimate to
require about a billion molecular dynamics steps. Unfortunately, this is again far beyond our
patience. As for the barrier hopping transitions, however, an approach similar to the nudged
elastic band method could be developed to find the lowest energy crossing point of the
two potential energy surface, in order to estimate the activation energy of the transition
process.
4.8 Density Functional Dependence
The PBE gradient corrected functional we employ in our work is quite popular in solid state theory.
However, adding a fraction of HartreeFock exchange to the density functional calculation improves
the description of band gaps and other properties of many atomic systems [192, 99]. We have
checked the dependence of the calculated line shapes on the employed functional by comparing our
PBE results with calculations using the PBE0 hybrid functional. Due to the high computational effort
that is required with hybrid functionals, we only investigated the primary states of our model
defects, i.e. the oxygen vacancy in its dimer configuration and the closed hydrogen bridge. A
comparison between the PBE and PBE0 prediction of the defect levels for the model defects
considered here has been done by Alkauskas et al. [204]. This comparison found that the
defect levels were largely shifted by about 1.45eV, but the relative positions were almost
the same. The findings of Alkauskas and coworkers agree well with our results shown in
Fig. 4.27, except for the strong shift which is absent in our work due to the marker based
alignment.
4.9 Energy Alignment
The biggest uncertainty in our calculations comes from the alignment of the energy scales between the
atomistic defect model and the free charge carrier states in the semiconductor device. Although the
marker method performs well for defects that are sufficiently similar to the marker defect [102], the
accuracy of the alignment depends heavily on the accuracy of the experimental reference values.
Blöchl’s placement of the (+∕) transition level 0.2eV above the silicon midgap is based on CV
measurements of MOS structures after hydrogen exposure [59, 81]. However, the distribution of
the states which have been measured in [81] spans from 0.1eV below silicon midgap to
0.3eV above, and the assignment of this transition level to the hydrogen interstitial is
debatable. As an alternative one can take the band edges of the KohnSham system as reference
levels. This leads to large differences in our calculations as the SiO2 valence band edge
of the KohnSham system lies 1.53eV above the level estimated from Blöchl’s alignment
scheme. Interestingly, this difference is reduced to 149meV for the PBE0 functional. Another
uncertainty concerns the SiSiO2 valence band offset. Different values have been given in the
literature, ranging from 4.2eV to 4.8eV [205]. Considering all these uncertainties, in a future
comparison of model defects to experimental data it may thus be favorable to leave some
freedom to optimize the energy alignment between the energy scale of the line shapes and the
device.
4.10 Discussion of the DFT Results in the Context of BTI
In the multistate multiphonon model for BTI, the more stable configurations in the neutral state and
the positive state are 1 and 2, the metastable configurations are 1′ and 2′, respectively. Based on the
energy differences of the configurations given in Sec. 4.1.2, the configurations 1 and 2′ correspond to
the neutral and positive state of the dimer configuration for the oxygen bridge and the closed
configuration of the hydrogen bridge. Consequently, the puckered/broken configurations are assigned
1′ in the neutral state and 2 in the positive state.
 1 → 1′  1′→ 1  2 → 2′  2′→ 2 
    
Oxygen vacancy  3.00eV  36.0meV  0.43eV  0.34eV 
Hydrogen bridge  1.18eV  1.16eV  0.46eV  0.18eV 
BTI defect [206]  ≈ 1.00eV  ≈ 1.2eV  ≈ 0.7eV  ≈ 0.3eV 
Table 4.6:  Activation energies for the structural reconfigurations as calculated from the
atomistic defect models using the nudged elastic band method. 
The reaction paths for the 1 ⇌ 1′ and 2 ⇌ 2′ transitions are shown in Fig. 4.8, the extracted
barriers are compared to estimates for the BTI defect in Tab. 4.6. The hydrogen bridge shows
reasonable agreement with the BTI defect values both in the neutral and the positive state. For the
oxygen vacancy, the 2 ⇌ 2′ barriers are also close to the experimental estimates, but a strong
discrepancy arises for the 1 ⇌ 1′ barriers. The barrier of about 36meV for the transition from the
puckered to the dimer configuration is much too small for BTI, as the experiments show defects that
can be switched several times between the positive and the neutral state in the secondary
configuration [37]. A 1′→ 1 barrier as small as the oxygen vacancy value would be overcome
instantaneously after neutralization of the defect. However, this issue may be a feature of the
crystalline host structure, as higher barriers have been reported [64] for the oxygen vacancy in
amorphous silica.
For the discussion of the charging and discharging reactions 1 ⇌ 2′ and 2 ⇌ 1′, we first take a look
at the f(0∕+)
v and f(+∕0)
v line shapes in Fig. 4.23. One can see a striking difference in the hole capture
(electron emission, f(0∕+)
v ) behavior of the dimer and the puckered state of the oxygen vacancy, while
this difference is much less pronounced for the closed and broken hydrogen bridge. The
position of the f(0∕+)
v maximum of the oxygen vacancy dimer is close to the SiO2 valence
level, in agreement with the switching trap level for that transition given in Tab. 4.4. A
f(0∕+)
v line shape at this energetic position makes a hole capture from the Si valence level very
unlikely, even under large bias. Gös et al. have calculated switching trap levels for oxygen
vacancies in amorphous silica [193] and have found a very small spread in the distribution of
these levels for the investigated samples. One can thus assume that the low position of the
hole capture line shape is only weakly influenced by the surrounding atoms. This makes
the oxygen vacancy dimer position a very unlikely candidate for the state 1 of the BTI
defect.
The maximum of the hole capture line shape of the closed hydrogen bridge on the other hand is too
close to the silicon valence band to show the BTI behavior. The widening due to an increase in the
temperature as well as the bias dependence are stronger for carriers further away from the line shape
maximum. In the shown configuration, both the bias and the temperature dependence are weak, which
unfortunately does not match the properties of the BTI defect. However, the calculations in [193]
showed a spread for the (0∕+) switching trap level of the hydrogen bridge in amorphous
silica that is larger compared to the same level of the oxygen vacancy. Thus, there may be
hydrogen bridges which match the hole capture properties of the BTI defect in amorphous
silica.
While the line shape of a defect includes the quantum mechanical effects, and the dependence on
the reservoir state energy and the band bending, it does not give any information about
the temperature behavior. However, as mentioned in Sec. 4.7, the activation barrier for
the charge state transition can be estimated from the intersection of the potential energy
surfaces. The parabolic approximations extracted from our density functional calculations are
shown in Fig. 4.20, where the energy of the involved free state is placed at the silicon
valence level. This corresponds to the most dominant energetic configuration for hole capture
with no field in the oxide. It can be seen that in the primary (dimer) state, the oxygen
vacancy has to overcome a large thermal barrier to capture a hole from the silicon valence
band edge, which corresponds to the small value of the line shape at this energy. The
temperature activation of the 1 → 2′ transition determines the temperature activation of the BTI
degradation, which is in the 0.2 – 0.3eV regime. The initial hole capture barrier of the
oxygen vacancy is ≈ 2.7eV at zero bias, and a reduction to the experimental BTI activation
energy requires a field induced shift of ≈ 2.6eV. At a typical NBTI field of 7MV∕cm, this
requires a distance of about 3.7nm between the defect and the interface, which is beyond the
thickness of most stateoftheart MOS oxides. The hydrogen bridge on the other hand
shows a transition barrier of ≈ 35meV for the initial hole capture. This value is too low
compared to experimental BTI data, considering that this value is without any applied
field.
Only minor changes to the line shapes are observed when moving from the PBE to the PBE0
functional, except for the (0∕+) line shape of the oxygen vacancy in the dimer configuration. This
line shape is moved further down on the energy scale, which makes a hole capture from
the silicon valence band even more unlikely. The predicted charge transition properties
depend strongly on the selected alignment scheme. Using the KohnSham valence level of
our density functional calculations as a reference instead of the (+∕) transition level of
the hydrogen interstitial shifts all trap levels and line shapes down by more than 1.5eV.
The line shape of the oxygen vacancy in this case is moved into the SiO2 valence band,
making its hole capture behavior even less compatible with the BTI defect. The hydrogen
bridge, on the other hand, assumes an energetic configuration that corresponds well with the
BTI experiments. Judging from the good agreement between the SiO2 valence level in
Blöchl’s marker method and the prediction of the PBE0 functional, the PBE valence level
does not seem to be a reasonable reference. Nevertheless, in the next section we take the
line shapes referenced to the PBE valence level as they give better results for the hole
capture rate in this case. The rate calculations presented here thus have to be considered a
proofofconcept of the calculation method rather than a presentation of rigorously calculated
results.
4.11 Hole Capture Rates
Once the line shapes are calculated and an accurate device simulation method is
selected, one can proceed to calculate capture and emission rates using (4.41)–(4.44). In
the present work, the carrier concentration of the MOS structure has been calculated
selfconsistently using the nonequilibrium Green’s function method [207]. The densities n(,E),
p(,E), D_{n}(,E) and D_{p}(,E) are readily obtained from the nonequilibrium Green’s
functions [208] .
The formalism assumes thermal equilibrium in the gate and bulk region where level broadening due to
scattering is modeled using an optical potential [207]. The oxide is treated as a nonequilibrium
domain with ballistic quantum transport. Fig. 4.28 shows the hole density in the band
diagram as obtained from the NEGF method. The device calculations are done by the Vienna
SchrödingerPoisson software package (VSP2) [179]. The device structure consists of a polySi gate
and an ndoped bulk separated by a 2nm SiO2 layer. For electrons the unprimed and primed valleys
with 0.19m_{e} and 0.91m_{e} electron mass are included. Holes were considered with 0.49m_{e} effective
mass.
We have implemented the rate calculation directly into the VSP2 simulator. The calculation of
the NMP hole capture rates proceeds in a two step process. First, the band bending is
calculated by solving the Poisson and the NEGF equations selfconsistently. Secondly,
the NEGF problem is again solved nonselfconsistently on a different energy grid that
accounts for high energy holes as these contribute considerably to the NMP transitions. The
integration is implemented as a postprocessing step using the numerical NEGF and line shape
data.
At the time this thesis is written, the transition rate calculations involving unoccupied states in the
device, corresponding to hole emission and excess electron emission, is still work in progress. Of the
transitions involving occupied states, corresponding to carrier capture, the hole capture transition is
most important for the present work. For this reason we concentrate on the calculation of hole capture
rates using (4.42). Due to the empirical factor λ in (4.42), all capture rates can be calculated up to a
constant factor, therefore all computed time constants in this section are given in arbitrary
units.
First, we investigate the effect of the bias induced shift of the relative positions of the line shape and
the hole states. As mentioned in Sec. 4.4, this shift is especially relevant for the trapping behavior of
oxide defects. Application of a gate voltage induces the largest electric field, and in consequence the
largest band bending, in the oxide due to the absence of native carriers there. As illustrated in
Fig. 4.29, the dependence of the line shape on the value of the reference energy at the defect site plays
a crucial role for the transition kinetics, as the band bending energetically shifts the relative position
of the line shape and the spectrum of the hole states, leading to large changes in the capture
rate.
In semiconductor device modeling, defects are primarily considered as recombination centers and
described using the ShockleyReadHall (SRH) theory [209], which characterizes a defect via the
capture cross section σ and the (thermodynamic) trap level E_{T}. The capture of particles
is modeled as a flux through the cross section at the speed v_{th}, which is the “thermal
velocity”
k_{n→d}  = σv_{th}n(_{d}).  (4.91) 
The emission of an electron into a state of energy E is calculated based on the occupancy in thermal
equilibrium as
k_{d→E}  = k_{n→d}e^{
}.  (4.92) 
The temperature and bias dependence of the capture transition in the classic SRH model comes solely
from the density of carriers at the defect site. NMP capture can be modeled within the framework of
the SRH theory as an energy dependent capture cross section σ(E) [117]. We compare the capture
rates calculated from the NMP line shapes to SRHtype rates that only depend on the density of holes
at the defect site to illustrate the impact of the vibrational influence on the transition
rate.
Fig. 4.30 left shows the gate voltage dependence of the NMP hole capture time constants for the
two model defects. The gate voltage dependence predicted from the NMP theory is much stronger
than the prediction of th SRHlike model. Further, the gate voltage dependence increases with
the distance of the defect from the silicon bulk. Both of these strong dependencies are
caused primarily by the energetic shift of the line shape functions relative to the holes in
the inversion layer, as illustrated in Fig. 4.29. This behavior is in qualitative agreement
with the capture rates observed in TDDS measurements. As explained in the previous
section, the agreement would not be as good in the originally selected energy alignment
scheme.
The temperature dependence of the hole capture is shown in Fig. 4.31. The calculation of
NMP hole capture rates for the given defects becomes numerically challenging for low
temperatures. As the temperature decreases, the line shapes become increasingly narrow,
thus making high energy holes the dominantly captured particles. Accurate representation
of high energy holes in the NEGF algorithm requires an improved refinement strategy
for the energy grid, which is currently under development. To overcome this limitation,
the hole distribution over energy was calculated from a classical density of states for the
Arrhenius plot, taking only the total hole concentration at the defect site from the NEGF
calculation.
NMP defects show a strong temperature activation, in contrast to capture rates whose temperature
dependence comes from the total carrier density alone, as in the SRH description. The
calculated activation energies for the hole capture are in excellent agreement with the activation
energies obtained from TDDS experiments in [32], which is again strongly influenced
by the selected alignment scheme. The difference in temperature activation between line
shapes calculated using the quantum mechanical formula of (4.68) and those calculated
based on classical statistical physics (4.84) are also compared in Fig. 4.31. For the closed
hydrogen bridge, this difference becomes visible only below 140K. For the puckered oxygen
vacancy, the classical formula reproduces the quantum mechanical behavior over the complete
temperature range investigated. Again, we take this as an indication that the quantum
mechanical effects of the nuclei are negligible for the calculation of rates of typical BTI
experiments.
4.12 Related Work
NMP theory has been used by several authors in the context of semiconductor devices
[144, 180, 181, 202, 32, 182]. The transition rate formulas employed are usually based on linear
electronphonon coupling. As shown in Fig. 4.20, this assumption is not fulfilled for our density
functional based defect models. Purely quadratic coupling has been investigated in the literature
[128, 129], but a linearquadratic coupling that changes both the equilibrium position and the
frequency of the coupling mode has not been considered up to now.
With respect to the modal spectrum, it is either assumed that the transition couples to an infinite
number of phonon modes — all having the same oscillator strength — where each can contribute only
one phonon [127, 202, 181], or that the transition couples to only one effective mode which receives
or emits an arbitrary number of phonons [144, 117, 131]. Interestingly, for linear coupling modes
both assumptions lead to essentially the same expression for the capture rates. For the situation we
find in the bias temperature instability, in our opinion it is more reasonable to assume that
the NMP kinetics are determined by a small number of local modes at the defect site. A
coupling to a large number of modes does not seem reasonable for the defects involved in BTI
and RTN considering the large variations in transition rates between the defects that are
observed in measurements. These variations can only be explained by differences in the local
environment of the defect structure which can only hold a small number of vibrational
modes.