4  Atomistic Modeling and the BTI Defect

The defect in the multi-state multi-phonon 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 pseudo-molecules 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 quantum-chemistry method. This quantum region is surrounded by the so-called classical region, which contains a large number of atoms whose interactions are described using empirical potentials [66]. Finally, the quantum-mechanical 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 pseudo-potential 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 set-up of an atomistic defect model using the embedded cluster approach requires a lot of fine-tuning of both the empirical potential and the pseudo-potential 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.

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 gradient-corrected functional of Perdew, Burke, and Ernzerhof [101], the core electrons are represented using the projector augmented wave method as implemented in the Vienna ab-initio simulation program (VASP)  [169, 170]. The electronic structure is expanded using a plane-wave 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 non-orthogonal unit cell, see Fig. 4.2.

Figure 4.2:

The α-quartz unit cell given in  [168]. The cell contains nine atoms and the cell angles are α = β = 90∘ and γ = 120∘

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.

Figure 4.3:

Construction of an orthogonalized α-quartz supercell. The unit cell (i) is repeated along axes and to give a 2 × 2 supercell (ii). Due to the geometry of the cell it is possible to change the cell shape as indicated by the red rectangle (iii). Finally, the internals of the cell are adjusted to fit the new shape (iv).

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.

Figure 4.4:

The orthogonalized 72 atom α-quartz supercell used in our calculations.

The energy cut-off 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 quasi-Newton algorithm [171, 96] until the forces are below 10^-2eV∕Å. Optimizations starting far from the optimum were pre-optimized using the more robust conjugate gradient algorithm.

Figure 4.5:

The Kohn-Sham eigenvalue spectrum of the α-quartz host cell. Just as the LDA, the GGA underestimates band gaps. The band gap in the present work amounts to 5.9 eV.

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 Kohn-Sham 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.

Figure 4.6:

Structures of the dimer (left) and the puckered (right) configuration of the oxygen vacancy as obtained from our DFT calculations. The neutral (top) as well as the positive (bottom) charge state are shown.

Figure 4.7:

Structures of the closed (left) and the broken (right) hydrogen bridge as obtained from our DFT calculations. The neutral (top) as well as the positive (bottom) charge state are shown.

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

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 pre-optimized using the conjugate gradient algorithm until forces are below 1eV∕Å. The final optimization to reduce the forces below 0.1eV∕Å is done using the quasi-Newton algorithm. The resulting energies on the optimized reaction paths for both defects are shown in Fig. 4.8.

4.3  Macroscopic Device Simulation

Figure 4.9:

In a band diagram, the local values of the conduction and valence band edges are drawn as a function of position. At a material boundary, the band edges abruptly change giving rise to energetic barriers or wells. The shown band diagram corresponds to the situation in an MOS structure.

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 in-depth 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 non-ground-state 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 semiconductor-oxide 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

Figure 4.10:

The mean potential Φ(), which models the coulomb interactions within the semiconductor device, can be described as a bending of the band structure. The band bending as obtained from device simulation is shown for three different gate voltages (from left to right 0V, -1V, -2V). The bending at zero gate voltage arises from the work function difference between the gate and the substrate.

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 state-of-the-art TCAD simulators, are based on semi-classical electrons and holes. These semi-classical 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 self-consistently with the Poisson equation  [179], hence these carrier models are sometimes called Schrödinger-Poisson models. Due to the assumed equilibrium occupancy of the particle states these models give good results for non-conducting 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ödinger-Poisson models have been employed in calculations of defect-device 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ödinger-Poisson calculations have to be placed either at the oxide-gate or the oxide-substrate 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 non-equilibrium 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]

Dn(x,E)	= ∫ wn2δ(E - E wn′)dwn′	(4.8)
Dp(x,E)	= ∫ wp2δ(E - E wp′)dwn′	(4.9)

4.4  Non-Radiative Transitions in the Device-Defect System

For the calculation of non-radiative 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 device-defect 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 Born-Oppenheimer states, which are mixed by the effective non-adiabacity 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 non-adiabacity operator mixes the electronic states. In the device-defect 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 non-adiabacity 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 substrate-oxide 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 first-order. 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 self-consistent calculation.

Figure 4.11:

The energy of the reservoir state from which the electron is captured as well as the mean potential in the device simulation contribute to the energetic shift between the potential energy surface of the initial and the final state. Thus, every reservoir state has a different configuration of the potential energy surfaces as indicated by the different red curves. Additionally, a change in the mean potential from the device also leads to an energetic shift, which is indicated as the blue dashed curve.

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′fkEk[Φ] + E+ d () + Ek′[Φ] + q0Φ(d))iI	= EiIiI ,	(4.22)
(N + ∑ k≠k′fkEk[Φ] + E0 d())fF	= EfF fF .	(4.23)

(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′fkEk[Φ] + Ek′[Φ] + q0Φ(d)	(4.26)
EfF	= E0F + ∑ k≠k′fkEk[Φ],	(4.27)

(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 Frank-Condon 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 non-radiative 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.

Figure 4.12:

The line shape is referenced to the local band edges. The figure shows the position of a line shape (yellow) relative to the valence band edge (red) and the conduction band edge (blue). Different gate bias voltages lead to different band bendings and different positions of the capture line shape relative to the free quasi-particle states. In the example, the negative bias shifts the line shape maximum closer to the states near the silicon valence band edge, which is the typical NBTI situation.

An analogous derivation can be done for the emission of an electron. In this case the electronic states are

i	= (ĉ† d + ∑ k≠k′fkĉ† k)∅	(4.33)
f	= (ĉ† k′ + ∑ k≠k′fkĉ† k)∅.	(4.34)

EiI	= E0I + ∑ k≠k′fkEk[Φ] and	(4.35)
EfF	= E+F + ∑ k≠k′fkEk[Φ] + Ek′[Φ] + q0Φ(d),	(4.36)

(4.37)

The resulting expressions for the electron trapping and detrapping rates with respect to a state k′ are

Wk′→d	= λ2f(+∕0) v (Ek′[Φ] - Ev(d)),	(4.38)
Wd→k′	= λ2f(0∕+) v (Ek′[Φ] - Ev(d)),	(4.39)

(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 quasi-particle 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 ² 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

kn→d	= λ∫ -∞∞n( d,E)f(+∕0) v (E - Ev(d))dE,	(4.41)
kp→d	= λ∫ -∞∞p( d,E)f(0∕+) v (E - Ev(d))dE,	(4.42)
kd→n	= λ∫ -∞∞(D n(d,E) - n(d,E))f(0∕+) v (E - Ev(d))dE,	(4.43)
kd→p	= λ∫ -∞∞(D p(d,E) - p(d,E))f(+∕0) v (E - Ev(d))dE,	(4.44)

(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:

1. Calculation of the line shape functions f(+∕0) and f(0∕+) from the atomistic model.
2. Calculation of Ev(), n(,E), D_n(,E), p(,E), and D_p(,E) form the macroscopic device model.
3. 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 NMP-transition 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 Born-Oppenheimer energy Ed^q(_q⁰) at the optimum configuration of the defect _q⁰ 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.

Figure 4.13:

Charge state formation energies versus electron chemical potential for a hypothetical defect that has five charge states ++,+,0,-, and --. The crossing points of the formation energy lines corresponding to different charge states are the thermodynamic transition levels.

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₁ to charge state q₂ 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 negative-U defects as illustrated in Fig. 4.15. In SiO2, a negative-U 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 negative-U 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 non-equilibrium behavior of the defect.

Figure 4.15:

Illustration of a hypothetical negative-U defect. The neutral state of the defect is less stable than the charged states for all values of the chemical potential. At μe = 0eV, the negative charge state becomes more stable than the positive charge state, which gives rise to a E+∕- transition level.

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 Franck-Condon 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 Born-Oppenheimer 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^q₁∕q₂ from the condition E_sw^q₁∕q₂ = 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⁰) 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.

Figure 4.17:

The relation of the switching and thermodynamic defect levels to the Born-Oppenheimer potential energy surfaces of an atomic model of a point defect. It is assumed that the defect structure is neutral with N electrons.

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 well-defined 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 so-called 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 Si-SiO2 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.

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 Hartree-Fock 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 Hartree-Fock 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 Kohn-Sham 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 Kohn-Sham 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₁,Q₂,Q₃,…}. 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) + Q2  and	(4.60)
E+ a (Q)	= E+ d (+) + (Q - Q′)2.	(4.61)

(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 I-th 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) + Eth(+∕0) - E,	(4.67)
f(0∕+) v (E)	= avg I ∑ F 2
	δℏω0(I + 1∕2) - ℏω+(F + 1∕2) + Eth(+∕0) - E,	(4.68)

(4.69)

as defined in Sec. 4.5. The Franck-Condon factors F,Mω0,0|I,Mω+,Q′ and
F,Mω+,Q′|I,Mω0,0 are overlaps of harmonic oscillator wave functions [127]

F,Mω0,0|I,Mω+,Q′	= ∫ -∞∞h Mω0F (Q)h Mω+I(Q - Q′)dQ	(4.70)
F,Mω+,Q′|I,Mω0,0	= ∫ -∞∞h Mω+F (Q - Q′)h Mω0I(Q)dQ.	(4.71)

Figure 4.19:

Comparison of the harmonic oscillator overlaps calculated by numerical integration to the recurrence relations of Schmidt  [200]. For large quantum numbers the roundoff errors in the numerical integration become dominant, while the recurrence relations are numerically stable.

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 Huang-Rhys 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]

Wp	= e-S(2n+1)p∕2I p(2S,
n	= , S = .

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 non-orthogonality 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 anti-binding  [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.

Figure 4.24:

Quantum mechanical line shapes of the oxygen vacancy dimer both before (top) and after (bottom) smearing for different temperatures. With increasing temperature, the probability of finding the vibrational system in energetically higher states increases, increasing the weight of the corresponding Dirac-peaks. In the smeared line shape this shows as a broadening. For low temperatures, the smeared line shapes show oscillations.

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 defect-device 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)	= ∫ ℝNZ0-1e- δd,	(4.73)
f(0∕+) v (E)	= ∫ ℝNZ+-1e- δd,	(4.74)

Z0	= ∫ ℝNe- d,	(4.75)
Z+	= ∫ ℝNe- d.	(4.76)

Inserting the one dimensional potentials (4.60) and (4.61) into (4.73) and (4.74) yields

f(+∕0) v (E)	= ∫ -∞∞Z +-1e- δdQ′	(4.77)
f(0∕+) v (E)	= ∫ -∞∞Z 0-1e- δdQ,	(4.78)

Z+	= ∫ -∞∞e- dQ =	(4.79)
Z0	= ∫ -∞∞e- dQ = .	(4.80)

(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

Q1,2(E) = ,

(4.82)

f(0∕+) v (E)	= Z0-1 +	(4.83)
f(+∕0) v (E)	= Z+-1 +	(4.84)

Oxygen vacancy dimer

Hydrogen bridge closed

Figure 4.25:

Comparison of the classical line shapes from (4.84) to their quantum mechanical counterparts that were calculated numerically from (4.68) at the indicated temperatures for the dimer configuration of the oxygen vacancy and the closed hydrogen bridge. Deviations between the two versions arise at low temperatures and in the weak coupling regime (energies below the peak of the line shape) due to the absence of tunneling in the classical version.

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 ₁,₂,₃,…,_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 isothermal-isobaric (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)

(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 cut-off 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.

Figure 4.26:

Comparison of the line shapes for the closed hydrogen bridge as calculated from the approximate harmonic potential energy surface in the classical and quantum mechanical form to the line shapes calculated from molecular dynamics. Excellent agreement between the different calculation methods is found around the maximum of the line shape.

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 squared-exponential 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

Figure 4.27:

Comparison of the line shapes extracted from density functional theory based on the PBE functional and the PBE0 functional. Except for the (0∕+) line shape of the oxygen vacancy in the dimer configuration, the line shapes calculated using PBE0 agree well with the PBE results.

The PBE gradient corrected functional we employ in our work is quite popular in solid state theory. However, adding a fraction of Hartree-Fock 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 Kohn-Sham system as reference levels. This leads to large differences in our calculations as the SiO2 valence band edge of the Kohn-Sham 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 Si-SiO2 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 multi-state multi-phonon 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.

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 state-of-the-art 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 Kohn-Sham 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 proof-of-concept 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 self-consistently using the non-equilibrium Green’s function method  [207]. The densities n(,E), p(,E), D_n(,E) and D_p(,E) are readily obtained from the non-equilibrium 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 non-equilibrium 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ödinger-Poisson software package (VSP2)  [179]. The device structure consists of a poly-Si gate and an n-doped 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 self-consistently. Secondly, the NEGF problem is again solved non-self-consistently 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 post-processing 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.

Figure 4.29:

Illustration of the bias induced shift of the relative position of the line shape and the free hole states. The f(0∕+) v line shape of the closed hydrogen bridge is used, the defect position is 2Å from the Si-SiO2 interface. With more negative gate voltage, the overlap between f(0∕+) v (E - Ev(d)) and h(d,E), and in consequence the transition rate, strongly increase.

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 Shockley-Read-Hall (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”

kn→d

= σvthn(d).

(4.91)

kd→E

= kn→de- .

(4.92)

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 SRH-like 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.

Figure 4.31:

Arrhenius plots of the capture rates for the defect types comparing the quantum mechanical NMP capture rates (4.68) for the puckered oxygen vacancy and the closed hydrogen bridge to the capture rates calculated using classical atoms (4.84) and to a SRH like model. The indicated activation energies are obtained by fitting to an Arrhenius law.

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 electron-phonon 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 linear-quadratic 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.