Thesis

Atomistic Modelling of Charge Trapping Defects
in Silicon Dioxide
Al-Moatasem Bellah El-Sayed
Department of Physics and Astronomy
University Collge London
2015
A thesis submitted in partial fulfilment of the requirements for the
degree of Doctor of Philosophy (Ph.D.) of the University College
London
Supervisor: Prof. Alexander L. Shluger

I, Al-Moatasem Bellah El-Sayed, confirm that the work presented in this thesis is
my own. Where information has been derived from other sources, I confirm that this
has been indicated in the thesis.
Al-Moatasem El-Sayed
June 2015
Abstract
This thesis focuses on the atomistic modelling of electron trapping sites in silicon
dioxide (SiO2 ) and interactions of hydrogen with the SiO2 network and the resulting
defects they produce. They are discussed in the context of electronic device reliability
issues.
The results presented here were calculated in both crystalline (c-) and amorphous
(a-) SiO2 . Due to its disordered nature, modelling a-SiO2 is challenging and required
the use of a melt-and-quench technique using a classical interatomic potential. All
models were evaluated against experimental data to ensure that they are indeed rep-
resentative of a-SiO2 .
Using density functional theory (DFT) and the models described above, extra
electrons were shown to trap in pure c- and a- SiO2 in deep band gap states for the first
time. They can trap spontaneously on pre-existing structural precursors in a-SiO2 .
The optical absorption spectrum of the intrinsic electron traps was calculated using
time-dependent DFT and shows a peak at 3.7 eV, which is in good agreement with a
previously unidentified experimental absorption peak measured at low temperatures.
Electronic device fabrication processes widely employ hydrogen due to its perceived
ability to improve their reliability. The results in this thesis show that both atomic
and molecular hydrogen are involved in defect generation processes. Atomic hydrogen
was found to interact strongly with strained Si-O bonds to form a stable defect. The
barriers to create and annihilate this defect as well as the distribution of its properties
were calculated. Hydrogen molecules were found to generate Si-H bonds which can
trap holes to form a similar defect and release a proton which can modulate the
defect’s properties.
These results provide insight into the atomistic nature of defects that can be
involved in electronic device reliability issues and help guide the design of reliable
fabrication processes.
Contents
1 Introduction 1
1.1 Background & Motivation . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Electronic Device Reliability Issues . . . . . . . . . . . . . . . . . . . 3
1.3 Atomistic Nature of Defects Implicated in Reliability Issues . . . . . . 6
1.4 Scope of this Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Theoretical Background and Methods 11

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.1 Gaussian Plane Waves Method . . . . . . . . . . . . . . . . . 18
2.3.2 Hybrid Functionals and the Auxiliary Density Matrix Method 21
2.4 Classical potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4.1 Charge variable potentials . . . . . . . . . . . . . . . . . . . . 24
2.4.2 Charge-optimized many-body potential (COMB) . . . . . . . 26
2.4.3 ReaxFF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.5 Embedded Cluster Method . . . . . . . . . . . . . . . . . . . . . . . . 30
2.6 Nudged Elastic Band Method . . . . . . . . . . . . . . . . . . . . . . 33
3 Modelling Amorphous Silicon Dioxide and the Si/SiO2 interface 37

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Choosing the Classical Potential . . . . . . . . . . . . . . . . . . . . . 39
3.2.1 Small Silicon and Silica Test Clusters . . . . . . . . . . . . . . 40
3.2.2 NVE Molecular Dynamics using ReaxFF and COMB . . . . . 44
3.2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.3 Making Amorphous Silica Models . . . . . . . . . . . . . . . . . . . . 47
3.3.1 Modelling a-SiO2 using the ReaxFF potential and the Melt-and-
Quench Method . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.2 DFT Optimisation of ReaxFF Structures . . . . . . . . . . . . 53
3.3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.4 Constructing the Si/SiO2 Interface . . . . . . . . . . . . . . . . . . . 58
3.4.1 Melt and Quench in situ . . . . . . . . . . . . . . . . . . . . . 59
3.4.2 a-SiO2 Layer Strained to Silicon Lateral Dimensions . . . . . . 61
3.4.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4 Defects in Bulk SiO2 66

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2 Quartz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2.1 Neutral Oxygen Vacancy . . . . . . . . . . . . . . . . . . . . . 67
4.2.2 Positive Oxygen Vacancy . . . . . . . . . . . . . . . . . . . . . 71
4.2.3 Negative Oxygen Vacancy . . . . . . . . . . . . . . . . . . . . 78
4.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3 Amorphous Silicon Dioxide . . . . . . . . . . . . . . . . . . . . . . . . 80
4.3.1 Neutral Oxygen Vacancy . . . . . . . . . . . . . . . . . . . . . 80
4.3.2 Positive Oxygen Vacancy . . . . . . . . . . . . . . . . . . . . . 82
4.3.3 Singly Negative Oxygen Vacancy . . . . . . . . . . . . . . . . 83
4.3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5 Electron Trapping at Impurities in Silicon Dioxide 88
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.2 Calculation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.3 Results of calculations . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3.1 Ge-doped α-quartz . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3.2 Li-doped α-quartz . . . . . . . . . . . . . . . . . . . . . . . . 96
5.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 100
6 Intrinsic Electron Trapping Sites in Silicon Dioxide 102

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.2.1 Classical Calculations . . . . . . . . . . . . . . . . . . . . . . . 104
6.2.2 Density Functional Theory Calculations . . . . . . . . . . . . 106
6.3 Results of Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.3.1 Electron Trapping in Bulk α-Quartz . . . . . . . . . . . . . . 107
6.3.2 Electron trapping in amorphous SiO2 . . . . . . . . . . . . . . 108
6.3.3 Concentration of electron trapping sites in a-SiO2 . . . . . . . 112
7 Optical Absorption Spectra of Trapped Electrons in SiO2 117

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.1.1 Experimental Observations . . . . . . . . . . . . . . . . . . . . 118
7.3 Results of Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7.3.1 Optical Absorption of the Ge Electron Trapping Centre . . . . 123
7.3.2 Optical Absorption of the Electron Trapping Centre . . . . . . 125
8 Hydrogen’s Interactions with Silicon Dioxide 129

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.3 Results of calculations . . . . . . . . . . . . . . . . . . . . . . . . . . 133
8.3.1 Atomic Hydrogen in a-SiO2 . . . . . . . . . . . . . . . . . . . 133
8.3.2 Molecular hydrogen in a-SiO2 . . . . . . . . . . . . . . . . . . 141
9 Concluding Remarks 151

9.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
9.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
A BKS Parameterisation for Si and O 155
B Formation Energy Calculation 156

B.1 Charge Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
List of Figures
1.1 A transmission electron micrograph (TEM) image of a silicon MOSFET. 3

1.2 Typical structure of a MOSFET. . . . . . . . . . . . . . . . . . . . . 4
1.3 Defect assisted tunnelling in electronic devices. . . . . . . . . . . . . . 5
1.4 Atomistic structure of the hydrogen bridge. . . . . . . . . . . . . . . . 7
2.1 Silicon-silicon bond-order dependence on interatomic distance. . . . . 29

2.2 Schematic of embedded cluster method. . . . . . . . . . . . . . . . . . 32
3.1 Two interlinked SiO4 units illustrating the flexible Si-O-Si angle. . . . 38
3.2 Atomic structures of Si3 Oy clusters . . . . . . . . . . . . . . . . . . . 40
3.3 Atomic structures of Si3 O6 optimized using COMB. . . . . . . . . . . 42
3.4 Atomic structures of Sin On clusters . . . . . . . . . . . . . . . . . . . 42
3.5 Atomic structures of Sin clusters, where n=2-6. . . . . . . . . . . . . 43
3.6 Plots of total energy from MD simulations using ReaxFF and COMB 46
3.7 Histograms of the short-range geometrical properties of a-SiO2 from
320 a-SiO2 models generated using the ReaxFF potential. . . . . . . . 49
3.8 Real space representation and spatial frequencies of crystalline and
amorphous solids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.9 Total neutron structure factors from a-SiO2 models and from neutron
scattering experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.10 Histograms of the short-range geometrical properties of 320 DFT opti-
mised models of a-SiO2 . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.11 Histogram of the atomic displacements of the ReaxFF a-SiO2 models
before and after DFT optimisation. . . . . . . . . . . . . . . . . . . . 55
3.12 Electronic densities of state from 320 models of a-SiO2 . . . . . . . . . 56
3.13 Silicon supercell with oxygen atoms inserted to form silica layer. . . . 60
3.14 Si/SiO2 stack after in situ melt and quench of the SiO2 layer. . . . . 61
3.15 Si/SiO2 interface made with a-SiO2 layer placed on top of silicon. . . 62
3.16 Geometrical properties of Si/SiO2 interface made with a-SiO2 layer
placed on top of silicon. . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.1 Crystal structure of α-quartz . . . . . . . . . . . . . . . . . . . . . . . 68

4.2 Atomic and electronic structure of neutral oxygen vacancy in α-quartz 70
4.3 Displacement field of neutral oxygen vacancy in quartz. . . . . . . . . 71
4.4 The dimer configuration of the positive oxygen vacancy in α-quartz. . 73
4.5 Displacement field of the positive oxygen vacancy. . . . . . . . . . . . 74
4.6 The puckered configuration of the positively charged oxygen vacancy
in α-quartz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.7 The negatively charged oxygen vacancy in α-quartz. . . . . . . . . . . 79
4.8 The neutral oxygen vacancy in a-SiO2 . . . . . . . . . . . . . . . . . . 81
4.9 One-electron level of the positive oxygen vacancy in a-SiO2 . . . . . . 83
4.10 One-electron level of the negative oxygen vacancy in a-SiO2 . . . . . . 84
4.11 Thermodynamics of the oxygen vacancy in α-quartz and a-SiO2 . . . . 86
5.1 Ge impurity in α-quartz . . . . . . . . . . . . . . . . . . . . . . . . . 93

5.2 Spin density of the Ge electron trap in α-quartz . . . . . . . . . . . . 94
5.3 Energies of the six different electron trapping configurations in Ge-
doped α-quartz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.4 Interstitial Li atom in α-quartz . . . . . . . . . . . . . . . . . . . . . 97
5.5 [SiO4 /Li]0 centre in α-quartz. . . . . . . . . . . . . . . . . . . . . . . 99
6.1 A schematic of one-dimensional diabatic potential energy surfaces . . 106
6.2 Intrinsic electron trap in α-quartz. . . . . . . . . . . . . . . . . . . . . 109
6.3 Bottom of the a-SiO2 conduction band. . . . . . . . . . . . . . . . . . 111
6.4 Spin density of the intrinsic electron trap in a-SiO2 . . . . . . . . . . . 114
7.1 Experimental optical absorption spectra of electron-irradiated a-SiO2 . 119

7.2 Growth kinetics of the optical absorption bands of electron irradiated
a-SiO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
7.3 Calculated optical absorption spectrum of the GEC in α-quartz. . . . 124
7.4 Calculated optical absorption spectrum of the intrinsic electron trap in
a-SiO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
8.1 A model of a-SiO2 containing an Si–H bond and an Si–O–H bond. . . 133
8.2 Atomic configuration and spin density of the hydroxyl E0 centre and
its precursor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
8.3 Histograms of the bond length distributions around the hydroxyl E0
center. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
8.4 Histogram of the one-electron level of 116 configurations of the hydroxyl
E0 centre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
8.5 Correlations between the one-electron levels and relative stabilities of
the hydroxyl E0 center with the non-bonding Si– –O interaction. . . . 137
8.6 Formation energy of the hydroxyl E0 centre in a-SiO2 . . . . . . . . . . 139
8.7 Histogram of the thermodynamic switching levels of the hydroxyl E0
centre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
8.8 Atomic structure and spin density of the neutral E0 centre. . . . . . . 142
8.9 The neutral E0 centre formed by trapping a hole at an Si–H bond. . . 144
8.10 A graph showing how the energetic position of the Si–H state against
the nearest non-bonded O distance. . . . . . . . . . . . . . . . . . . . 145
8.11 The effect of a nearby proton on the defect levels of the neutral E0 centre.146
List of Tables
3.1 Bond lengths of the Si3 Oy clusters . . . . . . . . . . . . . . . . . . . . 41

3.2 Bond lengths and angles of Sin On clusters optimised using COMB and
ReaxFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3 Bond lengths of Sin clusters optimised using COMB and ReaxFF. . . 44
4.1 Comparison of neutral oxygen vacancy geometries in different models 69

4.2 Comparison of positive oxygen vacancy dimer configuration geometries
in different models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.3 Geometry of the puckered configuration of the positive oxygen vacancy
in different studies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.4 Principal values of the hyperfine tensor for the E0 centre . . . . . . . 78
5.1 Hyperfine splittings of the Li electron trap in α-quartz. . . . . . . . . 99
6.1 Principal values of the hyperfine tensor of the electron trap in a-SiO2 113
8.1 Hydrogen’s reactions with a-SiO2 . . . . . . . . . . . . . . . . . . . . . 148
A.1 BKS Parameterisation for Si and O atoms used in this work. . . . . . 155
List of Publications
The following publications are derived from work presented in this thesis:
1. A. -M. El-Sayed, Y. Wimmer, W. Goes, T. Grasser, A. L. Shluger, “Theoretical

Models of Hydrogen-Induced Defects in Amorphous Silicon Dioxide”, Phys. Rev.
B, 92, 014107 (2015)
2. A. -M. El-Sayed, K. Tanimura, A. L. Shluger “Electron Localization in Amor-

phous SiO2 : Optical Signatures of Deep Traps”, J. Phys.: Condens. Matter 27,
265501 (2015)
3. A. -M. El-Sayed, M. B. Watkins, T. Grasser, V. V. Afanasev, A. L. Shluger,

“Hole Trapping at Hydrogenic Defects in Amorphous Silicon Dioxide”, Micro-
electron. Eng. 147, 141–144 (2015)
4. A. -M. El-Sayed, T. Grasser, V. V. Afanas’ev, A. L. Shluger, “Hydrogen Induced

Rupture of Si–O Bonds in Amorphous Silicon Dioxide”, Phys. Rev. Lett. 114,
115503 (2015)
5. A. -M. El-Sayed, M. B. Watkings, A. L. Shluger, V. V. Afanas’ev, “Nature of

Intrinsic and Extrinsic Electron Traps in Silicon Dioxide”, Phys. Rev. B 89,
125201 (2014)
6. A. -M. El-Sayed, M. B. Watkins, A. L. Shluger, V. V. Afanas’ev, “Identification

of Intrinsic Electron Trapping Sites in Bulk Amorphous Silica from Ab Initio
Calculations”, Microelectron. Eng. 109, 68–71 (2013)
7. S. Ling, A. -M. El-Sayed, M. B. Watkins, V. V. Afanas’ev, A. L. Shluger, “A

Computational Study of Si-H Bonds as Precursors for Neutral E0 centers in
amorphous silica and at the Si/SiO2 interface”, Microelectron. Eng. 109, 310–
313 (2013)
8. A. -M. El-Sayed, A. L. Shluger, “Atomistic Modelling of Defects Implicated in

the Bias Temperature Instability in Bias Temperature Instability for Devices
and Circuits, Chapter 12, Springer (2014)
Acknowledgements
I would like to take this opportunity to offer my heartfelt appreciation to a number of

people whose help and companionship was indispensable during four years over which
this work was accomplished. Firstly, I am hugely indebted to Alex Shluger whose
supervision provided unwavering support and guidance throughout my studies. I am
similarly indebted to Matthew Watkins, whose ready willingness to answer any and all
of my questions saved me from losing my mind. I would also like to thank the rest of
the Shluger group who were always there for a scientific or computational discussion.
For providing me with the opportunity to work collaboratively with them, I thank
Valery Afans’ev and Katsumi Tanimura. I also owe a huge amount of thanks to
Tibor Grasser for our fruitful collaboration and for subsequently welcoming me into
his group as a postdoctoral researcher.
Additionally, I extend my gratitude to the London Centre for Nanotechnology and
the Thomas Young Centre for creating an environment that encouraged and nurtured
my scientific curiosity with their regular seminars and paper discussions.
I thank the UK Doctoral Training Centre in Energy Demand Reduction and the
Built Environment as well as the EPSRC for funding this work. I also thank the Mate-
rials Chemistry Consortium and the London Centre for Nanotechnology for providing
me with the immense computational resources I needed to accomplish this work.
Finally, I thank my family and friends for their staunch support, understanding
and inspiration.
To my parents, family and friends . . .
Introduction
1
1.1 Background & Motivation
This thesis is focussed on the atomistic modelling of defective sites in amorphous sili-
con dioxide (a-SiO2 ) and their ability to trap charge carriers. Although SiO2 is known
as a common constituent of sands around the world, it is widely used in technolog-
ical applications. In particular, its amorphous form is used in electronic devices as
an insulating layer in metal-oxide-semiconductor field effect transistors (MOSFET),
an integral component of integrated circuits (IC). The effective performance of ICs
has been improved over the decades by a process known as scaling, whereby the con-
stituent MOSFETs’ dimensions are scaled down in order to fit more of them within
an IC of a given area. This process has been largely successful, allowing for the light-
Chapter 1. Introduction
ning fast improvements in performance we have seen in electronic devices since the
first realization of the transistor in the late 1950s. Nowadays, the insulating layer
of a-SiO2 in a typical MOSFET in a commercially available electronic device can
have a thickness as small as 1 nm. Due to this extreme scaling, a number of issues
which adversely affect the performance of the device have come to the fore. Col-
lectively, these are termed electronic device reliability issues and include phenomena
1
such as leakage current, f
noise, stress-induced leakage current, and bias tempera-
ture instabilities. Many of these reliability issues implicate defective sites in a-SiO2 .
Due to the very small dimensions of modern electronic devices and the complexity
of their fabrication processes, identifying defects experimentally in these devices is
highly challenging. Computational modelling of materials can offer a comprehensive
insight into the atomistic and electronic structure of defects which can complement
and enlighten experimental results. This thesis uses a number of computational mod-
elling techniques to characterize defects in a-SiO2 which can be involved in electronic
device reliability issues.
In 2009, the UK doctoral training centre in energy demand reduction and the
built environment, from whom I received my funding, was established. The aim of
this centre is to perform research that will result in the reduction of energy demands
in modern society throughout the lifetime of the researcher, fitting in with the En-
gineering and Physical Sciences Research Council’s strategic plan of reducing carbon
emissions and securing energy supplies. The production of electronic devices is an
enormous global industry. For example, according to the consumer electronics asso-
ciation (CEA), the average household in the United States had 24 electronic devices
in 2012. 1 Improving the reliability of the MOSFET, which is an integral component
to all of these 24 devices, reduces the amount of energy wasted while trying to power
the consumer electronics which are so ubiquitous in modern society and helps to
accomplish the goal of energy demand reduction.
In addition, the Modelling of the Reliability and Degradation of Next Generation
Nanoelectronic Devices (MORDRED) project funded by the European Union began
Page 2
in 2011 and was aimed at using modelling techniques to develop our understanding of
the reliability of electronic devices from an atomistic level and integrate these results
into simulations of large-scale electrical circuits. Due to my funding arrangement and
the chance to work with world-class researchers as part of the MORDRED project, I
was drawn into the exciting world of atomistic modelling of defects in order to improve
the reliability of electronic devices.
1.2 Electronic Device Reliability Issues
Figure 1.1: A transmission electron micrograph (TEM) image of a silicon MOSFET. 2
Silicon based transistors are ubiquitous in modern electronic devices, being de-
ployed as voltage-controlled switches. These switches are used to build up highly
integrated logic circuits. Silicon transistors available today are usually of the MOS-
FET variety. They are built upon a crystalline silicon substrate with a thin layer
of amorphous silicon dioxide sandwiched between the Si substrate and a metal. A
transmission electron micrograph (TEM) of a silicon MOSFET is shown in Fig. 1.1.
In addition, the silicon substrate has two terminals, the source and drain, separated
by a channel. The thin oxide layer acts as an insulating layer while the metal acts
as an electrode and is known as the gate. A schematic diagram of a MOSFET is
shown in Fig. 1.2. Depending on the magnitude of the voltage applied between the
gate electrode and the Si substrate, charge can be induced in the channel between
Page 3
Figure 1.2: Typical structure of a MOSFET. 4 The oxide is sandwiched between the
Si substrate and the gate electrode. There are two terminals in the Si substrate: the
source and drain.
the source and drain and the potential difference between them allows a current to
flow. The current flowing between the source and drain is effectively controlled by
the voltage applied across the gate and the Si substrate. 3 Hence the MOSFET can be
used as an amplifier, or more simply as a digital switch. This latter use has resulted
in MOSFETs being used to construct logic circuits which work in conjunction with
each other to form a typical IC.
The dominance of the Si/SiO2 system is due in large part to the extremely high
electrical quality of the interface and the ease with which silica can be grown on
silicon. In this context, high quality means that very few defects exist at the Si/SiO2
interface. Technologies today can produce Si/SiO2 systems with as few as 1010 cm−2
defects at the interface. 5 Other material systems have in the past been considered,
including systems with higher channel carrier mobility, but few systems so far can
match the high electrical quality of the Si/SiO2 interface. In many other systems, the
high concentration of electrically active defects at the interface between two otherwise
very suitable materials renders their use as a transistor ineffective. 6 This was in fact
the case with Ge, a material with a higher carrier mobility and the first material
studied for use as a transistor. In addition, SiO2 has a wide band gap (≈ 9 eV)
providing a barrier which prevents charge carriers from leaving the channel.
Page 4
Figure 1.3: a) Schematic band diagram of tunnelling of electrons through defects in

SiO2 from a metal electrode to a Si substrate. The vertical axis represents the energy
of the electrons. The shaded in region represents the electrons in the system. The
horizontal levels represent defect in the oxide. A single electron highlighted in red is
shown to tunnel from the metal electrode, known as the gate, into the SiO2 defects
and then into the Si substrate. b) Schematic of a defect in SiO2 and how it affects
the charge carriers in the Si substrate. The diagram shows SiO2 deposited on Si. A
negatively charged defect sits in the centre of the SiO2 layer. It has an electrostatic
effect whereby it repels the electrons in the channel of the Si substrate.
Tunnelling charge carriers interact with the oxide layer and the presence of defects
is correlated with the tunnelling of the charge carriers. 7 A simplified band diagram
shown in Fig. 1.3a shows a schematic diagram of charge carriers tunnelling from the
gate to the Si substrate through defects in the oxide. A model, known as percolation
theory, developed by Degrave et al., implicates electrically active defects in the break-
down of the gate oxide. 8 In this model defects are randomly generated in the oxide
when a voltage is applied across it for some period of time. The defects have energy
levels in the oxide band gap and carriers can tunnel through these defects. Breakdown
is achieved when enough of these defects exist and are close enough that they eventu-
ally provide a conduction path for the carriers from the channel to the gate electrode.
This model correctly predicts the dependence of oxide thickness on breakdown but
does not specify which defects are responsible for forming the conduction path.
Defects are also implicated in a number of other MOSFET reliability issues. In
particular, charge trapping at point defects in the oxide have been implicated as the
Page 5
cause of the negative bias temperature instability (NBTI). 9 NBTI is characterised by

a shift in an intrinsic property of a MOSFET known as the threshold voltage after the
application of a negative voltage across a MOSFET at elevated temperatures and for
extended periods of time. It has been known since 1966 10 but has recently come to the
fore as it became more pronounced with the smaller devices that resulted from scaling.
It is almost universally attributed to charge trapping oxide defects and the creation
of interface traps. It is thought that defects in the oxide can trap charge carriers from
the channel. This trapped charge has an electrostatic effect on the channel which
changes the properties of the MOSFET in a detrimental manner. A schematic of how
the defect in an oxide can affect the charge carriers in the semiconductor substrate is
shown in Fig. 1.3b. However, an understanding of its atomistic origins are still not
clear.
1.3 Atomistic Nature of Defects Implicated in Re-
liability Issues
To aid in the understanding of reliability issues, computational methods are used

to construct atomistic models of defects and assess what role they may play in these
issues. Here, we briefly discuss two studies in the literature that attempted to identify
defects which cause reliability issues in electronic devices using atomistic simulations.
Blöchl modelled a number of point defects in a-SiO2 to assess their role in leakage
currents. As mentioned in the previous section, leakage currents are a major concern
for the reliable design of MOSFET devices. 11 Experimental results show that atomic
hydrogen can induce a high concentration of defects in a-SiO2 12–14 suggesting that
H can be somehow involved in these leakage currents. Blöchle modelled hydrogen
related defects using density functional theory (DFT) in a periodic cell of α-quartz cell
containing 24 SiO2 units. They included interstitial hydrogen, oxygen vacancies and
their complexes. To assess their role in leakage currents, he evaluated experimental
data and showed that the defects involved should have a small relaxation energy upon
Page 6
trapping charge and that they should have a defect level where the Si band gap is
located relative to SiO2 . Atomic hydrogen interstitials were found to be unstable in
the neutral charge state, instead being thermodynamically favourable in either the
positive or negative charge states. Its defect levels were found to be too close to
the SiO2 valence band for it to have a significant role in leakage currents. Oxygen
vacancies were also found to have defect levels which render it useless for carrying
leakage currents. However, the hydrogen bridge was found to satisfy both the criteria
for leakage current: it displays small relaxations upon trapping and it has a switching
level in the middle of the Si band gap. His model of the hydrogen bridge is shown in
Fig. 1.4.
Figure 1.4: Atomistic structure of the hydrogen bridge in the neutral, negative and
positive charge states. The relaxations are relatively small compared to relaxations
of the oxygen vacancy defect in SiO2 . 11
Schanovsky et al. investigated the atomistic nature of the defects that could be
involved in NBTI. 15 Using ab initio calculations, they modelled the oxygen vacancy
and hydrogen bridge defects in periodic cells of α-quartz containing 72 atoms. DFT
was used along with the Perdew, Burke and Ernzerhof (PBE) functional and the wave
Page 7
functions were expanded in a plane wave basis set. The total energies calculated for
these defects were then used as parameters to calculate the capture and emission time
constants from a semiconuctor substrate into the defects and compared to experimen-
tal results for NBTI. They showed that neither the oxygen vacancy nor the hydrogen
bridge defects in α-quartz could be responsible for NBTI as their defect levels would
result in much slower capture times. However, it is important to note that the re-
ported results were made in crystalline SiO2 samples, while the measurements were
made on devices containing a-SiO2 .
1.4 Scope of this Work
This thesis presents atomistic calculations of a number of defects in silicon dioxide

which can capture charge. Each defect is extensively characterized and their formation
processes are also discussed. Analysis of these defects shows that they may be involved
in electronic device reliability issues.
We begin by discussing the methods used to perform these calculations. A brief
overview of molecular dynamics is presented. This is followed by an overview of
density functional theory (DFT) and a number of advanced techniques associated with
solving the DFT equations. Classical charge-variable potentials are then discussed.
They are a recent development in the field of computational materials science which
allow the charge of each particle in a system to vary dynamically with the charges
being dependent on the instantaneous geometry of the system. This should allow
one, in principle, to be able to study an element in varying charge states, such as the
transition of silicon from its elemental to fully oxidised form. The embedded cluster
method is then presented, and the methods chapter concludes with a discussion of
climbing-image nudged elastic band (CI-NEB), used to calculate the energetic barriers
of atomistic processes.
The following chapter discusses the process of modelling a-SiO2 and the Si/SiO2
interface. Making amorphous models is challenging as they have no long-range order
Page 8
and therefore its structure has not been atomistically resolved by experiment. In order
to construct a-SiO2 and the Si/SiO2 interface, a charge-variable potential was chosen
based on its stability in a molecular dynamics simulation and its ability to describe
a number of SiOx clusters with Si in varying charge states. Then, a melt-and-quench
method is used to generate 320 models of a-SiO2 . The models’ geometrical and elec-
tronic properties were characterized extensively and compared to experiment where
appropriate. This comparison to experiment is the strongest justification available
at present that these models are indeed representative of a-SiO2 , describing both the
short- and long-range order of the material. The same charge-variable potential was
then used to construct a model of the Si/SiO2 interface. This model was compared
to experiment and shown to agree with a number of experimental observations.
The oxygen vacancy was then calculated in both crystalline and amorphous SiO2
using DFT and the amorphous models made in the preceding chapter. The vacancy is
a well-known defect that has been characterized both experimentally and theoretically
in the literature. The results in this chapter are compared to previous studies and
show excellent agreement, indicating that the methods used can reliably predict defect
structures in SiO2 .
The next chapter discusses electron trapping at impurities in α-quartz. Models
for both Ge- and Li-doped quartz were characterized extensively showing that the
electron trapping at these centres is qualitatively very similar. Electron trapping
at Ge in quartz has been discussed in the literature. 16 However, this is the first
characterization of the atomistic structure of the electron trap in Li-doped quartz.
Electron trapping at intrinsic sites is discussed in the following chapter. Due
to the disorder of a-SiO2 , the bottom of the conduction band is actually partially
localized. It was found that when an extra electron is added to the bottom of the
conduction band, it relaxes into a trapped state deep in the a-SiO2 band gap. The
key to this relaxation was found to be a wide O–Si–O angle, which naturally exists in
a-SiO2 . Electron trapping at intrinsic sites was also studied in α-quartz. However, due
to all the O–Si–O angles being the same, the trapping process requires overcoming
Page 9
an energetic barrier. The character of localization at intrinsic sites was shown to

be similar to that at the extrinsic sites. The concentration of electron traps was
calculated in a model of a-SiO2 containing almost 500,000 atoms and was shown to
agree well with experimental estimates of electron traps involved in electronic device
degradation processes.
The optical absorption spectrum of the intrinsic electron trap was calculated using
the embedded cluster method and time-dependent density functional theory (TD-
DFT) for the first time. Characterization of the electronic structure of the electron
trap in the embedded cluster model agree very well with the previous chapter’s results
calculated in a periodic model. The calculated optical absorption spectrum shows a
strong excitation at 3.7 eV from the electron trapping state into the bottom of the
a-SiO2 conduction band. In addition, there is also a strong excitation from the top of
the a-SiO2 valence band into the electron trap. These excitations correlate strongly
with experimental optical absorption spectra measured at low temperature.
Interactions of hydrogen with a-SiO2 are presented in the final chapter. Both
atomic and molecular hydrogen were shown to interact with a-SiO2 to produce similar
defect states. For the first time, calculations show that atomic hydrogen interacts
with strained Si–O bonds, breaking them and resulting in a 3-coordinated Si facing a
hydroxyl group. A model of the formation and annihilation of this defect is presented
in this chapter. Molecular hydrogen was shown to create Si–H bonds when introduced
into an SiO2 melt. The states of these bonds sit near the top of the valence band
and were shown to trap a hole. This process releases a proton and leaves behind
a 3-coordinated Si whose state in the band gap is highly affected by how far it is
from the nascent proton. The role of the hydrogen related defects in electronic device
reliability issues is discussed.
Page 10
Theoretical Background and Methods
2
2.1 Introduction
The aim of this chapter is to give a detailed overview of the theoretical methods that
were used to obtain the results presented in this thesis. Modelling point defects in
crystalline and amorphous materials is a complicated affair due to their deviation
from an ideal model, ths requiring the use of advanced modelling techniques. The
chapter begins by covering the theory of molecular dynamics, which is used to evolve
the positions of particles in a many-body system. It then goes on to explain the basics
of density functional theory (DFT) and further approximations which were employed
in the calculations presented in this thesis. This is followed by an introduction to
classical potentials and the recently developed charge-variable potentials COMB and
Chapter 2. Theoretical Background and Methods
ReaxFF. The embedded cluster method is briefly described and the chapter concludes
with the nudged elastic band method (NEB) for calculating reaction pathways.
2.2 Molecular Dynamics
Molecular dynamics is a computer simulation method in which the force on each

particle in a system is calculated and integrated to allow each particle to evolve its
position and velocity, generating a trajectory for the system. Using some form for
the potential energy, which can be described classically (see section 2.4) or quantum
mechanically (see section 2.3), the force on each particle is calculated as the differential
of the vector sum of its interactions with all other particles in the system:
N
X dU
F (r) = − , (2.1)
i=1
dr i
where F (r) is the total force on the atom, U is the potential energy, N is the total
number of particles, and r i are the coordinates of atom i. The acceleration can be
determined from the force according to Newton’s second law of motion, generating
the trajectory of motion of the system:
d2 r
F = ma = m . (2.2)
dt2
Using this method with a potential energy surface which depends on each particle’s
position relative to all other particles constitutes a complex many-body problem.
The coupling of the particles’ motions makes it extremely difficult to calculate the
trajectory analytically, so one must use a numerical method to calculate the trajectory.
One of the most popular and widely used methods of integrating the equations of
motion when using pairwise continuous potential energy surfaces is the Verlet algo-
rithm. 17 The integration is broken down into small steps separated by a fixed time
which we call the time step. Writing out the position at the next and previous time
step as Taylor expansions of position up to third order, we obtain for the next time
Page 12
step:
1 1
r(t + δt) = r(t) + δtv(t) + δt2 a(t) + δt3 b(t) + . . . , (2.3)
2 6
where r(t) are the positions; δt is the time step; v(t), a(t), and b(t) are the first (also
known as velocity), second (known as acceleration), and third order differentials of
position, respectively. The position at the previous time step is:
1 1
r(t − δt) = r(t) − δtv(t) + δt2 a(t) − δt3 b(t) + . . . . (2.4)
2 6
Adding equations (2.3) and (2.4) gives:
r(t + δt) = 2r(t) − r(t − δt) + δt2 a(t). (2.5)
Using this equation, the new positions can be calculated from the acceleration and
each particle’s previous positions. It can be seen in equation 2.5 that the 3rd order
terms in equations 2.3 and 2.4 cancel out. This in fact makes the use of equation
2.5 an order more accurate than the use of equation 2.3 alone. The velocities are not
explicitly defined in the Verlet integration scheme but can be calculated by dividing
the difference between the positions at t+δt and t by 2δt:
r(t + δt) − r(t)

v(t) = . (2.6)
2δt
Using the information laid out in this section one can set up an algorithm to evolve
each particle’s position and obtain a trajectory from the particles’ forces. To obtain
the forces one can use a variety of methods; for example, a classical (see section 2.4)
or a first principles quantum mechanics approach (see section 2.3).
Page 13
2.3 Density Functional Theory
DFT is a theory providing information on the electronic properties and chemical

bonding of materials. It is, in principle, an ab initio method, meaning that it requires
no parameters except for the fundamental properties of the material being studied,
e.g. atomic mass, mass of an electron, etc. DFT follows from the Hohenberg-Kohn
theorems, 18 the first of which states that there is a one-to-one mapping of the density,
ρ(r), of a system to its external potential, vext ; in other words, the density of a system
uniquely determines its external potential. This means that all the properties of a
system can be determined from its density, and the total energy of a system can be
written as a functional of its density as:
E[ρ(r)] = E T [ρ(r)] + E V [ρ(r)] + E H [ρ(r)] + EXC [ρ(r)], (2.7)
where E T [ρ(r)] is the kinetic energy, E V [ρ(r)] is the potential energy, E H [ρ(r)] is
the Hartree energy, and E XC [ρ(r)] is the exchange-correlation energy of the system.
Writing the energy of a system as a functional of its density reduces a many-body
problem into an equation which depends on the three spatial variables of the density.
The final term, E XC [ρ(r)], is the exchange-correlation energy functional and contains
all the complicated many-body effects that have not been described within all the
other energy terms. The second Hohenberg-Kohn theorem states that the true ground-
state density of the system is the density which minimises the energy of the system,
i.e.:
E[ρ(r)] ≥ E0 for ρ(r) 6= ρgs (r),
(2.8)
E[ρ(r)] = E0 for ρ(r) = ρgs (r),
where ρgs (r) is the ground-state density and E0 is the ground state energy. Thus, by
minimising equation 2.7, the ground-state density and energy can be obtained.
Equation 2.7 means that the total energy of the system can be calculated from
the 3 dimensional density. However, due to a lack of approximations for E T [ρ(r)],
these calculations are, as yet, impractical for the accurate description of materials. In
Page 14
a landmark paper, Kohn and Sham introduced an auxiliary system of non-interacting

electrons which have the same ground-state density as the interacting system which
we can call the Kohn-Sham system (KS). 19 As the KS particles are non-interacting,
the wave function can be written as a Slater determinant:
ψ1 (r1 ) ψ2 (r1 ) . . . ψN (r1 )
1 ψ1 (r2 ) ψ2 (r2 ) . . . ψN (r2 )

Ψ(r1 , . . . , rN ) = √ .. .. .. .. , (2.9)
N! . . . .
ψ1 (rN ) ψ2 (rN ) . . . ψN (rN )
where ψi (rj ) are single particle orbitals. The electron density can then be obtained
from these orbitals using the Born interpretation of the wave function:
N
X
ρ(r) = |ψi (r)|2 = |Ψ|2 . (2.10)
i=1
The total energy of the system can be evaluated as the expectation value of the
Hamiltonian of this wave function. An external potential, vef f is chosen so that the
true ground-state density of the system is reproduced from the KS system. The total
energy of the system of non-interacting electrons is thus:
Egs = hΨgs [ρ]|Ĥ|Ψgs [ρ]i, (2.11)
where:
h̄∇2
Ĥ = − + vef f , (2.12)
2m
and:
Z
ρ(r) δExc [ρ]
vef f (r) = vext (r) + vHartree (r) + vxc (r) = vext (r) + dr 0 0
+ , (2.13)
|r − r | δρ(r)
where vext is the interaction of the electrons with the ions in the system, vHartree is
the electron-electron repulsion and vxc is the exchange-correlation potential. Using
the Kohn-Sham formulation of the Schrödinger equation one can calculate the single-
Page 15
particle orbitals: 19
h̄∇2

− + veff (r) ψi (r) = εi ψi (r), (2.14)
2m
where veff (r) is the effective potential, ψi (r i ) is the Kohn-Sham orbital and εi is the
eigenvalue of the Kohn-Sham orbital. Equation 2.14 is calculated self consistently
starting with an initial guess for the electron density which can be inserted into the
effective potential. Once the new density converges with the density from the previous
step, the Kohn-Sham equations have converged and the ground-state density has been
obtained.
Although DFT has had considerable success describing ground-state properties, it
is largely unsuitable for the calculation of excited state properties. However, the time-
dependent extension to DFT (TDDFT) has been used to calculate excited state prop-
erties rather successfully. 20–22 The theory behind TDDFT is, unsurprisingly, rather
similar to DFT. A one-to-one mapping between the time-dependent ground state and
a time-dependent external potential up to a time-dependent but spatially indepen-
dent function is established in the time-dependent extension to the Hohenberg-Kohn
theorem: the Runge-Gross theorem. 23 Total energy is not a conserved quantity in
TDDFT, thus there is no analogue of the variational principle, and instead the ener-
gies are defined as the stationary points on the quantum mechanical action integral:
Z t1
δ
A[Ψ] = dthΨ(t)|i − Ĥ(t)|Ψ(t)i, (2.15)
t0 δt
where Ψ(t) is the time-dependent wave function (which can be uniquely obtained from
the time-dependent density) and Ĥ(t) is the Hamiltonian of the system. The time-
dependent density of the system is then obtained by finding the stationary points on
the action potential, i.e. the points where A[Ψ] = 0. To obtain the time-dependent
density, a Kohn-Sham system of non-interacting particles in a time-dependent poten-
tial is constructed, analogous to the Kohn-Sham system in DFT:

δ 1 2
i ψi (t) = − ∇ + vef f (t) ψi (t), (2.16)
δt 2
Page 16
where ψi (t) are the time-dependent single particle orbitals and vef f (t) is the time-
dependent effective potential:
vef f = vext (t) + vHartree (t) + vxc (t), (2.17)
where vext (t) now includes a time-dependent external potential perturbing the system
(e.g. a laser) as well as the static potential generated by the ions, vHartree (t) is
the time-dependent Hartree potential, and vxc (t) is the time-dependent exchange-
correlation potential. The time-dependent density can then be constructed from the
time-dependent single particle orbitals, once again analogous to DFT:
N
X
ρ(t) = |ψi (t)|2 , (2.18)
i=1
where ρ(t) is the time-dependent density. By self-consistently calculating a stationary

point on the action integral, one obtains the time-dependent density of the system.
To calculate the excitation energies and their respective oscillation strengths, the
linear response of the time-dependent density has to be examined. After a number
of algebraic manipulations and application of the Tamm-Dancoff approximation, this
leads to the eigenvalue equation: 21
AX = ωX, (2.19)
where ω are the excitation energies, and the eigenvectors X can be used to calculate
the oscillator strengths. The elements of matrix A are given as:
Aia,jb = δij (a − b ) + (ij|jb) + (ia|fxc |jb), (2.20)
where (ab|bc) are two electron integrals represented in Mulliken notation, and fxc is
the time-dependent exchange-correlation kernel. Solving this eigenvalue problem is
equivalent to finding the excitation energies which are at the poles of the density
Page 17
response equations. 24
2.3.1 Gaussian Plane Waves Method
In a typical DFT calculation, the Kohn-Sham orbitals are described within a chosen
basis set of functions. Two popular functions typically used for solving the Kohn-Sham
equations are Gaussian and plane wave functions. Each type has its advantages and
disadvantages. Gaussian basis sets are atomically centred functions. They allow for
a compact description of the wave function and efficient algorithms exist to analyti-
cally calculate matrix integrals. However, Gaussian basis sets do not describe orbitals
properly near the nucleus and can suffer due to the incomplete nature of the basis
set. Plane waves are used as basis sets due to their indifference to the atomic coor-
dinates and species. Their periodic nature means they are a natural choice of basis
set when describing periodic systems, and plane waves, in principle, form a complete
basis. Fast Fourier transform techniques allow for a faster calculation of the Hartree
energy and checking the convergence of a calculation using a plane wave basis set is
a straightforward matter. However, describing the density properly in a plane wave
basis usually requires a larger number of functions relative to a Gaussian basis set
and plane waves represent regions of empty space with the same accuracy as a region
which contains an atom, leading to unnecessary computational expense. Although in
principle plane waves form a complete basis, it is necessary to truncate the basis set
in order to make a calculation feasible.
The Gaussian and plane waves method 25 (GPW) is designed so that the density
of the system is represented in both bases in order to use the more time-advantageous
basis set for the various elements of the calculation while maintaining accuracy. The
following describes both Gaussian and plane wave basis sets, followed by their com-
bination to make the GPW method.
Page 18
Gaussian Basis Sets
The density represented in a Gaussian basis set is written as:
N
X
ρ(r) = |ψi (r)|2 , (2.21)
i=1
where ρ(r) is the density, N is the number of electrons, and ψi (r) is the i’th molecular
orbital. The molecular orbitals can be approximated using the linear combination of
atomic orbitals (LCAO) method, where each molecular orbital is described as:
X
ψ= Cmi ψj , (2.22)
j
where Cmi is the mixing coefficient of the atomic orbital and ψi is the atomic orbital
which is a contracted Gaussian, made of:
X
ψ= Cµn gµ , (2.23)
µ
where Cµi is the mixing coefficient and gµ is a primitive Gaussian function (gµ =
2
. . . e−αr ). In a typical calculation, the coefficients of the molecular orbitals, Cmi , are
varied until the ground-state density is obtained while the mixing coefficients of the
contracted Gaussian, Cµn remain fixed.
Plane Wave Basis Sets
The density represented in a plane wave basis set is:
1X
ρ̃(r) = ρ̃(G)eiG·r , (2.24)
Ω G
where ρ̃(r) is the density, Ω is the volume of the cell, G are the reciprocal lattice
vectors and ρ̃(G) are the expansion coefficients. The density, ρ̃(r) is determined from
ρ(r) using Fourier transforms on a real-space grid.
Page 19
Combining Gaussians and Plane Waves to Make the Gaussian Plane Waves
Method
The electrons in the system can be split into those near the nucleus whose form does
not change much during a calculation, known as core electrons, and those which sit
away from the nucleus which are involved in chemical reactions, known as the valence
electrons. Replacing the core electrons with what is known as a pseudopotential 26
provides a smoothly varying density which can be easily mapped from the Gaussian
basis set to a plane wave basis set. Fewer Gaussian functions are required to produce
the characteristic cusp behaviour at the nucleus, and due to the smoothly varying na-
ture of pseudopotentials, less plane waves are required. This allows for easier mapping
of the density from a Gaussian to a plane wave basis set. The Kohn-Sham equation
(see equation 2.7) within these two basis sets is calculated as:
X 1
E[ρ] = P µv hψµ (r)| − ∇2 |ψv (r)i + hψµ (r)|Vloc
PP
|ψv (r)i + hψµ (r)|VnlP P |ψv (r)i+
µv
2
X ρ̃ ∗ (G)ρ(G) Z
2πΩ + εxc (r)dr,
G
G
(2.25)
where P µv hψµ (r)|− 12 ∇2 |ψv (r)i is the kinetic energy, hψµ (r)|Vloc
PP
|ψv (r)i+hψµ (r)|VnlP P |ψv (r)i
are the local and non-local parts of the pseudopotential, 2πΩ ρ̃∗(G)ρ(G)
P
G
is the Hartree
R G
energy and εxc (r)dr is the exchange correlation energy. The density within a Gaus-
sian basis set is used to calculate the kinetic and potential energy analytically. Use of
the Gaussian product theorem makes this a relatively trivial endeavour. Mapping the
density onto the plane wave basis allows one to use fast Fourier techniques to calculate
the Hartree energy. This method significantly reduces the computational cost of DFT
calculations allowing one to study much larger systems than if one were to use either
basis set on its own. The GPW method is implemented in the Quickstep 27 module
within the CP2K code. 28
Page 20
2.3.2 Hybrid Functionals and the Auxiliary Density Matrix
Method
A huge number of approximations exist in the literature and a judicious choice of func-
tional is necessary for the accurate calculation of a material’s properties. Exchange-
correlation functionals were initially based on local measurements of the homogeneous
electron gas 29 and gradient-corrected modifications. 30 These functionals, however, suf-
fer from two prominent shortcomings which are particularly relevant to the study of
defects in solid state systems. They tend to give wrong descriptions of the band gap,
with both LDA and GGA underestimating it by up to 40%. 31 In addition, these func-
tionals have trouble describing the localized defect states which are known to exist
experimentally. 32,33 Recent developments have resulted in more advanced approxima-
tions, such as the set of functionals known as hybrid functionals which use a portion
of Hartree-Fock (HF) exchange which is explicitly calculated for all electrons in the
system. The band gaps and localized defect states calculated using these hybrid func-
tionals show significant improve compared to the local and gradient-corrected approxi-
mations. The main results of this thesis have made use of hybrid functionals, including
the B3LYP and PBE0 TC LRC hybrid functionals 34 which both contain a portion of
HF exchange energy. It is calculated from the density matrix, P =< ψµ |ρ|ψλ > and
the two electron integrals as:
1 X µσ λv
ExHF [P ] = − P P hµv|λσi, (2.26)
2 λσµv
where P is the density matrix and hµv|λσi are the exchange integrals, described as:
Z Z
1
hµv|λσi = ψµ∗ (r 1 )ψv∗ (r 1 ) ψλ (r 2 )ψσ (r 2 ), (2.27)
r 12
where ψ are the single-particle orbitals. The HF exchange energy is computationally

expensive, scaling to the fourth power of the number of basis functions used in a
calculation as can be seen above. The form of the exchange-correlation energy in the
Page 21
PBE0 TC LRC functional is:
P BE0 T C LRC
Exc = αExHF,T C + αExP BE,LRC + (1 − α)ExP BE + EcP BE , (2.28)
where α is a parameter that controls the amount of HF exchange, typically defined

as 0.2 in the literature and is the value used in this thesis unless otherwise stated. 34
1
T C is a truncated exchange operator where the r12
term becomes zero if r12 exceeds a
pre-defined rc cut off value. The exchange-correlation energy in the B3LYP functional
is: 35
B3LY P
Exc = ExLDA + a0 (ExHF − ExLDA ) + ax (ExGGA − ExLDA ) + EcLDA + ac (EcGGA − EcLDA ),
(2.29)
where a0 = 0.20, ac = 0.72, and ac = 0.81. LDA represents exchange-correlation in
the local density approximation, 29 while GGA represents gradient-corrected exchange-
correlation. 30
To reduce the computational cost of calculating HF exchange integrals, the aux-
iliary density matrix method (ADMM) was employed in many of the calculations
presented in this thesis. It constructs a much smaller density matrix using an aux-
iliary basis set. The calculation of HF exchange can be performed using this basis
set and the smaller density matrix, thus reducing the number of calculations needed.
The auxiliary basis set is described as:
X
ψ̃ = C̃µi ψ̃i (r), (2.30)
µ
where ψ̃ is the wave function in the smaller auxiliary basis set, C̃µi is the Gaussian
orbital coefficient and g̃i (r) is the Gaussian orbital. The density matrix is constructed
from the orbital coefficients as:
X
P̃ µv = C̃µi C̃vi . (2.31)
i
Page 22
The orbital coefficients for the auxiliary basis set and auxiliary density matrix can
be obtained by minimising the square difference between the wave functions in the
original basis set and the auxiliary basis set. The HF exchange energy can then be
calculated in this basis set as:
ExHF [P ] = ExHF [P̃ ] + (ExHF [P ] − ExHF [P̃ ]),

(2.32)
≈ ExHF [P̃ ] + (ExDF T [P ] − ExDF T [P̃ ]),
where ExHF [P̃ ] is the HF exchange energy while ExDF T [P ] is the exchange energy
calculated using GGA. It is assumed that the difference in the exchange calculated
using different basis sets using GGA and HF is approximately the same. Using this
auxiliary density matrix and the bottom half of equation 2.32 can greatly reduce the
cost of a DFT calculation using a hybrid functional.
2.4 Classical potentials
The problem with ab initio calculations, like DFT discussed in section 2.3, is that
they are computationally expensive. The advantage of using a classical potential is
that one can model systems containing up to billions of atoms and simulations can
be run for up to a microsecond. 36,37 To calculate equilibrium energies and structures
one can use a classical interatomic potential along with some simulation method, e.g.
Molecular dynamics or minimisation algorithms. The disadvantages of using classical
potentials is that you don’t calculate the electronic structure of the system, which
is essential when modelling defects as their electronic structure can deviate strongly
from that of the ideal material. Interatomic potentials have to be parametrised,
usually to a particular quantity of interest. For instance, the BKS potential, 38 used to
model silica, can reproduce the structure and density of crystalline silica polymorphs
quite well but has trouble describing its phononic properties and certain crystalline
transitions. 39 Potentials are designed to be as transferable as possible but this is
Page 23
difficult to do as it is hard to capture the physics of a species through all possible

chemical environments using a simple potential. Parameterizing a potential is one
of the less glamorous aspects of modelling but is needed to obtain results that can
reliably explain physical phenomena.
Often times it is more convenient to start off with a well parameterized potential
and fine tune it to accurately describe the desired quantity. An interatomic poten-
tial proposed by Watanabe 40 was developed as an extension to the Stilinger-Weber
potential, which accurately models crystalline silicon, to model systems containing
both silicon and silica and has been used to study the structure of Si/SiO2 interface.
Many potentials exist which describe bulk silica and silicon but few exist which are
transferable enough to describe the transition of silicon from its elemental to its fully
oxidised silica form.
2.4.1 Charge variable potentials
A problem which can prevent a potential from being truly transferable is the varying
charge state of a species when it is in different materials, e.g. Si in bulk silicon and
SiO2 . Conventional potentials contain the charge of a species as a fixed parameter,
thus a particle can not be simulated undergoing a change in charge state. To overcome
this problem, Rappé and Goddard developed a scheme which allows a particle’s charge
to be calculated from the instantaneous geometries of the particles in a system. 41 The
scheme is based on the electronegativity equalisation principle which states that the
electronegativities of species coming together becomes constant.
χA = χB = . . . = χN . (2.33)
The effective electronegativity of a system of particles is the geometric mean of the

electronegativities of the species. 42 One can define an energy for each atom with
Page 24
respect to its charge as a Taylor expansion about the neutral atom:

δE 1 2 δE
Ei (q) = Ei0 + qi + qi0 . (2.34)
δq i0 2 δq i0
The energy required to make an ion of charge +1 is assigned to the ionisation potential
while a charge of -1 is termed the electron affinity. They can be calculated as:
Ei (0) = Ei0 , (2.35)

δE 1 δE
Ei (+1) = Ei0 + + = Ei0 + IP, (2.36)
δq i0 2 δq i0

δE 1 δE
Ei (−1) = Ei0 − + = Ei0 + EA. (2.37)
δq i0 2 δq i0
Mulliken defined electronegativity as the average of the electron affinity and ionisation
potential. Using this definition and the results above we can say:
δ2E

1 δE δE δE
= − EA = IP − →2 = IP + EA →,
2 δq 2 δq i0 δq i0 δq i0

δE 1
= (IP + EA) = χoi , (2.38)
δq i0 2
and:
1 δ2E 1 δ2E

δE
= IP − = EA + →,
δq i0 2 δq 2 2 δq 2
1 δ2E

2
= IP − EA = Jiio , (2.39)
2 δq
where χ0i and Jii0 are the Mulliken electronegativity and the self-Coulomb interaction
(also sometimes referred to as the idempotential). Both these quantities can be ob-
tained from experimental data. To determine the equilibrium charge of the species in
a system, we now write the total electrostatic energy of a system as:
N N
X 1 X
E(q1 , ..., qN ) = Ei0 + χoi qi + Jiio qi2 + qi qj Jij . (2.40)
i
2 i<j
Page 25
We can now minimize this energy with respect to charge and invoking the electroneg-
ativity equalisation principle and a charge totality condition (usually set to zero). We
can solve for N simultaneous equations:
N
δEi X
χi = = χoi + Jiio qi + qj Jij , (2.41)
δqi i<j
N
X
C= qi . (2.42)
i
We now have a system whereby the charge of a particle in a system can be calculated
from its geometry using only three parameters: electronegativity, self-Coulomb term
and the covalent radius. This method can be applied to an interatomic potential
so that charge can be calculated at each step of a simulation, allowing each parti-
cle’s charge to evolve dynamically. This method has been applied to a wide range
of materials, including lipid molecules, 43 proteins, 44 solid-state materials 45–47 and hy-
drocarbons. 48 Using this method, one can parameterize an atomic species so that it
can be studied across a range of charge states.
2.4.2 Charge-optimized many-body potential (COMB)
The COMB potential 45 was developed to be able to study Si/SiO2 systems on a large
length and time scale. Many interesting processes occur on a time- and/or length scale
that would be far too computationally expensive to calculate using ab initio methods
and, as discussed above, typical interatomic potentials can not describe an atomic
species in its various charge states. The COMB potential allows the charge of the
species to vary in order to overcome this problem. COMB is based on the Yasukawa
modification to the Tersoff extended potential for Si/SiO2 . The problem with the
Tersoff potential is that it assigns static charges, so Yasukawa coupled it with Rappé
and Goddard’s charge equilibration method to allow each particle’s charge to vary
over the course of a simulation. 41 Yasukawa’s potential 49 was used to study silicon,
silica and Si/SiO2 systems and these studies revealed challenges which the COMB
potential aims to overcome:
Page 26
• The lowest energy state of bulk silicon in the diamond lattice using the Yasukawa
potential results in charged Si atoms, i.e. the lowest energy state is not charge
neutral. The Si atoms are charged ±3.3 depending on the sublattice they sit
in. The charge neutral Si system is actually an energy maximum using this
potential.
• The relative stabilities of silica polymorphs do not match experimental data.

In particular, α-quartz is not predicted to be the lowest energy polymorph, in
contrast to what is known experimentally. 50
• The Yasukawa potential overestimates the O–Si–O units in the polymorphs α-

and β-cristobalite.
Functional Form
The interatomic potential energy in the Yasukawa and COMB potential is calculated
as:
Vij (r ij , qi , qj ) = UijR (r ij ) + UijA (r ij , qi , qj ) + UijI (r ij , qi , qj ) + UijV (r ij ), (2.43)
UijR (Rij ) = fsij Aij e(−λij rik ) , (2.44)
UijA (r ij , qi , qj ) = −fsij bij Bij e(−αij rij ) , (2.45)

fLij ηi ηj qi qj
UijI (r ij , qi , qj ) = , (2.46)
4πε0 r ij
21
fL CV DWi CV DWj
UijV DW (r ij ) = ij , (2.47)
r 6ij
where UijR is the repulsive energy, UijA is the attractive energy, UijI is the ionic bond
energy and UijV DW is the van der Waals energy. The charges used to calculate the
attractive and ionic energy are obtained using the charge equilibration scheme. Three-
body effects are incorporated into bond-order, bij , which is used as a weighting factor
Page 27
when calculating UijA :
" !ni #−1/2ni

X
bij = 1 + βi ζijk g(θjik ) . (2.48)
k6=i,j
The parameters used in the COMB potential are derived from experimental data
and ab initio calculations of α-quartz. The lowest energy silicon structure is charge
neutral using these parameters and the relative stabilities of the silica polymorphs
match experiment. Lattice constants of α-quartz match experiment but can still
be improved for other SiO2 polymorphs. The COMB potential was designed to be
extremely transferable with the ability to study Si/SiO2 interfaces and has also been
parametrised for other elements. 46
2.4.3 ReaxFF
ReaxFF is another potential that has been developed to study silicon in varying oxida-
tion environments. 51,52 It is an empirical potential whose parameters were determined
from B3LYP calculations with a 6-31G** basis set on SiOx clusters with silicon in
varying oxidation states. DFT calculations on crystalline Si and SiO2 polymorphs
were also included in the parameterisation training set. The parameters were opti-
mized using a successive one-parameter search technique. 53 Parameterising to clusters
which contain Si in varying oxidation states is meant to include structures which could
exist at the interface of an Si/SiO2 system and which is not represented by either bulk
crystalline silicon or silica. ReaxFF describes covalent interactions in terms of a dy-
namically changing bond order with the charges of each particle determined using
the Rappé-Goddard method, allowing the use of classical MD to study the Si/SiO2
interface.
No connectivities are defined with this potential but instead a bond-order term is
calculated which depends on the atoms’ geometries. The energy terms depend on the
Page 28
Figure 2.1: Silicon-silicon bond-order dependence on interatomic distance. 51
bond order which is calculated as:
0 0 0 0
BOij = BOijσ + BOijπ + BOijππ ,
pbo,2 pbo,4 p
r ij r ij r ij bo,6
= exp pbo,1 + exp pbo,3 + exp pbo,5 , (2.49)
r σ0 r π0 r ππ
0
0 0σ
where BOij is the total bond order, BOij is the bond order contribution from the
0π 0ππ
σ bonding, BOij is the bond order contribution from π bonding, BOij is the bond
order contribution from ππ bonding and pbo,x are fitting parameters. Fig. 2.1 shows
how the total bond order is affected by the contributions from different bondings, and
how the bond order for different bondings varies with interatomic distance.
The energy of the system is split into the following components:
ESystem = EBond +EOver +EU nder +ELP +EV al +EP en +ET ors +EConj +EV DW +ECoulomb .
(2.50)
The interactions ascribed to covalency (e.g. bond-stretching, angles, torsion) contain
contributions from σ, π and π −π bonds. For example, the bond stretching interaction
Page 29
is described as:
σ
Ebond = −Deσ BOij π
exp [pbe,1 (1 − (BOij )pbe,2 )] − Deπ BOij ππ
− Deππ BOij , (2.51)
where De and pbe are fitting parameters. Each bond order has a different dependence
on interatomic distance. Describing covalent bond interactions in terms of a bond
order (dependent on instantaneous geometry) allows for the dynamic formation and
annihilation of bonds. The non-bonding interactions (Coulomb and van der Waals
interactions) do not use the bond-order but use the charge calculated using the method
developed by Rappé and Goddard.
2.5 Embedded Cluster Method
Both classical potentials and ab initio methods were used to model materials through-
out this thesis. As mentioned in section 2.4, classical potentials can be used to model
a system containing up to a billion atoms. Well parametrised classical potentials are
sufficient to provide a good description of the atomistic structure and thermodynamic
properties of many materials. On the contrary, DFT can typically model systems con-
taining up to a few thousand atoms. However, it can provide information regarding
the electronic structure of the system which is inaccessible with classical potentials.
In addition, a charged point defect in a crystalline material is much more difficult to
describe using classical potentials due to the distortion they may introduce relative
to the ideal lattice to which the potential is usually parameterized.
To study the long range interactions of a point defect in an ideal lattice, it would be
useful to be able to use both classical and ab initio theories in the same calculation. A
system can be divided into two regions, with one described by a quantum mechanical
(QM) and the other by a classical description. Methods of this type are generally
referred to as quantum mechanical molecular mechanics (QM/MM) methods. The
embedded cluster method is a QM/MM method which is particularly suitable for point
defects in solid state systems due to approximations which provide correct Madelung
Page 30
potential variations in the QM region.

The starting point of the embedded cluster method is the creation of a finite
nano-sized cluster which describes the material. This nanocluster is divided into two
regions - Regions I and II. Region I is located at the centre of the nanocluster and
is the region of interest. It may, for example, contain a point defect and its local
surroundings. Region I is further subdivided into three regions:
1. A QM region which contains the point of interest and as much of its surroundings
as computationally permissible.
2. A region in which all the atoms are described using classical interatomic poten-
tials.
3. An interface region which is sandwiched between the QM and classical region

and interacts with them both.
Region II surrounds Region I and is made up of atoms whose position remains fixed
and which are described using a classical interatomic potential. Fig. 2.2 shows an
example of Region I in the embedded cluster model and highlights the different regions
within it. All atoms in Region I are allowed to move in an embedded cluster calculation
while atoms in Region II (which are not shown in Fig. 2.2) are static. The function
of the atoms in Region II is to ensure that all sites in the region of interest in the
QM cluster experience correct Madelung potential variations. This essentially allows
for Region I to experience forces as though it were embedded in an infinite crystal,
despite the entire system being finite.
To allow the QM region to interact with the classically polarizable region, an
interface region is defined between them as a set of pseudoatoms which bind both
regions. In the case of a mixed ionic-covalent material (such as SiO2 : the main focus
of this thesis), the pseudoatom is described by a single electron and a short-range
repulsive potential which are both parameterized to provide a good approximation for
the interatomic interactions in an ideal lattice. The pseudoatom also interacts with
the classically polarizable region. This interaction is described as a sum of Coulomb
Page 31
Figure 2.2: A cross-section through an example of Region I in an embedded cluster

model containing both QM and classical atoms. The quantum cluster atoms are
highlighted in blue, the pseudoatom interface atoms are highlighted in white while
the classical atoms are highlighted in red. The quantum cluster is embedded into the
classical region and all atoms in Region I are allowed to relax.
and short-range classical Morse interactions. The interactions between atoms in the
classical region and atoms in the QM region are included as Coulomb terms in the
Hamiltonian of the QM region. The total energy of the system is given as:
ET otal = EQM + ECl , (2.52)
where EQM and ECl are the energies of the QM and classical region respectively. EQM
is given as:
Coul Short
EQM = hψ|Ĥ0 + Venv |ψi + Venv , (2.53)
Short
where Ĥ0 is the Hamiltonian of the QM cluster only, Venv is the short range inter-
Coul
action between the classical and QM atoms, and Venv is the external potential which
Page 32
includes interactions between the classical cores and shells with the QM atoms:
 
N n NQM
Coul
Xcore X qicore X qicore ZjQM
Venv =  + +
i j
|ricore − rj | j
|ricore − rjQM |
 
NX n NQM
shell X qishell X qishell ZjQM
 + ,
i j
|rishell − rj | j
|rishell − rjQM |
where n and NQM are the numbers of electrons and QM atoms in the QM region; Ncore
and Nshell are the numbers of cores and shells in the classically polarizable region; Z,
q core and q shell are the charges of the QM nuclei, the classically polarizable cores and
the classically polarizable shells, respectively. This total energy can be minimised to
obtain the ground state structure of the system in the embedded cluster method.
2.6 Nudged Elastic Band Method
The activation energy of a reaction is an extremely important parameter in the context

of chemical dynamics, allowing accurate predictions of the rate of a reaction. The
concept of an activation energy stems from the Arrhenius equation and is essentially
the energy required in order for a reaction to proceed. The Arrhenius equation can
be written as:
EA
−K
kab = Ae BT , (2.54)
where kab is the rate of the reaction, EA is the activation energy, and KB T is the
Boltzmann constant scaled by temperature. Further research led to the development
of harmonic transition state theory, in which the activation energy also plays a central
role. 54 For these reasons, there has been a large research effort into methods that can
provide minimum energy pathways (MEP) of a reaction from which the activation
energy can be extracted. A MEP is the trajectory between reactants and products
that follows the lowest energy along the potential energy surface between them. The
activation energy(ies) can be determined from the MEP as the first order saddle points
(there may be more than one activation energy) along the MEP.
Page 33
A method that has been used to successfully calculate MEPs for many reactions
in various materials is the climbing image nudged elastic band method (CI-NEB).
Although a number of methods are available in the literature for the calculation
of MEPs, such as the linear synchrous transit method 55 and eigenvector following
methods, 56 the use of CI-NEB is advantageous as it is a highly transferable method
applicable to many reactions in a wide range of materials. Compared to other meth-
ods, it is known to provide a good prediction for the MEP as well as the saddle point
which other methods struggle with. However, CI-NEB requires that the initial and
final states are known a priori, which can be troublesome when researching new and
complex multi-step reactions. One also has to be aware that CI-NEB finds the ‘local’
MEP; there may be a lower energy pathway that is inaccessible due to the choice of
initial pathway. The following is a brief overview of CI-NEB which has been used
throughout this thesis to extract reaction pathways and their activation energies.
The starting point for CI-NEB is the plain elastic band method (PEB), where
several images between an initial and a final state are connected by springs forming
a band whose energy can be optimised to provide the MEP. The images can be
tentatively constructed using a linear interpolation between the initial and final states.
However, one has to be careful that the linear interpolation could create a path that
is far from the MEP, e.g. the substitution of two lattice sites where atoms would
overlap in some of the image configurations causing unphysically high energies. A set
of springs is then connected between the images, and the force felt by each image is
thus:
FiP EB = Fi + Fspring,i , (2.55)
where Fi is the force due to the potential energies of each configuration, and Fspring,i
is the force due to the springs:
Fspring,i = ki+1 (ri+1 − ri )2 + ki (ri − ri−1 )2 . (2.56)
After creating the images and connecting them with springs, an objective function
Page 34
can then be defined as:
P P
X X k
S(r1 , r2 , . . . ) = V (ri ) + (ri − ri−1 )2 , (2.57)
i=0 i=1
2
where the first term is a sum of the potential energies of each of the images and the
second term is a sum of the potential energies caused by the springs connecting the
images. One can then use a suitable algorithm to minimize the objective function in
equation 2.57 with respect to the atomic coordinates of the images, thus minimising
the energy of the band and leading to an optimised path between the initial and final
configurations.
The PEB method suffers from two salient shortcomings: 1) the optimised band
can cut corners and; 2) the images can slide down the band leading to an uneven
spacing between them. These shortcomings are due to the image spring forces that
are not aligned along the direction of the band and the potential forces that do not
act perpendicular to the direction of the band.
The nudged elastic band method (NEB) attempts to fix the shortcomings in the
PEB method by only using the spring force components parallel to the band and
potential force components perpendicular to the band. These corrections ‘nudge’ the
band onto the true MEP. The force felt by each image on a NEB pathway is:
k
FiN EB = Fi⊥ + Fspring,i , (2.58)
where the first term is the potential force component perpendicular to the direction of
the band and the second term is the spring force component parallel to the direction
of the band. They can be obtained trivially during the calculation as:
Fi⊥ = −V (ri ) + V (ri ) · τî τî , (2.59)
k
Fspring,i = (ki+1 (|ri+1 − ri |) + ki (|ri − ri−1 |)) · τî , (2.60)
where τi is the direction of the band at image i and k is the force constant of the
Page 35
spring. Instead of using the objective function introduced for the PEB method, the
NEB method finds the MEP by minimising F N EB directly.
However, the activation energy obtained from the NEB MEP may be smaller than
the true value, as the springs will pull the highest energy images down toward a lower
energy. This issue can be fixed by removing the springs connecting the highest energy
image, and inverting the force felt by the atoms in this image to push the image onto
the saddle point:
CI
Fhighest = Fhighest − 2fhighest · τ̂highest τ̂highest . (2.61)
This manipulation allows the highest energy image to be pushed up the potential
energy surface to find the true first order saddle point, providing a more accurate
estimation of the activation energy. This correction combined with NEB is known as
climbing image nudged elastic band (CI-NEB). This method has been used throughout
this thesis to obtain accurate estimates of the activation energies of various reactions.
Page 36
Modelling Amorphous Silicon Dioxide and
3
the Si/SiO2 interface
3.1 Introduction
The ubiquity and technological importance of a-SiO2 has led to a huge experimental
and theoretical research effort to understand its atomic structure over the decades.
It is all too easy to forget that almost nothing was known about the atomistic na-
ture of amorphous solids until just under a century ago, with competing theories
proposing that amorphous solids were either supercooled liquids or composed of solid
crystallites. (Ultimately, amorphous materials were shown to not flow, definitively
determining that they are solids. 57 ) Initial X-ray diffraction (XRD) studies showed
Chapter 3. Modelling Amorphous Silicon Dioxide and the Si/SiO2 interface
that amorphous solids do not exhibit the sharp, intense spots which are character-
istic of a crystal, instead displaying a diffuse signal over the entire pattern, leading
researchers to propose that the underlying solid is made up of crystallites; in partic-
ular, cristoballite crystallites in the case of a-SiO2 . A revolutionary advance in our
understanding of the atomic nature of amorphous solids came with the introduction of
the continuous random network model (CRN) by Zachariasen in 1932. 58 He ruled out
crystallite theory by comparing the mechanical strengths of a-SiO2 and its crystalline
analogues, arguing that their similarities indicate that the interatomic forces hold-
ing a-SiO2 should be similar to those holding together crystalline SiO2 (c-SiO2 ). He
then defined a-SiO2 as a 3D structure with a lack of symmetry and containing similar
bonding to c-SiO2 . XRD studies of the atomic structures of the c-SiO2 variants reveal
that they are made up of SiO4 tetrahedra. It then follows that one expects a-SiO2
to be made up of SiO4 tetrahedra. Thus, Zachariasen proposed that the structure of
a-SiO2 is made up of corner-sharing SiO4 tetrahedra with a varying intertetrahedral
angle, as shown in Fig. 3.1.
Figure 3.1: Two interlinked SiO4 units illustrating the flexible Si-O-Si angle. Si atoms
are coloured in yellow while O atoms are coloured red.
Page 38
The main focus of this thesis concerns the properties of charge trapping defects
in amorphous SiO2 . In order to be confident that the results of the calculations
will be truly representative of the properties of defects in a-SiO2 , the initial a-SiO2
model must be consistent with our experimental perceptions of the atomic structure
of a-SiO2 discussed above. This chapter discusses the techniques used to model and
characterise the atomic structures of a-SiO2 studied in this thesis. Models of a-
SiO2 were created using classical potentials and a melt-quench technique, similar to
what has been reported in the literature. 38,59 However, the choice of potential was
constrained to charge variable potentials in order to be able to model Si in varying
oxidation states. This would, in principle, allow for the modelling of an Si/SiO2
system and the creation of point defects. The first part of this chapter discusses
the choice of interatomic potential to be used in the modelling of a-SiO2 structures.
This is followed by a technical discussion of the melt-quench method used to generate
a-SiO2 and the methods used to characterise the a-SiO2 models generated. Finally,
this chapter concludes by modifying and extending the melt-quench method to model
Si/SiO2 systems and characterise their structures in a similar manner to the a-SiO2
structures generated using the melt-quench method.
3.2 Choosing the Classical Potential
To create models of a-SiO2 , we will be using molecular dynamics simulations. This

requires a choice of potential which can accurately describe its structure. The clas-
sical potentials considered for modelling bulk a-SiO2 and the Si/SiO2 interface are
ReaxFF 51 and COMB. 45 To model both these systems accurately, the potential must
be able to reproduce the structures of Si, SiO2 , and the intermediate states which
can exist at the interface. In order to choose the potential, a library of test clus-
ters containing Si in varying oxidation states was chosen whose structures have been
studied experimentally or at higher levels of theory than classical potentials. Both
ReaxFF and COMB were then used to optimize the structures of these small silica-
Page 39
like clusters so that they can be evaluated against the literature. The clusters test
both ReaxFF and COMB’s ability to reproduce the structure of a cluster containing
Si in an intermediate oxidation state.
3.2.1 Small Silicon and Silica Test Clusters
Wang et al. studied molecules of Si3 Oy (y=1-6) as structures representative of those

which may exist during the oxidation of silicon, 60 representing the sequential oxida-
tion of Si3 . The clusters were studied using photoelectron spectroscopy and ab initio
calculations. Geometry optimizations were performed at the Hartree-Fock level using
a 6-31G* basis set. The resulting geometries were further optimised using second
order Møller-Plesset perturbation theory and a 6-311G* basis set. For y=1-3, each
additional oxygen atom bridge bonds between a Si–Si bond, producing cyclic struc-
tures (see Fig. 3.2). For y=4-6 the cluster begins to adopt a tetrahedral structure
with a central Si atom bonded to four oxygen atoms.
Figure 3.2: Atomic structures of Si3 Oy clusters as reported by Wang et al. 60 From
left to right and top to bottom, y=1-6.
To test the abilities of ReaxFF and COMB to model intermediate oxidation

states of Si, these clusters were optimized using each potential as implemented in
the LAMMPS code 61 and compared to Wang’s results. 60 Starting from their ideal
structures according to the study by Wang, the Si3 Oy clusters were optimized using
the ReaxFF and COMB potentials with a conjugate gradient optimiser. The bond
Page 40
lengths obtained with both potentials are in very good agreement with the results
obtained by Wang, as can be seen in table 3.1. Both potentials have a minimum
associated with these clusters apart from Si3 O6 , for which the COMB potential was
not able to find one. Instead, the central Si is bonded to three oxygens instead of
four and the terminal Si is bonded to two oxygens as seen in Fig. 3.3 . This clearly
does not match the structure in Fig. 3.2. These results give an early indication that
the COMB potential does not provide the right potential energy surface for a Si/SiO2
interface model.
Wang paper 60 ReaxFF COMB

Si3 O
Si-Si 2.30 2.32 2.46
Si-O 1.75 1.69 1.72
Si3 O2
Si-Si 2.41 2.16 2.15
Si-O (1) 1.67 1.66 1.81
Si-O (2) 1.67 1.67 1.65
Si3 O3
Si-O 1.69 1.66 1.73
Si3 O4
Si-O (1) 1.72 1.76 1.63
Si-O (2) 1.67 1.62 1.70
Si3 O5
Si-O (1) 1.52 1.45 1.47
Si-O (2) 1.67 1.60 1.66
Si-O (3) 1.68 1.68 1.60
Si-O (4) 1.66 1.77 1.73
Si-O (5) 1.73 1.60 1.50
Si3 O6
Si-O (1) 1.53 1.46 N/A
Si-O (2) 1.69 1.56 N/A
Si-O (3) 1.68 1.76 N/A
All bond lengths are in Å
Table 3.1: Bond lengths of the Si3 Oy clusters optimised using COMB and ReaxFF
and compared to the results by Wang et al. 60
The electronic and structural properties of neutral and charged Sin On clusters
with n=3-5 were calculated using DFT within the local density approximation by
Chelikowsky et al. 62 These small silica-like clusters were investigated to extract bond
stretching and angle bending forces. The pseudopotential was tested against previous
Page 41
Figure 3.3: Atomic structures of Si3 O6 optimized using COMB.
Figure 3.4: Atomic structures of Sin On clusters studied by Chelikowsky et al. 62
theoretical work, agreeing closely with previous results in the literature despite a
systematic underestimation of the Si–O bond length. The geometry of the Si3 O3
cluster agrees with the planar ring structure calculated by Wang, 60 but the geometries
of the Si4 O4 and Si5 O5 resemble a buckled ring, seen in Fig. 3.4.
Chelikowsky paper 62 ReaxFF COMB

Si4 O4
Si–O 1.65 1.63 1.65
Si–O–Si 127 158.2 111.5
O–Si–O 104 109.7 131.63
Si5 O5
Si–O (1) 1.65 1.54 1.65
Si–O (2) 1.65 1.62 1.73
Si–O–Si 122 149.9 123.5
O–Si–O 104 93.9 119.9
◦
All bond lengths given in Å; Bond angles in .
Table 3.2: Sin On clusters optimised using COMB and ReaxFF.
Table 3.2 shows that the Si–O bond lengths of the Si4 O4 cluster calculated with
both ReaxFF and COMB are in reasonable agreement with Chelikowsky’s results.
However, the bond angles are not reproduced correctly using these potentials. The
same conclusion can be drawn from the results of the Si5 O5 cluster; the bond lengths
Page 42
seem to be reproduced satisfactorily whereas the bond angles are not reproduced as
accurately. This leads one to believe that should one of these potentials be used
to model the interface, the structures which may exist at the interface may not be
reproduced with the correct bond angle. However, it is clear that these structures
are qualitatively reproducible with both potentials and that potential energy minimas
exist for these structures.
Chelikowsky et al. also studied the structure of small Si clusters (Sin : n=2-6) using
the same DFT method. 63 They were obtained from a random atomic configuration
at an initial temperature of 3000 K, followed by an anneal down to 300 K. The final
structure is obtained from a steepest-descent optimization of the geometry obtained
at 300 K. The ReaxFF and COMB clusters both optimize to the same qualitative
Figure 3.5: Atomic structures of Sin clusters, where n=2-6.
structure. As can be seen in table 3.3, both potentials generally overestimate the
Si–Si bond lengths and seem to have trouble describing the bond angle in Si3 , but
other than that they both have a similar description of small Si clusters.
Summary
This library of small Si and SiOx clusters show that both ReaxFF and COMB can
describe Si in intermediate oxidation states reasonably well. Potential energy minima
exist for almost all the clusters studied with Si in varying oxidation states. From
these calculations, it is proposed that adopting a scheme where these structures may
be obtained using a charge variable force-field followed by a further optimization
performed at a higher level of theory would lead to quantitatively accurate geometries
for structures which exist at the interface. As both potentials have performed similarly
Page 43
in describing these small clusters, further investigation is needed so as to choose the

force-field to model the Si/SiO2 interface.
Chelikowsky paper 63 ReaxFF COMB

Si2
Si-Si 2.23 2.14 2.30
Si3
Si-Si 2.16 2.20 2.28
Si-Si-Si 80 60 128
Si4
Si-Si (1) 2.31 2.23 2.37
Si-Si (2) 2.40 2.42 3.36
Si5
Si-Si (1) 2.28 2.40 2.55
Si-Si (2) 3.04 3.40 3.59
Si6
Si-Si (1) 2.35 2.44 2.43
Si-Si (2) 2.71 2.44 2.80
All bond lengths are given in Å
Table 3.3: Sin clusters optimised using COMB and ReaxFF.
3.2.2 NVE Molecular Dynamics using ReaxFF and COMB
The chosen potential will be used in molecular dynamics simulations. Thus, the
implementation of the potential in the LAMMPS code was investigated by running
NVE molecular dynamics. An NVE simulation fixes the number of atoms, the volume
and the total energy of the system. These preliminary simulations were also used to
adjust the time step to be used in further simulations.
Using the LAMMPS code, both potentials are used to run an NVE molecular
dynamics simulation on a 2x2x2 periodic supercell of β-cristobalite. The simulation
is run for 5 ps and the charge on each particle is calculated at every time-step. Total
energies calculated in an NVE simulation fluctuate, which is to be expected from using
numerical integration techniques for the evolution of the atoms’ motions. However, on
closer inspection of the COMB total energy graph, there appear to be discontinuities
at certain time-steps, as illustrated in Fig. 3.6a. These discontinuities do not exist
in the ReaxFF total energy curve. Upon further analysis it was found that these
Page 44
discontinuities stem entirely from the potential energy component of the total energy,
appearing at the same time step and with the same energy difference. The kinetic
energy does not suffer from this discontinuity. The charges calculated by COMB were
analysed and were shown to suffer a similar discontinuity at the same time step as
the energy discontinuities, shown in Fig. 3.6b. Analysis of the trajectory at the time
steps associated with these discontinuities shows a movement of oxygen atoms across
a boundary into the adjacent image. The implementation of the COMB potential in
the LAMMPS code has trouble dealing with periodic boundary conditions, leading
to discontinuities in charge which are the cause of the discontinuity in the total and
potential energy discontinuities. The ReaxFF implementation in LAMMPS does not
suffer from this problem. Due to ReaxFF’s ability to describe the silica and silicon
clusters and its proper implementation in LAMMPS, it was decided that the ReaxFF
potential would be used as the potential to construct models of a-SiO2 and the Si/SiO2
interface.
Due to the use of the Verlet algorithm in molecular dynamics, the total energy in an
NVE simulation fluctuates, as mentioned earlier. This is a well known issue associated
with the truncation of the Taylor expansion of the positions which introduces an
error in their calculation. 64 This truncation manifests itself as fluctuations in the
total energy and is dependent on the time step used in the simulation: the larger
the time step chosen, the greater the fluctuation. Choosing the time step is a trade-
off between minimising the total energy fluctuation and being able to achieve the
desired time of simulation. The root mean square fluctuation (RMSF) is a measure
of the deviation of the total energy against the average total energy and is calculated
according to equation 3.1. A generally accepted value for the RMSF limit is one part
in ten thousand. 65 The RMSF was calculated as:
T
1X
RM SF = (E(t) − Ē)2 , (3.1)
T t=1
where T is the time of the simulation, E(t) is the total energy at time t and Ē is the
average total energy. Two MD simulations were performed under NVE conditions
Page 45
Figure 3.6: a) Total energy and b) Average charges of the Si (black crosses, left hand
side axis) and O atoms (red crosses, right hand side axis), all plotted against time step
and extracted from the same NVE molecular dynamics simulation using the COMB
potential. Both the total energy and the average charges show discontinuities at the
same time step.
Page 46
for a 2x2x2 periodic supercell of β-cristobalite using the ReaxFF force field for 50 ps
each. The only difference between the two simulations was the time step chosen; one
simulation was run using a time step of 0.5 fs, the other was run using a time step of
0.1 fs. The RMSF for the 0.5 fs time step was 10.77 × 10−3 while it was 3.05 × 10−4
for the 0.1 fs time step. The RMSF value for the 0.5 fs time step is an order of
magnitude larger than the generally accepted limit while it is a little larger for the 0.1
fs. A time step smaller than 0.1 fs would greatly increase simulation times making it
infeasible to choose a smaller time step. A trade off has to be made between accuracy
and computational cost, therefore the time step that will be used in the molecular
dynamics simulations in this thesis is 0.1 fs.
3.2.3 Summary
Short molecular dynamics simulations were run to assess the computational stability
of the charge equilibrating potentials ReaxFF and COMB. The COMB potential
was found to show discontinuities in the species’ calculated charges which causes
discontinuities in the potential energy surface. In contrast, the ReaxFF potential
showed no discontinuities and was found to be stable for use in MD simulations.
In addition, a time step for the MD simulations was chosen for the remainder of
this thesis based on convergence of the root mean squared fluctuation below 0.0001.
Thus the potential that is used to model a-SiO2 throughout this thesis is the ReaxFF
potential.
3.3 Making Amorphous Silica Models
Amorphous silica is a metastable state of SiO2 which is defined by its disorder and lack
of well defined long-range geometrical properties. It is made up of SiO4 interlinked
tetrahedra with a flexible Si-O-Si angle. The same interatomic forces that act on
crystalline polymorphs of silica act on a-SiO2 , leading to a material with a fairly well
defined short-range but lack of long-range order.
Page 47
The short range order of a-SiO2 results in a number of geometrical properties

being described by a statistical distribution, as opposed to the discrete values found
in its crystalline analogues. As the SiO4 unit maintains its tetrahedral shape, the
O-Si-O angle has a sharp peak around 109.5◦ . The Si-O bond length also forms a
sharp peak at around 1.6 Å. The Si-O-Si angle is quite flexible, leading to a rather
broad distribution, between 120◦ and 180◦ , with the most frequent value occurring at
148◦ . 66,67
Amorphous silica is made experimentally by melting a crystalline polymorph of
silica followed by a rapid cooling known as a quench. Molecular dynamics can sim-
ulate this process, albeit with much faster quenching rates, and has previously been
used to investigate a-SiO2 using various interatomic potentials as discussed in the
introduction of this chapter. 59,68,69 In this section, the ReaxFF potential was used
to model amorphous silica within this general scheme which shall be referred to as
the melt-and-quench method. Starting from β-cristobalite, a crystalline polymorph
of SiO2 , a thermostat will be applied to raise the temperature of the system to form
a liquid melt which will then be rapidly cooled down to form an amorphous silica
model. The resulting structures are then characterised and compared to experiment.
All the structures are then further optimised at the DFT level using the CP2K code.
The electronic structures of the models are extracted and analysed to understand the
electronic properties of a-SiO2 .
3.3.1 Modelling a-SiO2 using the ReaxFF potential and the
Melt-and-Quench Method
The silica polymorph chosen to melt and quench was β-cristobalite. It was chosen for
its density (2.33 g cm−3 ), which is the closest of the silica polymorphs to the average
amorphous silica density of 2.2 g cm−3 and due to its cubic cell shape. A Berendsen
thermostat and barostat was applied to a periodic supercell of β-cristobalite contain-
ing 216 atoms. A time constant has to be chosen for the thermostat and barostat.
Initially, a number of NVT (constant number of atoms, volume and temperature)
Page 48
simulations were run on β-cristobalite at 300 K for 50 ps to adjust the time con-
stants of the Berendsen thermostat so that the temperature oscillated within 3% of
its desired value and so that the system was able to equilibrate within a few hundred
time steps. The time constant for the Berendsen barostat was adjusted in a similar
manner, by fixing the pressure at 1.0 bar and running a number of simulations using
different time constants. The time constant chosen for the Berendsen barostat was
the smallest value which did not cause large fluctuations in volume. These criteria
led to a time constant of 5 fs−1 for the thermostat and 1000.0 fs−1 for the barostat.
Figure 3.7: Histograms of the short-range geometrical properties of a-SiO2 from 320 a-
SiO2 models generated using the ReaxFF potential and the melt and quench method.
a) Si–O bond lengths. b) O–Si–O bond angles. c) Si–O–Si bond angles.
With the time step and the thermo- and baro-stat’s time constants chosen, sim-
ulations were then performed to heat the system up. Initially, each particle in the
216 atom system of β-cristobalite was assigned a random velocity from a Gaussian
distribution centred around 300 K and with a constraint that the sum of all velocities
result in the system’s temperature being 300 K. With the temperature maintained at
300 K and the pressure at 1 bar, a simulation was run for 20 ps to equilibrate the
Page 49
system. The temperature of the system was then linearly increased to 5000 K over
200 ps. By analysing the gradient of the potential energy curve as the temperature
increased and finding the point where the gradient changes, the melting point can
be identified. This is because bonds are broken as the material melts point which
increases the rate of potential energy change. The melting temperature was found at
4500 K and by 5000 K the system was found to be a liquid. The system was main-
tained at 5000 K for 1 ns, then cooled down at a rate of 6 K ps−1 to 0 K. The cooling
rate was chosen as the studies by Vollmayr et al. suggest that the positions of the
peaks in the short-range geometrical properties of the amorphous structure converge
at a cooling rate of 8 K ps−1 . 68 To study the geometrical statistics of a-SiO2 gen-
erated using this procedure, 320 models were generated and the resulting structures
were then optimised and analysed.
The average density of the a-SiO2 samples obtained were 2.16 g cm−3 and ranges
from 1.99 to 2.27 g cm−3 , which is in excellent agreement with experiment (2.2 g
cm−3 ). 70 The mean Si–O bond length is 1.58 Å, slightly shorter than the experi-
mental value of 1.61 Å. 67 The average Si-O-Si angle obtained was 154.75◦ which is
slightly larger than the experimental value of 148.3◦ . The O-Si-O angle obtained is
109.30◦ , which is in excellent agreement with the experimental value of 109.47◦ . The
graphs shown in Fig. 3.7 demonstrate the lack of order in the systems and shows
the distribution of bond lengths and angles, as opposed to the sharp peaks which are
characteristic of ordered crystalline systems. Interestingly, the extracted distributions
display a normal distribution.
Although it is useful and intuitive to characterise amorphous structures in terms
of their short-range geometrical properties, there are no experimental data available
against which one could directly compare. Unlike crystalline materials, the bond
lengths cannot be directly extracted by means of X-ray or neutron diffraction exper-
iments. This is due to the crystalline structure being interpreted using Bragg’s law,
which is entirely inappropriate for amorphous solids. Instead of the characteristic
bright spots, shown in Fig. 3.8a, and halos seen in X-ray and powder X-ray diffrac-
Page 50
Figure 3.8: Real space representation and spatial frequencies of crystalline and amor-
phous solids. a) The left image shows the atomic structure and electron density
of α-quartz and the right image shows the electron density’s Fourier transform. The
periodicity is obvious in real space and is reflected as ordered, discrete spatial frequen-
cies in the right image. b) The left image shows the atomic structure and electron
density of an a-SiO2 model, while the right image is the electron density’s spatial
Fourier transform. The real space structure is disordered, although it is ordered at
short-range. As a result, its Fourier transform shows no obvious order.
tion experiments on crystalline materials, the diffraction pattern of an amorphous

solid consists of a bunch of diffuse halos, similar to extracting a plane from Fig. 3.8b.
Information on the structure of an amorphous solid can still be extracted from X-ray
and neutron diffraction experiments by measuring what is known as the structure
factor, which is defined as: 71
Page 51
1 X
S(Q) = bj bk heiQ[rj −rk ] i, (3.2)
N j,k
where S is the structure factor, Q is the scattering vector, bj and bk are the scattering
lengths (X-ray of neutron) while rj and rk are the atomic coordinates of atoms j and
k. Assuming that the material is isotropic, which is normally valid for amorphous
materials, the vector [rj −rk ] can take on any orientation with respect to the scattering
vector. Therefore it is necessary to sum over all orientations in order to obtain the
orientational average:
π
sin Q[rj − rk ]
Z
iQ[rj −rk ] 1 2 sinφdφ
he i= e(iQ[rj −rk ] cos φ)2π[rj −rk ] = , (3.3)
4π[rj − rk ]2 φ=0 Q[rj − rk ]
so that:
1 X sin Q[rj − rk ]
S(Q) = bj bk . (3.4)
N j,k Q[rj − rk ]
Thus, the structure factor can be efficiently extracted from the atomic coordi-
nates of a calculation. It can also be obtained from X-ray and neutron diffraction
experiments, as it is related to diffraction as:
I(Q)
S(Q) = , (3.5)
Nf2
where I is the intensity of diffraction, N is the number of atoms and f is the total
atomic scattering factor. The structure factor is a 1-dimensional representation of the
3-dimensional amorphous structure, providing information on the short-, medium-
and long-range order. Fig. 3.9 shows the average total neutron structure factor
calculated from 320 a-SiO2 models which were made using the ReaxFF potential as
well as an experimentally measured signal. 72 The structure factor of the ReaxFF
models agrees very well with the experimental data, showing the same number of
peaks despite a few of the peaks being shifted slightly. The peaks agree well even
at large Q values, indicating that the medium- and long-range order of the models
matches that of experiment. These results indicate that the melt-and-quench method
Page 52
using the ReaxFF potential provides a-SiO2 models which are representative of the
experimental samples.
Figure 3.9: Total neutron structure factors from ReaxFF a-SiO2 models, before and
after DFT optimisation, and from neutron scattering experiments. The solid black
line shows the average total structure factor of a-SiO2 from 320 models made using
ReaxFF and the melt and quench method. The dashed blue line is the average total
structure factor of these ReaxFF models after they have been optimised using DFT.
The red circles are experimentally measured data. 72
3.3.2 DFT Optimisation of ReaxFF Structures
After obtaining the a-SiO2 structures using ReaxFF, they were optimised further and
their electronic structures were extracted using DFT as described in chapter 2.3. To
describe exchange and correlation, the non-local PBE0 TC LRC functional 34 was used
along with the auxiliary density matrix method. 73 This functional uses a truncated
Coulomb operator to calculate the Hartree energy for which a truncation radius of 2.0
Å was used. This functional was used as it contains a portion of Hartree-Fock exchange
which improves the description of the band gap as well as describing localised states
better than traditional local and gradient-corrected density approximations. Using
Page 53
the Gaussian Plane waves method, the Gaussian basis set used was a double-ζ basis
set with polarisation functions for both the Si and O atoms, while the plane wave
basis set was truncated at 5440 eV (400 Ry). All of the systems were optimised
individually using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. The
structures were considered optimised when the forces on the atoms were converged to
within 37 × 10−12 N (2.3 × 10−2 eV Å−1 ).
Figure 3.10: Histograms of the short-range geometrical properties of 320 DFT opti-
mised models of a-SiO2 . a) Si–O bond lengths. b) O–Si–O bond angles. c) Si–O–Si
bond angles.
The optimisation resulted in a slight change in the models’ atomic structures. As

can be seen from Fig. 3.10, the Si–O bond lengths show similarly shaped distributions
before and after optimisation, but the average bond length shifts up to 1.61 Å, in
much closer agreement with the experimental value of 1.62 Å. 74 Similarly, the O–Si–O
and O–Si–O bond angles, shown in Fig. 3.10, show similar shapes before and after
optimisation. The average of the O–Si–O angle after optimisation is 109.47◦ , very
similar to the unoptimised structure. However, the O–Si–O angle averages at 146◦ ,
very different from the unoptimised structure and much closer to experiment. The
Page 54
total structure factor of the DFT optimised structures is plotted in Fig. 3.9 along
with the original ReaxFF structures and experimental data. The structure factors
are very similar at shorter ranges (< 7.5 Å−1 ) before and after DFT optimisation,
with the optimised structures showing better agreement at longer ranges. The error
compared to the experimental data is reduced after optimisation, indicating that the
optimisation is a useful step in calculating the atomic structure of a-SiO2 .
Figure 3.11: Histogram of the atomic displacements of the ReaxFF a-SiO2 models
before and after DFT optimisation.
To obtain a measure of the difference between the ReaxFF structures before and
after DFT optimisation, their displacement fields have been calculated as:
1
D = [J3,1 · (RDF T − RReaxF F )◦2 ]◦ 2 , (3.6)
where D is a matrix in which element i is the displacement of the ith atom, J3,1 is
a 3 × 1 matrix of ones, ◦ indicates the Hadamard product, RReaxF F and RDF T are
Page 55
matrices of the atomic coordinates before and after DFT optimisation, respectively.
The displacements of all atoms are plotted in Fig. 3.11 for both Si and O atoms. The
displacements peak at well below 0.1 Å and extend to 0.4 Å and 0.8 Å for O and
Si atoms, respectively. These results indicate that the ReaxFF structures are very
close to a DFT minimum, but that the parameters for the Si atoms could be better
optimised.
Figure 3.12: Electronic densities of state from 320 models of a-SiO2 . Each state is
broadened by a Gaussian with a sigma value of 0.1 eV. The black curve is the total,
while the red and blue curves are the densities of state projected onto Si and O
atoms, respectively. The areas of the curve which are filled in to the baseline indicate
occupied states, while the curves which are not filled in indicate unoccupied states.
The states’ energies are normalised so that the Fermi level lies at 0.0 eV.
To assess the electronic structures of the a-SiO2 models, their electronic densities
of state were calculated. Fig. 3.12 shows a summation of these extracted densities
of state from the 320 models of a-SiO2 . Each state is broadened by a Gaussian
with a sigma value of 0.1 eV to simulate homogeneous broadening. In addition,
the densities of state were projected onto Si and O atoms in order to assess each
Page 56
atom type’s contribution to the total density of states. All states are aligned so that
the Fermi level lies at 0.0 eV and the occupied states lying below the Fermi level
are filled in to the base line. Immediately, one can see that a-SiO2 is an insulator
with a very wide band gap. The average band gap extracted from these calculations
is 8.1 eV and ranges from 7.1 to 8.4 eV. This is a slight underestimation of the
experimental SiO2 band gap, which has been measured optically as 8.8 eV. 75 The
traditional exchange-correlation functionals, such as the local density and gradient
corrected approximations, are known to underestimate the band gap of materials, with
the band gaps of these models calculated with PBE (a gradient corrected functional)
averaging at 6.4 eV . However, the inclusion of Hartree-Fock exchange is known to
improve the calculated band gap. Increasing the amount of Hartree-Fock exchange
increases the band gap, but it also changes the chemical properties of the material
in an unpredictable manner. In fact, in the literature many studies have used the
proportion of Hartree-Fock exchange as an effective tuning parameter to replicate
the band gap. However, in this thesis the amount of HF exchange will remain at
0.2, the value quoted in the original PBE0 TC LRC functional’s paper. It is also
clear to see that the top of the a-SiO2 valence band is dominated by contributions
from O atoms. In fact, these are contributions from O non-bonding ‘p’ states. The
bottom of the conduction band is made up of hybrid Si ‘sp’ states as well as non-
bonding O ‘p’ orbitals as well. An interesting property that emerges from these
calculations is that the extrema of the bands are not entirely delocalised, but are
instead partially localised at particular locations in these structures. The disorder of
the amorphous SiO2 structures leads to deviation from band theory, with the states
induced by disorder known as the band’s Urbach tails in the literature. 76 The ability
of these Urbach states to act as precursors to strong localisation is known as Anderson
77
localisation and these states in a-SiO2 are explored in detail in chapter 6.
Page 57
3.3.3 Summary
A molecular dynamics procedure is established to generate a-SiO2 structures by ad-

justing the time step, thermostat and barostat and following a melt and quench pro-
cedure. The ReaxFF potential is used to generate 320 a-SiO2 samples. The resulting
amorphous structures are characterised and found to agree very well with experiment;
the density is very close to the experimental value as are the total structure factors.
These structures were then optimised using DFT leading to an even better agreement
of the atomic structure with experiment. Their electronic structures were extracted
showing that the calculated band gaps are in reasonable agreement with experiment,
albeit a little too low. These results indicate that the use of molecular dynamics and
ReaxFF can provide an accurate description of the structure and disordered nature
of amorphous silica. Further optimisation using DFT is not necessary to improve
the description of the atomistic structure, but is necessary to extract the electronic
structure.
3.4 Constructing the Si/SiO2 Interface
Making a credible model of the Si/SiO2 interface has proved to be extremely chal-
lenging due to a lack of information of the microstructure of the interface. From
transmission electron microscope (TEM) images it is known that silicon evolves from
its crystalline elemental form into the totally disordered and fully oxidised amorphous
silica over a 5 Å abrupt transition region. 78,79 It is also known that Si–O–Si angles
are more compressed at the interface than in bulk silica. 80 X-ray reflectivity (XRR)
experiments reveal that the density of the SiO2 layer in the Si/SiO2 system is higher
than that of bulk SiO2 , measuring between 2.36 g cm−3 and 2.41 g cm−3 . X-ray pho-
toelectron spectroscopy (XPS) measurements again show that the transition region is
approximately 5 Å and also give a measure of the concentration of partially oxidised
silicon atoms across the interface. From this data, Si+1 and Si+2 species appear right
at the interface on the bulk silicon side. On going further into the transition region
Page 58
the Si+1 and Si+2 species subside and Si+3 appears. The XPS data indicates that the
ratio of the partially oxidised silicon atoms is 1:2:3 for Si+1 :Si+2 :Si+3 . 81 This experi-
mental data sets out the criteria for an ideal model of the interface. In summary, the
model must satisfy the following criteria:
• Contain an a-SiO2 layer which matches the geometrical description of bulk amor-
phous silica past the interface and a substrate which matches elemental Si before
the interface.
• The model must contain a transition layer which neither matches bulk silicon
nor bulk silica.
• Contain an a-SiO2 bond length distribution which is slightly longer than that
of bulk a-SiO2 in the interface region.
• Possess an Si-O-Si angle distribution which is more compressed with respect to

bulk silica.
• The density of the amorphous silica layer should be higher than that of bulk
amorphous silica.
• The correct distribution of silicon in intermediate oxidation states across the

interface, obtained from XPS experiments.
3.4.1 Melt and Quench in situ
A model of the Si/SiO2 interface was prepared starting from a 3x3x4 periodic supercell
of crystalline silicon. Oxygen atoms were then inserted in the 001 direction between
silicon atoms over 4 layers to make a layer of silica. This initial stack is shown in
Fig. 3.13. Inserting oxygen atoms between silicon atoms results in the SiO2 layer
being very similar to the crystal structure of β-cristobalite; however, the lateral cell
vectors of this silica layer, which are the same as the silicon unit cell (5.28 Å), are
much more compressed than the β-cristobalite cell vectors (7.12 Å). To remove some
of the stress that comes from compressing this β-cristobalite layer, the silica layer
Page 59
Figure 3.13: Silicon supercell with oxygen atoms inserted to form silica layer.
Note that the atomic structure of the SiO2 layer is almost identical to that of
β-cristobalite.
was then extended in the direction of the interface so that the density of the silica
layer corresponded to that of β-cristobalite. Random velocities from a Gaussian
distribution were then assigned to all atoms so that the temperature of the system was
300 K and a molecular dynamics simulation was then run for 20 ps. The temperature
and pressure were maintained at 300 K and 1 bar using a Berendsen thermostat and
barostat whose time constants are obtained from section 3.3.1. The temperature of
the silicon substrate was then fixed at 300 K while the temperature of the SiO2 layer
was brought up to 5000 K to form a liquid melt of the silica layer. The temperature of
the silica layer was then fixed at 5000 K, while the temperature of the silicon substrate
remained fixed at 300K, and a simulation was run for 1ns. The silica layer was then
cooled down to 300 K at a rate of 6 K ps−1 .
The resulting structure is shown in Fig. 3.14. Upon visual inspection, it appears
that many oxygen atoms have moved out of the silica region and penetrated the silicon
substrate. Due to the high temperature required by the ReaxFF potential to melt
SiO2 (see section 3.3.1), the oxygen atoms have a very high velocity which allows
them to move through the silicon substrate easily. The system prepared this way also
Page 60
Figure 3.14: Si/SiO2 stack after in situ melt and quench of the SiO2 layer.
SiO2 layer does not form an abrupt interface with oxygen atoms penetrating deep
into the silicon substrate. In addition, there is a large number of coordination
defects both in the Si substrate and the SiO2 layer.
generates a wide range of coordination defects, including triply coordinated oxygen

atoms, singly coordinated oxygen atoms and overcoordinated silicon atoms. Although
this analysis is not particularly extensive, it already indicates that this method is
unsuitable for the creation of an ideal Si/SiO2 interface system. For the purposes of
investigating the defects which may exist at the interface, the stack prepared in this
way cannot be used as there are far too many coordination defects which distort the
physics and chemistry of the ideal Si/SiO2 interface.
3.4.2 a-SiO2 Layer Strained to Silicon Lateral Dimensions
To overcome the limitations of the melt and quench method and use it to model the
Si/SiO2 interface, the procedure needs to be improved in a way that would eliminate
coordination defects. A similar method was then investigated in which the amorphous
silica layer is made in an independent system so that it can be inserted between the
silicon substrate. This a-SiO2 layer is made starting from a periodic model of β-
cristobalite whose lateral cell dimensions have been fixed to that of silicon (in this
case 5.28 Å for both the a and b cell parameters). The size of the box in the 001
direction is allowed to move by using a Berendsen barostat that adjusts the pressure of
the system by scaling only the z dimension of the simulation box. Random velocities
Page 61
are assigned to all atoms in the system so that the system’s temperature is 300 K
and a molecular dynamics simulation with 3D periodic boundary conditions is run
for 20 ps. The temperature and pressure are again controlled using the Berendsen
thermostat and barostat. The temperature of the system is then raised linearly to
5000 K and the temperature is then maintained at 5000 K for 1 ns. The system is
then cooled down to 300 K at a rate of 6 K ps−1 and minimised to give an amorphous
silica sample sandwiched between silicon.
Figure 3.15: Si/SiO2 interface made with a-SiO2 layer placed on top of silicon.
The resulting amorphous silica structure is a continuous random network which

contains no defects. It has the same geometrical distributions as the amorphous silica
described in section 3.3.1. The amorphous silica layer was then heated to 600 K
and equilibrated at that temperature for 20 ps. In a separate system, a 3 × 3 × 3
supercell of crystalline silicon was heated to 600 K and equilibrated also for 20 ps.
The atomic coordinates and velocities of the silica system were then copied into the
crystalline silicon system, with the coordinates of the silica atoms translated in the
001 direction so that the silica layer was placed over the silicon substrate separated by
a 1.6 Å vacuum gap on either side of the silica layer. This vacuum gap is introduced
so that the atoms at the surfaces of the silicon substrate and the silica layer will have
freedom to try and position themselves in energetically favourable positions and is
also the reason for heating the system up to 600 K. The box will be allowed to relax
Page 62
in the dimension perpendicular to the interface and it is hoped that new bonds will
be formed at the interface. The 1.6 Å gap was chosen after several simulations were
run with vacuum gaps ranging from 1.6 Å to 2.5 Å. The simulations which contained
the 1.6 Å vacuum gap produced the least overcoordinated oxygens while the larger
vacuum gaps resulted in unphysical surface terminations for both the Si and SiO2
layers. A molecular dynamics simulation is then run on this Si/vacuum/SiO2 stack
under periodic boundary conditions for 1ns with the temperature maintained at 600K
and the pressure maintained at 1 bar. The system was then cooled down to 300 K
at a rate of 6 K ps−1 and Fig. 3.15 shows the interface formed using this method.
Immediately upon visual inspection, it can be seen that the interface formed is abrupt
and amorphous.
The geometrical properties of this interface model are presented in Fig. 3.16.
The Si–O bond length distribution peaks at around 1.65 Å. This is slightly longer
than that of bulk silica and is to be expected according to infra-red absorption (IR)
experiments. 2 However, the bond lengths range from 1.55 Å to almost 2.0 Å in contrast
to IR experiments which do not find bond lengths longer than 1.75 Å. The mean
O–Si–O angle is 109.22◦ , which is to be expected as the SiO4 units maintain their rigid
tetrahedral shape. As can be seen in Fig. 3.16, the Si–O–Si angle peak has shifted
to the left with respect to the peak in bulk amorphous silica (see Fig. 3.7). The
mean Si–O–Si angle from this interface model is 136.60◦ , in excellent agreement with
the conclusion reached by Hirose et al. (drawn from XPS experiments and molecular
orbital calculations) that the Si-O-Si angle is reduced to 135◦ . 82 The oxidation state
of silicon across the interface is also shown in Fig. 3.16. It shows that Si+1 and Si+2
appear first at the interface while Si+3 appears deeper into the interface. The profile
shows a gradual change over a 5 Å region from elemental silicon (Si0 ) to fully oxidised
silicon (Si+4 ), qualitatively similar to what is seen in XPS experiments. The ratio of
partially oxidised silicon atoms in this model is calculated as 1:0.4:0.3 (Si+1 :Si+2 :+3 ).
This ratio does not match XPS 83 data which shows a ratio of Si+1 :Si+2 :+3 of 1:2:3.
The density of the silica layer is 2.37 g cm−3 which is also in good agreement with
Page 63
Figure 3.16: Geometrical properties of Si/SiO2 interface made with a-SiO2 layer placed
on top of silicon. a) The Si–O bond length distribution in the stack. b) The Si–Si
nearest neighbour distance in the stock. Note that it displays a bimodal distribution;
one for the bonds in the Si substrate, the other for the Si nearest neighbours in the
a-SiO2 layer. c) The Si–O–Si bond angle distribution. d) The Si–O bond length
distribution as a function of its position along the z axis. e) The Si–Si distances as
a function of position along the z axis. f) The ReaxFF charges calculated for Si as a
function of position on the z axis.
experiment. 81
3.4.3 Summary
Using this method, one can reproduce bond angle statistics, density of amorphous sil-
ica and the qualitative profile of silicon oxidation states across the interface. However,
the bond length distribution contains lengths which are far too long. The concentra-
tion of silicon in intermediate oxidation states in the transition region does not match
the experimental ratio seen by XPS experiments. This needs to be corrected if the
model is to be used to study defects at this interface as the defect’s properties will be
Page 64
strongly influenced by the chemistry. The wrong chemistry at the interface will lead
to a wrong description of the defect.
3.5 Conclusions
An interatomic potential known as ReaxFF was chosen to model amorphous silicon

dioxide. The potential allows the charges of the individual atoms to be calculated
based on the system’s instantaneous geometry. ReaxFF was used to construct 320
models of a-SiO2 using a melt and quench method showed excellent agreement with
experimental data. The melt and quench method was extended to the Si/SiO2 system.
However, it requires additional modification before it can be used to confidently study
the atomic structure of the Si/SiO2 stack.
Page 65
Defects in Bulk SiO2
4
4.1 Introduction
Point defects can strongly influence a material’s physical and chemical properties.
The properties of silicon dioxide have long been exploited in a number of technological
applications. In particular, SiO2 is used in electronic devices due to its extremely high
quality interface with the semiconducting Si substrate – high quality in this context
means very few defects at the interface between Si and SiO2 . Defects in SiO2 have
been implicated in many electronic device reliability issues, such as random telegraph
noise (RTN), bias temperature instability (BTI), and stress-induced leakage current, 11
as well as in degrading the device performance by reducing the channel mobility. 84 It
is thought that defects in the oxide can capture and emit the available carriers in the
Chapter 4. Defects in Bulk SiO2
semiconductor channel of electronic devices, leading to a strong electrostatic effect on

the carriers flowing through the device.
This chapter serves two purposes. Firstly, it presents a brief overview of defects in
bulk silica, particularly those that have been implicated in electronic device reliability
issues. Secondly, a number of defects that have been studied in the literature have
been re-calculated to assess the validity and accuracy of the computational methods
used in this thesis, discussed in detail in chapters 2 and 3.3.1. The defects studied
in this chapter are the oxygen vacancy and its charged variants. Initially, the elec-
tronic structure of oxygen vacancies in α-quartz are investigated as it is instructive to
investigate defects in a crystalline environment with minimal structural disorder. In
addition, many crystalline point defects are thought to have analogues in amorphous
silica. 85 The positively charged vacancy, a variant of which is widely known as the
E0 centre, and the negatively charged vacancies are then studied. Finally, the oxygen
vacancy and its charged configurations are studied in a-SiO2 .
4.2 Quartz
4.2.1 Neutral Oxygen Vacancy
The neutral oxygen vacancy – also known as oxygen deficient centre I (ODC I) –
has been modelled previously using various methods, including ab initio and semi-
empirical calculations, and modelled within an embedded cluster scheme and periodic
boundary conditions. Blöchl calculated the structure of the neutral vacancy in a 72
atom periodic cell of α-quartz. 86 His calculations were performed using DFT with
the exchange and correlation described within the generalised gradient approxima-
tion (GGA). Sulimov et al. 87 also modelled the neutral oxygen vacancy in quartz
using the embedded cluster method 88 which allowed them to model larger systems
than if they were to use quantum mechanical methods on their own. In Sulimov’s
investigation, a cell of quartz containing 9270 classical atoms and up to 67 quantum
atoms was used initially. An oxygen atom was removed from the region which is
Page 67
calculated quantum mechanically to form the oxygen vacancy. Various calculations

were performed with the quantum mechanical region varying in size from a cluster
containing two silicon atoms (and four oxygen atoms) to eighteen silicon atoms (forty-
nine oxygen atoms). The quantum mechanical region was allowed to relax while the
classical region remained fixed and the structure of the neutral oxygen vacancy was
elucidated.
Results from both Sulimov and Blöchl’s calculations provide a similar picture on
the structure of the neutral oxygen vacancy in α-quartz. Removal of the oxygen atom
from between the two silicon atoms leaves behind a dangling bond on each of the
silicons. The silicon atoms, initially ≈ 3.2 Å apart, then displace toward each other
to form a fairly strong Si–Si bond which measures 2.4 Å. Sulimov also investigated
the displacement field of the neutral oxygen vacancy which was found to be very long
ranged indeed, extending to as far as 13 Å away from the vacancy site. The long range
network relaxation has been shown to be asymmetric due to the absence of inversion
symmetry in quartz. 87
Figure 4.1: a) Crystal structure of α-quartz and b) asymmetry of an SiO4 tetrahedron

in α-quartz.
In this investigation, the neutral oxygen vacancy in α-quartz was modelled using
periodic DFT as described in chapter 2. Initially, a 3 × 3 × 3 cell of α-quartz con-
taining 243 atoms was constructed using coordinates obtained from X-ray diffraction
studies. 89 To describe the electrons around the Si and O atoms, the GTH pseu-
dopotentials 90 and double-ζ basis sets (including polarisation functions) 91 were ap-
Page 68
plied to the atomic coordinates. The exchange and correlation were described by the
PBE0 TC LRC functional 30 with a cutoff radius of 2.0 Å for the truncated Coulomb
operator. The electronic structure was then calculated within the DFT formalism
using the CP2K code and the system’s geometry was optimised with respect to its
energy so that the maximum force per atom in the system is less than 37.1 pN.
Optimisation of this structure lead to the peculiar α-quartz structure; for each SiO4
tetrahedron, there are two long (1.61 Å) and two short bonds (1.60 Å). Fig. 4.1 shows
both the α-quartz structure and a single SiO4 tetrahedron with the short and long
bonds highlighted.
To model the neutral oxygen vacancy, an oxygen atom was then removed ran-
domly from the optimised α-quartz cell, leaving behind a cell containing 242 atoms.
All oxygen atoms have one short and one long bond to silicon resulting in the same
chemical environment, so it does not matter which oxygen atom is removed. The
system was then optimised, resulting in the structure of the neutral oxygen vacancy,
which is shown along with its highest molecular orbital (HOMO) in Fig. 4.2. Opti-
mising the vacancy displaced the silicon atoms surrounding the vacancy toward each
other, forming a Si–Si bond of 2.44 Å. This is in excellent agreement with the bond
length calculated by Blöchl in a smaller periodic model containing only 72 atoms and
calculated using the PBE functional using a plane wave code. The newly formed bond
is a σ bonding state which is high in Si ‘p’-orbital character. The one-electron defect
level associated with the vacancy is a deep state, lying 0.8 eV above the valence band,
shown in the calculated electronic density of states in Fig. 4.2.
86
This investigation Blöchl paper
SiL -SiS 2.439 2.445
SiL -O(3+) 3.884 3.893
SiS -O 1.653 1.642
SiL -O 1.648 1.641
All bond lengths are measured in Å
Table 4.1: Comparison of neutral oxygen vacancy geometries in different models
A comparison of this geometry with Blöchl’s model is presented in table 4.1,

showing excellent agreement. The displacement field has also been calculated and
Page 69
Figure 4.2: Atomic and electronic structure of neutral oxygen vacancy in α-quartz
The graph shows the electronic density of states for a single NOV in a 242 atom cell
of α-quartz. An occupied state breaks off the top of the valence band and is highly
localised at the vacancy, as seen in the inset which visualises this state as an
isosurface (with an isovalue of 0.09) alongside its atomic structure. The Si atoms
surrounding the vacancy displace towards each other forming a Si–Si bond.
Page 70
is shown in Fig. 4.3. It is interesting to see that although the local structure of
the neutral oxygen vacancy had converged in Blöchl’s smaller periodic model, in
these calculations, containing 242 atoms and a cell 5 Å bigger in each dimension, the
displacement induced by this defect persists right up to the edge of the simulation
box. It would be interesting to investigate whether the displacement field is even more
long ranged than this or whether the long ranged displacement is an artefact of some
parameter of the calculation.
4.2.2 Positive Oxygen Vacancy
The positive oxygen vacancy (POV) has been studied extensively both experimentally
and theoretically as it is implicated in various device reliability issues. 84 In particular,
ab initio calculations of the POV have been used to explain the origins of 1 /f noise
and thermally stimulated current in MOS devices, 84,92,93 with suggestions that these
Figure 4.3: Displacement field of neutral oxygen vacancy in quartz.
Page 71
issues are caused by thermally activated capture and emission of holes at neutral
oxygen vacancy sites. A similar mechanism has been proposed for explaining NBTI;
however, the atomic nature of the hole traps involved has not been fully resolved. 94–96
Various structures have been proposed for the POV, assuming that it is formed
directly from the neutral oxygen vacancy. A hole is introduced into the neutral oxygen
vacancy which localises at that site giving the structures of the POV. The most
prominent of the positive oxygen vacancies are categorised into two kinds: The first
is a “dimer” vacancy configuration, with an analogue known as the E0δ centre in a-
SiO2 , 86,97 while the second is a centre in which one of the Si atoms relaxes through
the plane of its oxygen neighbours forming a ‘puckered’ configuration known as the
E01 centre in α-quartz with its analogue known as the E0γ centre in a-SiO2 .
Blöchl investigated the positive oxygen vacancy using the same calculation details
as for his neutral oxygen vacancy calculations discussed earlier in this chapter. After
introducing a hole into the optimised NOV in quartz, the atomic coordinates are re-
laxed resulting in the dimer configuration of the positive oxygen vacancy. Mysovsky
et al. also calculated the structure of the positive oxygen vacancy using the em-
bedded cluster method. 98 The quantum mechanical cluster used in this study was a
Si11 O33 cluster, while the whole system (quantum mechanical cluster embedded into
classically polarisable system) contained 9270 atoms. The electronic structure of this
quantum mechanical cluster was calculated using DFT with the exchange-correlation
described using the GGA functional. The orbitals were described using a local orbital
basis set which consisted of norm-conserving pseudopotentials for the core orbitals
and a double-ζ basis set for the valence orbitals. The structure of the dimer configu-
ration was then obtained by putting a hole into the optimised neutral oxygen vacancy
structure and optimising the atomic coordinates.
Both Mysovsky and Blöchl’s studies reach similar conclusions, the dimer configu-
ration is highly similar to the structure of the neutral oxygen vacancy but the Si–Si
bond is relatively weak. Blöchl finds that the vacancy bond expands by 23%, so that
the Si–Si bond is comparable in length to an oxygen bridge in quartz.
Page 72
The puckered configuration is a relaxation of one of the silicon atoms through the
plane of its oxygen neighbours, forming a bond with an oxygen atom behind it. This
happens due to the asymmetry of the tetrahedral unit, such that there is always an
oxygen atom behind the silicon atoms associated with the long Si-O bonds while there
is no oxygen atom behind the silicon atoms associated with the short Si-O bonds.
Here, both the dimer and puckered configurations of the POV in quartz were
investigated. The dimer configuration was obtained spontaneously upon the addition
of a hole into an NOV structure. However, to obtain the puckered configuration,
a perturbation to the structure is required to overcome the energetic barrier to its
formation. Despite there being a barrier, the puckered configuration was found to be
more thermodynamically favourable than the dimer configuration.
Figure 4.4: Atomic and electronic structure of the dimer configuration of the positive
oxygen vacancy in α-quartz. The graph shows the electronic density of states with
the filled in section indicating the occupied states. The densities in the positive y-axis
are from α spin states while the β-spin states occupy the negative y-axis. The highest
occupied state in the α-spin and lowest unoccupied state in the β-spin are visualised
in the insets along with the atomic structure of the dimer configuration.
Page 73
Dimer configuration
Using the same calculation parameters described in chapter 4.2.1, the positive oxygen
vacancy was investigated. Starting from the neutral oxygen vacancy structure, a hole
was added to the system and its geometry was optimised. The atomic structure is
presented in the insets of Fig. 4.4 (note that both insets show the same atomistic
structure). Adding a hole into the system displaces the two Si atoms surrounding
Figure 4.5: Displacement of the atoms from their lattice sites when the POV is
introduced. The black squares are Si atoms and the red squares are O atoms.
the vacancy from their position in the bulk. A longer range distortion can also be
seen in the displacement field shown in Fig. 4.5, although it is not as strong as that
of the neutral oxygen vacancy. The Si–Si bond in the dimer configuration is 2.93 Å,
showing that it is significantly weaker than that in the neutral oxygen vacancy. This
can be understood as an electron taken out of the stabilising bonding orbital, leading
to a weaker bond between the two Si atoms and causing them to displace away from
each other. This idea can be confirmed by visualising the highest occupied and lowest
Page 74
86 98
This investigation Blöchl paper Mysovsky paper
SiL -SiS 2.932 3.011 2.890
SiL -O(3+) 3.608 3.554
SiS -O 1.603 1.600
SiL -O 1.598 1.597
Table 4.2: Comparison of positive oxygen vacancy dimer configuration geometries in

different models
unoccupied states in the α and β spin channels, respectively. They can both be seen
in the insets of Fig. 4.4. The HOMO in the α channel shows a bonding orbital
between the two Si atoms, while the LUMO of the β channel, the state from which an
electron was removed, is clearly localised on just one of the Si atoms. The relaxation
of the dimer configuration moves an occupied defect state further into the band gap,
so that it lies 1.2 eV above the valence band, shown in Fig. 4.4. Calculating the spin
moment reveals that the spin density is shared between the Si atoms surrounding the
vacancy, although it is shared asymmetrically. We can categorise the two Si atoms
based on whether the oxygen removed formed a long bond or a short bond, which we
will refer to as SiL and SiS , respectively. The spin moment of SiL , the Si on which the
hole is localised (see LUMO in β-channel in Fig. 4.4) is 0.36 while it is slightly higher
on SiS at 0.41. We also find that SiL has a slightly more positive Mulliken charge
of +0.96|e| compared to the +0.94|e| of SiS . This follows from the spin moment
calculation, as the unpaired spin is due to an electron which is less localised on SiL
and therefore that centre is more positively charged. However, both Si atoms are
much less positively charged than the average charge of +1.09|e| for the remaining
bulk Si atoms. The geometrical properties of this model and those obtained from
the literature are presented in table 4.2. The structure calculated in this periodic
model containing 242 atoms agrees excellently with previous results, suggesting that
the dimer configuration is well described in cells containing 72 atoms.
Page 75
Figure 4.6: Atomic and electronic structure of the puckered configuration of the
positively charged oxygen vacancy in α-quartz. The graph shows the electronic density
of states, with the positive and negative sections of the y-axis showing the density
of α and β electron states, respectively. The HOMO of the α electrons and LUMO
of the β electrons are shown in the inset, with the states visualised as red and blue
polyhedra for the occupied state, while they are green and orange for the unoccupied
state.
Puckered Configuration
As discussed earlier, the puckered configuration of the positive oxygen vacancy is a

relaxation of the silicon atom associated with the long Si–O bond through the plane
of its oxygen neighbours to benefit from a weak interaction with the oxygen behind it.
There is a potential energy barrier between the dimer configuration and the puckered
configuration which prevents the puckered configuration from being spontaneously
generated by hole trapping at the neutral oxygen vacancy structure. In order to model
the puckered configuration and overcome this energetic barrier, a Si atom in the dimer
configuration of the positively charged vacancy was manually translated through the
plane of its oxygen neighbours, followed by an optimisation using the calculation
Page 76
details presented in chapter 4.2.1. The structure of the puckered configuration is

presented in the insets of Fig. 4.6. The optimisation resulted in a weak bond forming
between the Si which has relaxed through the plane of its neighbours, which shall be
referred to as SiB , and there is now no bonding state at all between the two Si atoms.
The Si atom which hasn’t moved through the plane of its neighbours and which is
now 3-coordinated is referred to as the 3c-Si. The calculated spin moments reveal
that there is almost no spin density on SiB which has a spin moment of 1.7 × 10−3 ,
but is instead almost entirely localised on the 3c-Si with a spin moment of 0.85. The
3c-Si is also much less positive with a Mulliken charge of +0.86|e|, while the SiB has
a Mulliken charge of +1.14|e|, comparable to the remaining bulk Si atoms’ average of
+1.10|e|, but slightly more positive. The geometry of this configuration is presented
in table 4.3 along with a comparison with previous results. The geometry calculated
86 98
This investigation Blöchl paper Mysovsky paper
SiL -SiS 4.328 4.358 4.320
SiL -O(3+) 1.881 1.852 1.910
O(3+)-Si 1.812 1.793 1.790
SiS -O 1.650 1.643
SiL -O 1.602 1.599
Table 4.3: Geometry of the puckered configuration of the positive oxygen vacancy in
different studies.
agrees excellently with previously published results. The density of states of this
puckered configuration shows that the occupied and unoccupied states associated
with this defect sit 2.8 and 6.9 eV above the valence band, respectively, revealing
that this configuration acts as an effective hole trap. The hyperfine interactions of
this unpaired electron’s spin with the nuclei around the defect have been calculated.
The hyperfine interaction depends upon the spin density at the nucleus, therefore the
use of a pseudopotential in place of the core electrons is not applicable. An accurate
description of the electrons at the nucleus is required, so it is therefore necessary to
describe the core electrons with a basis set. In this calculation, a 6-311G** basis
set 99,100 was used on the Si and O atoms and the defect structure was reoptimised
Page 77
86 98
This investigation Blochl paper Mysovsky paper Experiment 101
SiS -46.45 -45.4 -49.93 -45.31
-40.31 -39.6 -43.6 -39.07
-40.30 -39.6 -43.55 -39.06
SiS1 -1.18 -1.4 -1.08 -0.92
-0.97 -1.2 -0.92 -0.75
-0.97 -1.2 -0.91 -0.75
SiS2 -1.45 -1.5 -1.16 -0.98
-1.22 -1.3 -0.99 -0.79
-1.21 -1.3 -0.98 -0.79
All results are measured in mT
Table 4.4: Principal values of the hyperfine tensor for the puckered configuration in
this model and from previous theoretical and experimental studies. The experimental
results are for what is known as the E01 centre in α-quartz.
within this basis. This basis led to minimal geometric changes; none of the bond
lengths changed by more than 0.05 Å. The hyperfine coupling tensor was calculated
and diagonalised to extract the principal values and compare to previous theoretical
values and experiment. This is presented in table 4.4. The calculated hyperfine
values are in excellent agreement with both previous theoretical calculations and
experimental results for what is known as the E01 centre in α-quartz.
4.2.3 Negative Oxygen Vacancy
The NOV introduces an unoccupied state which sits below the bottom of the α-quartz
conduction band. We investigate whether an electron can be localised in this state.
Starting from the optimised configuration of the NOV, we added an electron to the
system and optimised its geometry. This resulted in weak atomic displacements in the
vicinity of the defect, with the Si–Si bond length remaining at 2.4 Å. However, analysis
of the electronic structure reveals that the electron has indeed localized on the Si–Si
bond, as shown in the inset of Fig. 4.7. The Mulliken charges of the two Si atoms
surrounding the vacancy are +0.62 and +0.74 |e|, clearly much more negative than the
average of +1.08 |e| for the remaining bulk Si atoms in α-quartz. The extra electron
occupies a state which sits 3.0 eV below the bottom of the α-quartz conduction band,
as shown in Fig. 4.7. The negative oxygen vacancy’s atomic structure is very similar
Page 78
Figure 4.7: Atomic and electronic structure of the negatively charged oxygen vacancy
in α-quartz. The highest occupied state in the α channel of electrons is visualised
in the inset using an iso-value of 0.1. The defect levels are plotted against the Si–Si
length around the negative oxygen vacancy.
to that of the NOV and the dimer configuration of the positively charged oxygen
vacancy, but the extra electron is localised in an anti-bonding state between the two
Si atoms surrounding the vacancy.
4.2.4 Summary
The oxygen vacancy in α-quartz has been calculated in a number of charge states.
Our results for the neutral and positive charge states of the vacancy agree well with
previous theoretical and experimental studies in the literature. The unoccupied state
of the neutral oxygen vacancy can be occupied by an electron to make the negatively
charged oxygen vacancy. All of the vacancy’s charged variants are accompanied by
strong lattice distortions.
Page 79
4.3 Amorphous Silicon Dioxide
We now turn to exploring the oxygen vacancy defects in amorphous SiO2 (a-SiO2 ).
Due to the similarity in short-range order between α-quartz and a-SiO2 , one can
initially explore analogous structures to those presented in the previous section. In
fact, the structures of vacancies have indeed been shown to be analogous in a-SiO2 as
well as for other defects such as interstitials. 102,103 In addition, the vacancies in a-SiO2
have been well studied due to the implication of the positively charged vacancy known
as the E0 centre in various reliability issues. 9,84,104,105
Here, 20 of the a-SiO2 models generated in chapter 3 were used to study the
charged variants of the oxygen vacancy. Unlike crystalline α-quartz where only one
model was examined, the vacancy in a-SiO2 was calculated in 20 different models
due to the structural disorder of the amorphous samples. This leads to the defect’s
properties being spread over a distribution rather than being discrete, requiring it to
be studied in a number of different configurations. Initially, the neutral vacancy was
studied, followed by its positive and negatively charged variants.
4.3.1 Neutral Oxygen Vacancy
The neutral oxygen vacancy (NOV) in a-SiO2 was modelled using periodic DFT as
described in chapter 2 and using 20 of the a-SiO2 models constructed in chapter 3.
The electrons are again described using the GTH pseudopotentials 90 and double-ζ
basis sets (including polarisation functions) 91 while the PBE0 TC LRC functional 30
exchange-correlation functional 30 was applied with a cutoff radius of 2.0 Å for the
Coulomb operator. In each of the 20 models, an O atom was randomly removed from
the structure followed by a geometry optimisation. In all of the models, this resulted
in the Si atoms surrounding the vacancy displacing towards each other to form an
Si–Si bond, analogous to the NOV in α-quartz. However, unlike the NOV in α-quartz,
in a-SiO2 the Si–Si bond length is a distribution which averages at 2.53 Å and has
a width of 0.69 Å. The distribution of bond lengths can be seen in Fig. 4.8. These
Page 80
Figure 4.8: Atomic structure and electronic structure plotted against Si–Si bond
of the neutral oxygen vacancy in a-SiO2 . The inset shows the HOMO of a single
configuration of the NOV. The graph shows the one-electron states associated with
the defect plotted against its Si–Si bond length. The shaded in regions indicate the
average valence and conduction band edges in a-SiO2 , while the red and green circles
represent the occupied and unoccupied defect states, respectively.
results are in good agreement with the results from embedded cluster calculations by
Mukhopadhyay et al. 59,106
The NOV introduces a deep occupied state in the a-SiO2 band gap, similar to
that of α-quartz. This state sits on average 0.81 eV above the a-SiO2 valence band
and ranges over 1.02 eV. An unoccupied state associated with the defect breaks off
from the bottom of the a-SiO2 conduction band, sitting on average 0.92 eV below
it and ranging across 2.83 eV. The distribution of these states is plotted against the
Si–Si vacancy bond length. The defect’s occupied states show a weak correlation with
the bond length, while the unoccupied states show a weak anticorrelation. Mulliken
charge analysis reveals that the Si atoms surrounding the vacancy average at +1.05
Page 81
|e| and range from +1.0 to +1.10 |e|, much more negative than the +1.43 |e| of bulk Si
in a-SiO2 . It is interesting to note that the charges of the two Si atoms surrounding
the vacancies display a random asymmetry, that is there is no predictable way of
determining whether one Si will be more or less negatively charged than the other Si
making up the vacancy.
4.3.2 Positive Oxygen Vacancy
The positively charged oxygen vacancy was studied in 20 different a-SiO2 models.
Starting from the NOV configurations described in the previous section, a hole was
added to each of the 20 systems and their geometries were optimised. In all con-
figurations, this resulted in a significant displacement of the Si atoms surrounding
the vacancy so that the average distance between them is 3.16 Å, ranging between
2.7 and 4.3 Å away from each other. The Si–O bonds of these Si atoms averages at
1.58 Å, slightly shorter than the average Si–O bond length of 1.61 Å in a-SiO2 . The
atomic structure of these configurations are analogous to the dimer configuration of
the positive oxygen vacancy in α-quartz. The unpaired spin is localised between these
two Si atoms, with their Mulliken spin moments averaging at 0.43 and ranging from
0.34 to 0.59. Their average Mulliken charge is +1.49 |e| and ranges from +1.47 to
+1.51 |e|. The narrow range of the Mulliken spin moments and charges indicates that
both Si atoms surrounding the vacancy are very similar.
The positive oxygen vacancy possesses an occupied one-electron state in the α-
channel of electrons that sits on average 1.0 eV above the valence band. However, it
shows a wide distribution ranging over 2 eV with a strong dependence on the Si–Si
distance of the two Si atoms surrounding the vacancy, shown in Fig. 4.9. The hole
introduced into the β-spin channel of electrons localises on a state which sits deep in
the a-SiO2 band gap, on average 3.8 eV above the valence band, and its position also
shows a strong dependence on the Si–Si distance.
Page 82
Figure 4.9: One-electron level of the positive oxygen vacancy in a-SiO2 plotted against
the Si–Si length of the Si atoms surrounding the vacancy. The left panel shows the
states from the α spin channel while the left panel shows the β spin states. The red
circles are the occupied and the green circles are the unoccupied states. The average
of the bands and their edges are coloured in orange, with the valence band at the
bottom and the conduction band at the top.
4.3.3 Singly Negative Oxygen Vacancy
To model the negatively charged vacancy, a single electron was added to each of the
20 neutral oxygen vacancy configurations described in section 4.3.1. This resulted
in a strong local distortion around the vacancy, with the Si–Si bond extending to
an average of 3.1 Å, similar to the positive oxygen vacancy. The range of bond
lengths, however, is much narrower for the negatively charged vacancy, ranging over
0.3 Å. The average Mulliken charge of the Si atoms surrounding the vacancy drops
to an average of +0.98 |e|, ranging from +0.8 to +1.1 |e|, suggesting that these Si
sites are more negatively charged than the average Si in bulk a-SiO2 . The unpaired
spin is highly localised between these two Si atoms, with their average Mulliken spin
Page 83
moment averaging at 0.46 and ranging from 0.24 to 0.72. The broad distribution of
Mulliken spin moments suggests that the electron has a tendency to localise at one
of the Si sites around the vacancy. However, there is no clear correlation to suggest
at which Si site the electron would localise. The calculated one-electron levels for the
Figure 4.10: One-electron level of the negative oxygen vacancy in a-SiO2 plotted
against the Si–Si length of the Si atoms surrounding the vacancy. The left panel
shows the states from the α spin channel while the left panel shows the β spin states.
The red circles are the occupied and the green circles are the unoccupied states. The
average of the bands and their edges are coloured in orange, with the valence band at
the bottom and the conduction band at the top.
20 configurations of the negative oxygen vacancy are plotted against their respective
Si–Si bond lengths in Fig. 4.10. Unlike the positive oxygen vacancy, there is no strong
correlation of the defect’s one-electron level with the Si–Si bond length. The addition
of an electron into the α-channel of electrons localises an occupied state which sits
deep in the a-SiO2 band gap, averaging at 4.1 eV above the valence band and shown
on the left panel of Fig. 4.10. Visualisation of this occupied state indicates that the
Page 84
electron now occupies an anti-bonding orbital between the Si atoms surrounding the
vacancy, explaining the extension of the Si–Si bond relative to the neutral oxygen
vacancy.
4.3.4 Summary
The oxygen vacancy has been calculated in 20 models of a-SiO2 . The structures of
these vacancies are analogous to their counterparts in crystalline α-quartz. However,
due to the disorder inherent to the amorphous structures, the defect’s structural
and electronic properties can be highly dependent on their respective local chemical
environments.
4.4 Discussion
The oxygen vacancy has been modelled in both α-quartz and a-SiO2 in a number
of charge states. The first thing to note is that the vacancies modelled here in α-
quartz show analogous structures in a-SiO2 . Although this is true of these models,
this result should not be extrapolated to the conclusion that all defects in α-quartz
have analogues in a-SiO2 . Due to the disorder inherent to the amorphous structures,
the defect’s structural and electronic properties can be highly dependent on their
respective local chemical environments, unlike their counterparts in α-quartz.
Having calculated the vacancy in its various charge states, we can examine their
thermodynamics. The formation energy of the vacancy from bulk SiO2 was calculated
according to equation B.1 in the appendix. Fig. 4.11a shows the formation energies
of the vacancy in α-quartz plotted against the Fermi Level of the system, with zero
set as the top of the α-quartz valence band. We find that the neutral configuration is
thermodynamically stable over a limited range of Fermi energies, between ≈ 4.0 and
5.5 eV. The formation energies of the vacancy in a-SiO2 have also been calculated.
Instead of plotting the formation energies of each vacancy, it can be instructive to
inspect their thermodynamic switching levels; the Fermi level at which the formation
Page 85
(a) Formation energy plotted against the Fermi level for the oxygen vacancy in α-quartz.
The straight black line is the neutral charge state, the red and green lines are the dimer and
puckered configurations of the positively charged vacancy, respectively, while the blue line
is the negative charge state.
(b) Histogram of the thermodynamic switching levels of the oxygen vacancy in a-SiO2 plotted
against the Fermi level. The (+/0) switching level is coloured in green while the (0/-)
switching level is coloured in blue.
Figure 4.11: Thermodynamics of the oxygen vacancy in α-quartz and a-SiO2 .
Page 86
energy of a defect in two different charge states is equal, as indicated by the circles in
Fig. 4.11a. These are plotted in a histogram for all 20 configurations of the a-SiO2
vacancies studied and is shown in Fig. 4.11b. It is clear to see that the switching levels
are distributed over a range of Fermi levels, in contrast to the discrete switching levels
of the vacancy in α-quartz. The (+/0) switching levels occur at a much lower Fermi
level in a-SiO2 than in quartz, but the (0/-) level occurs at around the same Fermi
level. These switching levels indicate that the vacancy can be thermodynamically
stable in the neutral charge state over a wide range of Fermi levels ranging from 0.6
to 6.1 eV, although it is not necessarily the same defect configuration that is stable
over this range. In contrast, the vacancy in α-quartz is thermodynamically stable in
the neutral charge state over a relatively limited range of Fermi levels, from 4.0 to
5.6 eV. This would make the neutral oxygen vacancy dominant in an Si/SiO2 system
(c.f. Si/SiO2 band offset 107 ), regardless of whether the SiO2 is in the crystalline or
amorphous state. However, when the Fermi level is set to the top of the valence band,
the positive charge state is dominant by far while the negative charge state is most
stable when the Fermi level is set to the bottom of the conduction band.
Page 87
Electron Trapping at Impurities in Silicon
5
Dioxide
5.1 Introduction
Having reviewed the well known oxygen vacancy defects and their charge trapping
states in SiO2 , we now turn our attention to electron trapping at impurities. Silicon
dioxide (SiO2 ) is a material that is of fundamental importance to a number of tech-
nologies. Quartz, the most stable crystalline polymorph of SiO2 , is widely utilized in
optical devices due to its high melting temperature and transparency to visible and
UV radiation. Electron trapping is known to have a dramatic effect on the perfor-
mance and reliability of these devices. For example, quartz, which is normally clear
Chapter 5. Electron Trapping at Impurities in Silicon Dioxide
and transparent to visible radiation, becomes coloured upon irradiation at higher

energies. 108 The colouration is thought to be due to electrons trapping at localized
impurity sites throughout the structure. These trapped electrons absorb throughout
the electromagnetic spectrum and will therefore affect any radiation passing through
the sample. These interactions can attenuate the signals passing through the structure
and are therefore usually undesirable properties.
So far, the dominant electron traps in quartz have been associated with impurity-
related centres. It has been well established, both experimentally and theoreti-
cally, that electrons can be trapped by Ge impurities substituting for Si in both
α-quartz 16,109 and in a-SiO2 , 85 with experimental models of these centres recently
revisited by Griscom. 110 Pacchioni et al. 16 calculated a germanium substitution in
a silica cluster (with fixed hydrogen atoms terminating the cluster) representative
of α-quartz; the cluster used was a H12 O16 Si4 Ge cluster. DFT was employed with
the exchange and correlation described by the hybrid B3LYP functional. 35,111 The
orbitals on Ge and the first O and Si neighbours were described by a 6-31 G basis
set while the rest of the Si and O atoms were described by a 3-21G basis set. The
cluster was terminated by hydrogens described by a MINI-1 basis set. They showed
that electrons can indeed be trapped at Ge impurities in SiO2 in a deep, localized
state. The electron trapping is accompanied by a strong structural relaxation.
Ge dopants in Si-based technologies have come back into the fore due to a number
of technological advances. The use of SiGe alloy materials has allowed the discovery
of new, advanced electronic devices and has improved the performance of existing
Si devices. 112 The structure of a SiGe device is formed by implantation of Ge into
Si followed by oxidation. The native oxide formed is mainly SiO2 due to the highly
negative free energy of formation of SiO2 , which is much more negative than that
of GeO2 , rejecting the Ge and forming a sharp Ge rich layer at the interface. 113
Nonetheless, several studies have investigated the concentration of Ge in the native
oxide under various processing conditions. 114–116 It is found that Ge persists in the
oxide, and that the Ge present in the oxide exists as GeO2 , not as elemental Ge
Page 89
clusters which were originally thought to exist due to the processing conditions. 117
Therefore it is important for us to understand the charge trapping properties of Ge
impurities in SiO2 .
In addition, a trapped electron has been observed in α-quartz containing Li impu-
rities using ESR. 118 It is formed after a double irradiation, with the second irradiation
thought to trap an electron at the free Li ion. Analysis of the experimental conditions
and the ESR parameters obtained indicates that the electron is in fact trapped at a
four-coordinated Si site while the Li ion resides nearby and stabilizes the trapped elec-
tron. Although this defect has been identified experimentally, there is no theoretical
model of this defect centre reported in the literature.
In this chapter, we investigate electron traps at Ge and Li impurity sites in α-
quartz. An atomistic defect model is developed for each impurity and they are con-
trasted against each other. The structure of these electron traps is found to be
qualitatively very similar. Finally, each model is compared to experimental data.
5.2 Calculation Details
All calculations discussed in this chapter were performed in a 3 × 3 × 3 supercell of

α-quartz whose structure was obtained from a crystal database. 119 This supercell was
then optimized in the neutral, spin singlet charge state using DFT implemented in
the CP2K code, 27 as discussed in chapter 2. The non-local functional PBE0 TC LRC
was used in all calculations with a cutoff radius of 2.0 Å for the truncated Coulomb
operator. 34 Inclusion of Hartree-Fock exchange provides a more accurate calculation
of the band gap compared to the semi-local functionals such as PBE. It also better de-
scribes localized states that may be involved in charge trapping processes such as those
studied in this chapter. The CP2K code uses a Gaussian basis set with an auxiliary
plane-wave basis set. 25 The Gaussian basis set employed for all atoms was a double-ζ
basis set with polarization functions 91 in conjunction with the Goedecker-Teter-Hutter
(GTH) pseudopotential. 90 Calculating hyperfine interactions necessitated the use of
Page 90
all electron basis sets using the Gaussian and augmented plane-wave (GAPW) ap-
proach. The basis sets with contraction schemes of (8831/831/1),(8411/411/11) and
6-311G** were used for silicon, 120 oxygen 121 and Li, 122 respectively. The plane wave
cut-off was set to 5440 eV (400 Ry). In order to reduce the computational cost of
calculating the Hartree-Fock exchange terms, an auxiliary basis set was used as de-
scribed in chapter 2.3.2. This auxiliary basis set is a much sparser Gaussian basis
set containing less diffuse and fewer primitive Gaussian functions than the full basis
set employed in the rest of the calculation. This allows the Hartree-Fock exchange
terms, whose computational expense scales to the fourth power of the number of basis
functions, to be calculated on a much smaller basis set than the rest of the calculation
and therefore much faster. The pFIT3 auxiliary basis set was used for the Li, Si and
O atoms, with the Hartree-Fock exchange terms calculated using a full basis set on
the Ge atoms.
All geometry optimizations were performed using the Broyden-Fletcher-Goldfarb-
Shanno (BFGS) optimizer to minimize forces on atoms to within 37 pN (2.3 ×10−2
eV Å−1 ). After optimization in the neutral charge state, a single Si atom in the cell
was substituted for a Ge atom and optimized again to create the Ge-doped α-quartz
model. To create the Li-doped quartz models, no Si atoms were substituted but an
Li interstitial was instead added to the system. One extra electron was then added to
the Ge-doped system along with a background positive charge that is uniform across
the entire cell while the cell containing the Li interstitial was charge neutral. The
optimization resulted in the electron trapping centres in these systems.
In order to check the convergence of each defects’ properties, the electron trapping
centres were also modelled in α-quartz cells containing 576 atoms. However, the
results reported in this chapter are for the cells containing 243 atoms. The results for
the larger cells have not been included as they are consistent with those obtained in
the smaller cells, indicating that the 243 atom cells of α-quartz are sufficient for this
study. Barriers between configurations were calculated using CI-NEB (see chapter
2.6). Linear interpolation was used to generate 10 images which made up the band
Page 91
to be optimized with each of the images connected by a spring with a force constant
of 2 eV Å2 .
The calculated energies were corrected, where necessary, for the interaction be-
tween the charged defects using the method of Lany and Zunger, 123,124 described in
equation B.2 in the appendix.
5.3 Results of calculations
5.3.1 Ge-doped α-quartz
The first centre that was studied was the electron trapping at Ge impurities in α-
quartz. 109 The model of this centre has been developed both experimentally and the-
oretically, with it first being characterized by Weil et al. using ESR, 109 and reviewed
recently by Griscom. 110 The experimental characterization of electrons trapped in
Ge-doped α-quartz reveals two distinct defects: the Ge(I) and Ge(II) centres. 109,125
It has been suggested that both these defects reside on the same GeO4 tetrahedron
in α-quartz and a-SiO2 , with the Ge(I) centre assigned as the ground state. 109,110
The energy difference between these two defect centres was measured as 0.0025 eV
at 15 K and 0.0078 eV at 300 K. 109 On the theoretical side, cluster calculations by
Pacchioni et al. of what he called the Ge electron centre (GEC) have demonstrated
that a four-coordinated Ge atom in silica can trap an electron. This is accompanied
by an orthorhombic distortion of the pseudo-tetrahedral Ge centre which results in
two short and two long Ge–O bonds. 16
Here, we use the same terminology as Pacchioni and calculate the GEC in α-quartz.
Calculations were performed in 243 atom periodic cells of α-quartz with one Si atom
substituted for Ge. Geometry optimization in the neutral charge state maintains the
tetrahedral coordination of the Ge impurity, similar to the remaining Si atoms in bulk
α-quartz. However, the Ge–O bonds are slightly longer than its Si counterpart, with
two long (1.74 Å) and two short (1.73 Å) Ge–O bonds and O–Ge–O angles of ≈ 109◦ .
The Ge introduces an unoccupied state which sits below the bottom of the α-quartz
Page 92
Figure 5.1: a) The molecular orbital associated with the lowest unoccupied state in
Ge-doped α-quartz. It is strongly localized on the Ge impurity and its neighbours.
The isovalue of the density surfaces is 0.02. b) The density of states of neutral Ge-
doped quartz is shown underneath the structure. An unoccupied state which is high
in Ge and O character sits below the bottom of the conduction band.
conduction band and is highly localized on the Ge atom and its neighbours, as can
be seen in Fig. 5.1a. The calculated one-electron band gap of α-quartz is 8.6 eV.
The empty state sits 0.8 eV below the bottom of the conduction band and projection
of the density of states onto the individual atoms further confirms that the state is
highly localized on Ge and O atoms (see Fig. 5.1b). The neutral Ge atom has a
Mulliken charge of +1.76 |e| in α-quartz.
An extra electron was then added to this cell to make the GEC. It initially occupied
the localized, empty state on the Ge atom. The Coulombic interaction between the
extra electron and its slightly negatively charged O neighbours repels them away in
order to lower the total energy of the system. This results in an opening of the
O–Ge–O angle formed by the longer Ge–O bonds which in turn localizes the electron
further on the Ge atom. The two long Ge–O bonds elongate by 0.2 Å, becoming
1.9 Å, and the O–Ge–O angle opened up to 150◦ . The remaining two short Ge–O
bonds extended slightly, up to 1.8 Å. The bond lengths and angles of the negatively
charged GEC centre are asymmetric, similar to the Ge impurity in the neutral charge
state. The spin density is strongly localized on the Ge atom and its oxygen neighbours
(see Fig. 5.2a), with the Ge atom possessing a Mulliken spin moment of 0.72. The
Mulliken charge of the Ge atom is now +1.33 |e|, about 0.43 |e| lower than in the
Page 93
Figure 5.2: a) The spin density of the negatively charged Ge-doped α-quartz with
the spin density localized on the Ge centre and its neighbours. The isovalue of the
density surfaces is 0.02. b) The density of states of the negatively charged Ge centre
in quartz. An occupied state sits deep in the band gap and is high in Ge character.
neutral system. The optimization results in the occupied one-electron defect level
sitting 4.1 eV below the conduction band, as can be seen in Fig. 5.2b. The relaxation
energy of the system (i.e. the total energy difference between the system with the
electron at the bottom of the conduction band and the final structure) is 1.51 eV.
As mentioned above, it has been suggested that the Ge(I) and Ge(II) electron
spin resonance (ESR) signals can result from the electron trapping on the same GeO4
tetrahedron. 110 Calculations were made to check whether an electron could localize
in different configurations on the same GeO4 tetrahedron by creating initial configu-
rations that would favour these metastable states. In particular, the five remaining
O–Ge–O angles were opened on the same GeO4 tetrahedron described above to see
if an electron can localize on each O–Ge–O different angle in the tetrahedron. This
was accomplished by manually displacing the two neighbouring O ions of the origi-
nal, ideal GeO4 tetrahedron to create an electron trapping precursor. The two Ge–O
bonds associated with the displaced O atoms were 1.9 Å and the angle between them
was somewhere between 150◦ and 160◦ . The geometry of these systems, each con-
taining a newly opened O–Ge–O angle, were then optimized in the negative charge
state. Energy minima associated with all 5 remaining combinations of Ge–O bond
pairs were found. Their relaxation energies range between 1.36 eV to 1.51 eV. The
Page 94
Mulliken charges of all the Ge atoms are lower than in the neutral system, ranging
from +1.34 to +1.38 |e|. These results suggest that an electron can indeed localize
on the same GeO4 tetrahedron Ge doped α-quartz in six different configurations.
All six configurations show qualitatively similar geometries: elongation of two
Ge–O bonds and the opening of O–Ge–O angle between them. However, there are
quantitative differences in the extent of O–Ge–O angle opening which ranges from
132◦ to 150◦ . The electronic density of states of all six configurations shows an
occupied one-electron level that sits 4.31 eV on average below the bottom of the α-
quartz conduction band with a distribution of 0.1 eV. The relative energies of the six
different configurations are plotted with respect to the O–Ge–O angle in Fig. 5.3. A
general trend can be seen; the lowest energy configuration has the widest O–Ge–O
angle and smaller O–Ge–O angles result in higher total energies. The smaller O–Ge–O
angles are due to the different local environments that each O–Ge–O angle exists in.
The α-quartz unit cell is not highly symmetric with each O–Ge–O angle being in a
unique chemical environment. Due to this, some of the O–Ge–O angles have larger
space to relax into.
The isotropic hyperfine splitting constants were calculated for all six configura-
tions of the Ge electron centre. The splitting due to the interaction between the
Ge73 nucleus ranges from -26.325 mT to -31.19 mT, with the lowest energy structure
possessing an isotropic hyperfine constant of -29.17 mT. The calculated values are in
the range of the experimental values of -28.47 mT and -28.956 mT reported by Isoya
et al. for the Ge(I) and Ge(II) centres. 109 Although it was found that an electron
can localize on the same GeO4 tetrahedron, the structure of the Ge(II) centre cannot
be correlated with any of the structures. More accurate calculations of the hyperfine
splittings would be needed in order to conclusively identify the atomic structure of
the Ge(II) centre.
Page 95
Figure 5.3: Energies of the six different electron trapping configurations in Ge-doped
α-quartz plotted alongside the O–Ge–O angle. The energies are shown on the left
hand side y-axis and marked as black lines on the graph, while the O–Ge–O angles
are shown on the right hand side y-axis and marked on the graph as red crosses.
5.3.2 Li-doped α-quartz
Jani et al. experimentally studied the effect of Li impurities in α-quartz, 118 particu-
larly properties of the [AlO4 /Li]0 centres. A Li electron centre in quartz was found to
be formed by a two-step irradiation process. The first irradiation step performed at
150-300 K moves the Li away from its Al counterpart. 126 The second step, performed
after cooling the quartz sample down to 77 K, forms a [SiO4 /Li]0 centre. The ESR
spectrum of this centre shows splittings of 0.09 mT from a 7 Li and 40.47 mT from a
29
Si. This centre was found to be stable below 180 K and has been characterized by
Jani et al. as an extra electron trapped at a four-coordinated Si site with an adja-
cent Li+ ion providing stability. 118 This model has been supported by early molecular
cluster calculations by Wilson et al., 127 but no further calculations have been carried
out to confirm the structure of this defect.
In this chapter, we attempt to further elucidate the structure of the electron trap
in Li-doped α-quartz. In these calculations, a Li atom was added to a 3x3x3 supercell
of α-quartz and the geometry of the system was optimized in the neutral charge state
Page 96
Figure 5.4: a) Spin density of a Li interstitial atom in α-quartz. The spin density
is highly localized on the Li atom and its surrounding Si neighbours. The isovalue
used to plot this surface was set to 0.02. b) The density of states associated with
the interstitial Li in α-quartz. An occupied state sits just below the bottom of the
α-quartz conduction band and is highly localized on the Li atom as well as the Si and
O atoms.
initially. The Li atom occupies an interstitial position in the α-quartz lattice, as can
be seen in Fig. 5.4a, with a one-electron level sitting ≈ 1.5 eV below the bottom of the
α-quartz conduction band shown in Fig. 5.4b. It sits in an α-quartz channel between
four O atoms, the nearest neighbours sitting 1.96 or 1.98 Å away. The Li remains
rather neutral with a Mulliken charge of +0.2 |e|. The spin density is highly localized
on the Li centre with a Mulliken spin moment of 0.75. These results indicate that
the interstitial Li remains atomic in nature with no electron transfer taking place.
The average Mulliken charge of the Si atoms in this model is +1.44 |e|, ranging from
+1.40 to +1.49 |e|.
Inspired by the opening of the O–Ge–O angle in Ge-doped α-quartz, perturbation
of the lattice in a similar manner to the Ge electron trap was investigated to see if
it would induce electron transfer from the Li atom to nearby O–Si–O angles. An
O–Si–O angle of one the nearest Li interstitial’s neighbours was opened in a manner
similar to the O–Ge–O described above. Relaxing the structure results in an electron
transferring from the Li atom and localizing on the perturbed O–Si–O angle, as shown
Page 97
in Fig. 5.5a. The relaxed structure has an extended O–Si–O angle of 150◦ while the
Li ion is located 2.62 Å away from the Si centre, bound to the two O neighbours,
which are not associated with the O–Si–O angle opening. The spin density plot of the
system in Fig. 5.5a shows that the unpaired spin is mostly localized on the Si atom
(in the open O–Si–O angle) and two of its oxygen neighbours. Mulliken population
analysis reveals that the Li now has a charge of +0.49 |e| while the Si at the centre
of the wide O–Si–O angle has a Mulliken charge of +1.1 |e|. The Mulliken charge of
Si ions in quartz is +1.43 |e|, indicating that the Si has gained negative charge. The
Mulliken spin moment of the Li atom is now reduced to 0.04 while the spin moment
on the Si is 0.86.
The Li+ ion is stabilized by its interaction with the lone pairs on its oxygen
neighbours. The total energy of the Li stabilized electron centre is 0.28 eV lower
than that of the Li interstitial atom in α-quartz, i.e. the trapping energy of the Li
stabilized electron centre is 0.28 eV. The small trapping energy reflects the fact that
the initial state is already a fairly deep electron trap.
The barrier for transferring an electron from the atomic Li to the Si ion was
calculated using the CI-NEB method, described in chapter 2.6. Starting with the
interstitial as the initial configuration and the electron trapping state as the final
configuration, a band was set up by interpolating between these two points and was
subsequently optimized. The barrier was found to be 0.68 eV. De-trapping from this
state requires overcoming a barrier of 0.96 eV. The occupied one electron state of the
unpaired electron is located 3.1 eV below the bottom of the quartz conduction band
(see Fig. 5.5b). This demonstrates the stabilizing role of the Li+ ion in creating the
[SiO4 /Li]0 electron centre.
The calculated hyperfine splittings due to the interaction of the unpaired electron
with the surrounding nuclei were calculated and are shown in Table 5.1. They are also
compared with the experimental results by Jani et al. 118 in that table. The strongest
hyperfine interaction is with the 29 Si ion, with an isotropic hyperfine splitting of 43.07
mT. Smaller hyperfine splittings come from the 7 Li ion and on the 17
O neighbours.
Page 98
Figure 5.5: a) Spin density of the electron trapped at an Si atom from an Li impurity.
The spin is highly localized on the Si atom and its O neighbours. Isovalue of the
density plot is 0.015. b) Density of states of the electron trap. An occupied state sits
deep in the band gap and is high in Si and O character.
Signal Theory Experiment 118

aiso Si 43.07 40.47
aiso Li 0.11 0.09
Principal values 0.089 0.088
0.096 0.098
0.16 0.15
Table 5.1: Hyperfine splittings and principle values of the hyperfine tensor (in mT)
of the Li electron trap in α-quartz. The experimental values of hyperfine interactions
for the Li-doped quartz are shown for comparison.
These results are in excellent agreement with experiment, as can be seen in table 5.1.
To summarize, an Li impurity in α-quartz can lead to an electron being trapped
at an Si atom in a similar manner to the Ge electron trap in α-quartz. A common
feature of both centres is that the electron is highly localised on either Ge or Si ion
and is accompanied by an energy gain, elongation of two metal-oxide bonds and a
significant opening of the –O–(Ge)Si–O– angle.
Page 99
5.4 Discussion and Conclusions
These calculations demonstrate the qualitatively similar character of electron local-

ization for two impurities in crystalline SiO2 . In α-quartz, a substitutional Ge atom
provides a local structural perturbation which facilitates the localization of an extra
electron at the Ge site. For the first time, the atomistic structure of the [SiO4 /Li]0
centre has been calculated, showing that a Li atom in α-quartz donates an electron
to a neighbouring Si ion and further stabilizes the defect state by the Coulomb in-
teraction between the trapped electron and the Li+ ion. In both cases, the electron
localization on Ge and Si ions is facilitated by the opening of the O–Si(Ge)–O angle.
The electron localization in pure bulk α-quartz requires opening O–Si–O angle from
109◦ to 134◦ , but introducing this distortion costs 0.57 eV.
It is interesting to note that the experimentally measured ESR signal of the
[SiO4 /Li]0 centre in Li-doped α-quartz is reduced to zero at around 180 K. 118 The
calculated barrier for de-trapping is rather high at 0.96 eV, suggesting that it is per-
haps not the electron transferring back to the Li atom which is responsible for the
disappearance of the ESR signal. As mentioned earlier in this section, the interaction
between the trapped electron and the Li ion provides a stabilizing Coulomb potential
well for the trapped electron. As all the Si ions in α-quartz are equivalent, diffusion of
the Li ion into a nearby equivalent position provides an equally stabilising Coulomb
potential well for another equivalent Si atom. This will lead to electron tunnelling
from the original Si site to a new site. This may offer a possible alternative explana-
tion to the temperature dependence of the ESR signal. The calculated barrier for Li+
ion diffusion in pure α-quartz is 0.4 eV. 128 Using NEB, the barrier for diffusion of Li+
of the [SiO4 /Li]0 centre between equivalent sites across a ring in α-quartz was calcu-
lated as 0.56 eV. This Li+ ion displacement is accompanied by an electron transfer
to another Si site and is equivalent to diffusion of the whole centre. At low temper-
atures the Li ions cannot overcome this diffusion barrier, hence the strong, discrete
ESR signal. However, at increased temperatures Li ions start moving around and the
ESR signal should decrease until it vanishes completely due to rapid electron transfer
Page 100
between Si sites.
To summarize, these results demonstrate that the character of electron trapping
for the Ge and Li impurities in α-quartz are qualitatively similar. The structure of
both electron traps involves an opening of an O–Si(Ge)–O angle. In chapter 6, we
discuss whether electrons can be trapped at intrinsic sites in SiO2 in a similar manner
as we have seen in this chapter.
Page 101
Intrinsic Electron Trapping Sites in Silicon
6
Dioxide
6.1 Introduction
As mentioned in the introduction of chapter 5, electron trapping is known to have

a dramatic effect on the performance and reliability of electronic devices employing
SiO2 as a gate insulator and charge trap flash memory devices. 84,129 In that chapter,
we focussed on electron trapping at impurities in SiO2 and saw that electrons could
be trapped in Ge- and Li- doped SiO2 in deep states in the band gap. We now turn
to whether electrons can be trapped at intrinsic sites, i.e. unaided by impurities or
coordination defects.
Chapter 6. Intrinsic Electron Trapping Sites in Silicon Dioxide
Little is known regarding the possibility of intrinsic electron trapping in the a-

SiO2 network. Bersuker et al. used molecular models to suggest that electrons can
be trapped by Si–O bonds in a-SiO2 , leading to their weakening and thus facilitating
bond dissociation. 130 Using a F3 Si–O–SiF3 cluster, chosen to simulate the structure of
two SiO4 tetrahedra facing each other, they showed that an extra electron introduced
into this cluster may localize on a Si–O bond, causing the other Si–O bond associated
with the oxygen to contract and strengthen. The extra electron localizes in an oxygen
‘p’ state, significantly weakening one Si–O bond. The O–Si–O angle after the electron
has been localized on the Si–O bond is 144◦ .
Further calculations by Camellone et al. have shown that electrons can be trapped
in non-defective continuum random network models of a-SiO2 . 131 In their study, sev-
eral structures of 72 atoms of a-SiO2 were generated using classical MD simulations.
The electronic structures of these models were then calculated using DFT, utilizing
the generalized gradient approximation (GGA) for the exchange-correlation functional
and including a self-interaction correction. The potential energy surface of the sys-
tem with an extra electron was explored along a reaction coordinate defined as the
elongation of one Si–O bond. The global energy minimum corresponds to the neu-
tral equilibrium geometry where the electron is delocalized over the system, but a
metastable state was found where the Si–O bond extended to 1.83 Å. In this state,
an electron is localized on the Si atom. Extension of the Si–O bond and electron
localization also resulted in expansion of the O–Si–O angle up to 156.44◦ . The barrier
from the delocalized state to the metastable localized state was found to be 0.23 eV,
with the localized state higher in energy by 0.17 eV. Recent calculations have also
demonstrated that silicon dangling bonds at SiO2 surfaces are deep electron traps and
can form corresponding negatively charged defects. 132 However, these theoretical pre-
dictions have not yet been confirmed experimentally due to challenges in identifying
and characterizing electron trapping centres in SiO2 .
Unlike in optical fibres and other optical devices, where electrons and holes are
created by electronic excitation, in MOS devices they are often injected from the Si
Page 103
substrate. For example, electron trapping at an energy of 2.8 eV below the conduction
band of a-SiO2 has been observed using photon-stimulated tunnelling experiments in
device-grade oxides grown on Si and SiC crystals in a series of papers. 133–136 Further-
more, low-temperature capacitance 137 and Hall effect measurements 138,139 on 4H-SiC
MOS devices revealed that the density of these electron trapping states can be as
high as 1014 cm−2 eV−1 . The trap density of 1013 cm−2 measured inside a 2-nm thick
near-interface SiO2 layer 133,136 corresponds to ≈ 5 × 1019 cm−3 in terms of volume
concentration. This is much higher than observed densities of the established intrinsic
defects in thermally grown a-SiO2 . The absence of a comparable density of electron
traps in bulk a-SiO2 suggests that electron trapping at 2.8 eV deep centres takes
place not on pre-existing defects but rather at intrinsic sites in the oxide network it-
self. Whether the substrate plays any role in stabilizing these traps remains unclear.
These results, as well as the insight gained from the electron trapping at impurities
described in chapter 5, motivate further investigation of the possibility of electron
trapping in amorphous silica network.
In this chapter, electrons are shown to trap in continuous non-defective a-SiO2
networks, forming deep electron states in the gap. The geometric structure of these
centres is highly similar to that of electrons trapped by Ge, where the key to the
electron trapping is the wide opening of the O–Ge–O angle, or Li centres in quartz,
where it is facilitated by the opening of the O–Si–O angle. It turns out that precursor
Si sites with wide enough O–Si–O angles, which occur naturally in a-SiO2 structures,
can facilitate spontaneous electron trapping at Si sites. Using this fingerprint we
estimate the concentration of intrinsic electron trapping sites in a-SiO2 .
6.2.1 Classical Calculations
The calculations presented in this work make use of both classical force-fields and ab
initio theory. The ReaxFF 51 force-field was used to generate 20 models of amorphous
Page 104
SiO2 , each containing 216 atoms, modelled within periodic boundary conditions and
as described in chapter 2. ReaxFF was parametrized to reproduce the properties of
various silica polymorphs, small silica clusters and silicon polymorphs. 52 This force-
field allows one to calculate Si and O atoms in varying oxidation states based on the
instantaneous geometry of the system. This is accomplished by assigning a charge
dependent atomic energy and exploiting the electronegativity equalization principle. 42
The extended bulk silica structures used in this study - containing up to 401,760
atoms - were generated using the potential created by van Beest, Kramer, and van
Santen (BKS). 38 It takes the functional form:
qi qj −6
V (r) = + Aij e(−bij r) − Cij rij (6.1)
r
The parameters for Si and O are taken from Ref. 38 and are included in table A.1
in the appendix. This Buckingham-type potential allows one to perform calculations
much faster than the ReaxFF potential and is more suited to creating large a-SiO2
structures. As can be seen in the results, comparing results obtained with two very
different force-fields gives more confidence in our predictions. All classical atomistic
simulations were performed using the LAMMPS code. 61 To generate amorphous struc-
tures, classical molecular dynamics simulations were run using ReaxFF and BKS to
melt and quench crystalline SiO2 structures into an amorphous state, as described
in chapter 3. However, for the BKS simulations, the temperature was raised to 7000
K instead of the 5000 K melting temperature used for ReaxFF. The differences in
melting temperature stem from the differences in parametrization and use of different
functional forms in the potentials. The BKS structures have a higher density of 2.37
g cm−3 . The Si–O bond lengths of the BKS structures average at 1.61 Å, while the
O–Si–O angles average at 108◦ and the Si–O–Si angles average at 142 ◦ .
Page 105
Figure 6.1: A schematic of one-dimensional diabatic potential energy surfaces corre-

sponding to an initial electronic state and two final electronic states of the system
with extra electron. The energy labelled EB is the thermal barrier to electron trap-
ping. The energies labelled ET are the trapping energies, calculated as the total
energy difference between the initial and final state. There are two trapping energies
shown in the figure. The physical meaning of the negative E1T is that the final state
is thermodynamically unstable with respect to the initial state. The final state 2 is
thermodynamically stable with respect to the initial state.
6.2.2 Density Functional Theory Calculations
Density functional theory (DFT), implemented in the CP2K code, was used to further
optimize geometries of the ReaxFF structures and calculate their electronic struc-
tures, 27 as described in the calculation details section of chapter 5.
To discuss the electron trapping by impurities in perfect crystalline or amorphous
SiO2 structures, the total energies of the initial and final states with the extra electron
in the system were compared, as illustrated in Fig. 6.1. The left diabatic curve
(labelled as ‘initial state’ in Fig. 6.1) represents the system with an extra electron
in the initial state while the right curves (labelled as ‘final state 1 or 2’) represent
the final state with the electron localized on, for example, a Ge impurity or on a Si
ion in the pure matrix after full geometry relaxation. EB is the thermal barrier for
electron trapping from the initial state to the localized state and ET is the trapping
energy, calculated as the total energy difference between the initial and final electronic
states. In this chapter, two different scenarios are discussed. Trapping from the initial
state to the ‘final state 1’ in Fig. 6.1 requires a thermal barrier to be overcome.
Page 106
The final state is higher in energy, corresponding to a negative trapping energy, and
is thermodynamically less favourable than the initial electronic state. The second
scenario corresponds to electron trapping from the initial state to the ‘final state 2’ in
Fig. 6.1. This electron trapping is barrier-less or has a small barrier, the final state is
thermodynamically more favourable and corresponds to a positive ET . The physical
meaning of the initial and final states is discussed below for each particular system.
6.3 Results of Calculations
6.3.1 Electron Trapping in Bulk α-Quartz
An extra electron was added into the perfect α-quartz lattice and the geometry was
optimized using DFT, resulting in the electron becoming fully delocalized at the
bottom of the conduction band with no change in the structure of the lattice. This is
in good agreement with previous theoretical calculations. Using the insights gained
from electron trapping at impurities in SiO2 , perturbations to the α-quartz structure
were investigated to see if they could lead to localization of an electron. The key to
electron trapping in Ge- and Li-doped α-quartz seemed to be the wide O–(Ge)Si–O
angle. Therefore, the O–Si–O angle was perturbed in order to see whether this would
allow an electron to localize at that site.
The SiO4 tetrahedra in α-quartz are made up of two shorter Si–O bonds and two
longer Si–O bonds and a discrete O–Si–O angle which measures ≈ 109.5◦ . These
calculations started from the neutral optimized structure of α-quartz. A random Si
atom was chosen and the two O atoms associated with the two longer Si–O bonds
were displaced so that an O–Si–O angle became greater than 135◦ . The O–Si–O
angle associated with the two longer Si–O bonds was chosen by analogy with the
lowest energy electron trap in Ge-doped α-quartz. An extra electron was then added
to the now perturbed system and the geometry was optimized. This lead to the extra
electron localizing on the Si atom at the centre of the Si site around which the O atoms
had been displaced. The geometry optimization further opens the O–Si–O angle to
Page 107
161◦ and the two Si–O bonds elongate from 1.61 to 1.74 Å, while the other two bonds
of the tetrahedron both elongate to 1.69 Å. This structural relaxation is similar to the
one observed for both the Ge and Li electron centres in chapter 5. The spin density
of the system is now highly localized on the Si site, as can be seen in Fig. 6.2a. The
relaxation is highly localized around the defect, and as can be seen in Fig. 6.2a the
surrounding α-quartz lattice remains highly crystalline. Mulliken population analysis
reveals that the charge of the Si ion, on which the electron is localized, is +1.04 |e|,
significantly lower than the +1.43 |e| average charge of Si in α-quartz, and that the
Mulliken spin moment is 0.86, indicating that the spin is highly localized on it. The
electron trap creates a one-electron state 2.5 eV below the bottom of the α-quartz
conduction band which is principally Si and O ‘sp’ in character. The position of this
level can be seen in the density of states plotted in Fig. 6.2b.
The barrier to self-trapping an electron into this state from a delocalized conduc-
tion band state was calculated using CI-NEB as 0.57 eV, while the (self)-trapping
energy of this system was found to be −0.32 eV. This indicates that the self-trapped
electron polaron state in pure α-quartz is thermodynamically unstable with respect
to the delocalized state (see Fig. 6.1). De-trapping from the localized state into the
delocalized state requires overcoming a barrier of 0.25 eV. The O–Si–O angle at the
maximum of this barrier is 134◦ .
To summarize, these results show that an electron can be trapped in deep states at
Si sites in α-quartz. The trapping is driven by a wide O–Si–O angle. In bulk α-quartz,
all bonds have an equilibrium value of 1.61 Å, therefore to create an electron trapping
site, an O–Si–O angle has to be opened. This requires overcoming an energetic barrier
of 0.25 eV.
6.3.2 Electron trapping in amorphous SiO2
Electron trapping in a-SiO2 was studied using twenty periodic models of bulk a-
SiO2 . Each model contained 216 atoms. The geometries of the ReaxFF generated
amorphous structures were optimized in the neutral charge state within DFT and then
Page 108
Figure 6.2: The square modulus of the wave function of an extra electron occupying
the lowest state at the bottom of the conduction band of a-SiO2 . The bigger spheres
connected to four atoms are Si atoms and the smaller spheres connected to two atoms
are O atoms. The darker blobs represent the magnitude of the modulus of the wave
function. The isovalue used to represent the square modulus of the wave function is
0.005.
an extra electron was added to each model. In chapter 3, the geometrical distributions
of a-SiO2 were discussed with Fig. 3.10 showing the distributions of Si–O bond
lengths, O–Si–O and Si–O–Si angles obtained after DFT geometry optimization of the
neutral cells. The geometrical properties of the optimized structures change slightly
with respect to those obtained with ReaxFF, whose distributions can be seen in Fig.
3.7. The Si–O bond lengths of these 20 models after DFT optimization average at
1.62 Å, ranging from 1.58 Å to 1.67 Å. The Si–O–Si angles average at 147◦ , ranging
from 112◦ to 179◦ , while the O–Si–O angles average at 109◦ , ranging from 95◦ to
137◦ . Analysis of the ring size distribution after DFT optimization shows the 4- and
5-member rings to be dominant with smaller contributions from 3- and 6- member
rings. The electronic structure calculations predict an average one electron band
gap of 8.1 eV, ranging from 7.7 to 8.3 eV over all 20 models. For comparison, the
one-electron band gap of α-quartz is 8.6 eV.
These calculations start from 20 models of a-SiO2 optimized in the neutral charge
state. In each model, an extra electron was added so that each system is negatively
Page 109
charged. Along with the extra electron, an extra uniform positive charge was added
to each system. The extra electron initially occupies the lowest unoccupied state. In
all structures, this state is partially localized on several Si and O ions, as illustrated
in Fig. 6.3a for one of the structures. In comparison, a state which sits deep in the
conduction band is shown in Fig. 6.3b. This band-like state is delocalized over all
atoms in the system, in contrast to the partially localized lowest unoccupied state.
The energetic position of the lowest band-like state is highlighted as the mobility
edge on the density of states plotted in Fig. 6.3c. Each system’s geometry was then
optimized with the extra electron, resulting in spontaneous electron localization in
four out of the twenty models. This localization was accompanied by a strong local
distortion around a single SiO4 tetrahedron, similar to the trapped electron centre
relaxation in α-quartz. In each of the four structures, the extra electron is localized
on one Si ion, with the two O neighbours repelled so that an O–Si–O angle is opened
to ≈ 172◦ , as shown in Fig. 6.4a. The Si–O bonds making up this O–Si–O angle
extend from 1.63 and 1.64 Å to 1.78 and 1.82 Å, respectively (see Fig. 6.4). Mulliken
population analysis shows that as a result of the electron localization, the Si ion
charge decreases by about 0.25 |e| compared to its charge in the neutral system. The
average gain in energy resulting from the barrier-less electron localization (ET in Fig.
6.1) in the four models is 1.25 eV, ranging from 0.72 to 1.71 eV. The electron state
occupied by the extra electron is located at ≈ 3.17 ± 0.05 eV below the bottom of
the SiO2 conduction band, indicating a deep electron trap. This was calculated as
the difference in energy between the highest occupied and lowest unoccupied single-
particle states. The single-particle state of the electron trap from the four models is
plotted in Fig. 6.4.
The hyperfine splittings induced by the localized electron have been calculated and
are shown in Table 6.1. The strongest hyperfine interaction is with the Si ion; however,
there is a significant interaction with the nearby oxygen atoms. Interestingly, some
of the hyperfine interaction values are similar to the experimental values of a defect
known as the E0 centre in amorphous silica. 101 This is not surprising considering the
Page 110
Figure 6.3: The square modulus of the wave function of an extra electron occupying
the lowest state at the bottom of the conduction band of a-SiO2 . The bigger spheres
connected to four atoms are Si atoms and the smaller spheres connected to two atoms
are O atoms. The darker blobs represent the magnitude of the modulus of the wave
function. The isovalue used to represent the square modulus of the wave function is
0.0005.
Page 111
strong electron localization on one Si ion. However, more models would be needed to
give a more accurate distribution of the hyperfine values accessible to experiment.
In all four cases, the Si ion on which the electron traps initially forms the widest
O–Si–O angle in the sample, exceeding 132◦ . In the sixteen remaining a-SiO2 samples,
where the distribution of O–Si–O angles was slightly narrower, the extra electron
remained delocalized. To investigate whether this wide O–Si–O angle held the key
to the electron localization, structural distortions were introduced to make two other
random O–Si–O angles the widest in two separate systems. An angle in one of the
systems was increased from 120.3◦ to 132.1◦ . Adding an extra electron into this
system and optimizing the geometry results in the electron localizing on the Si ion
within the perturbed angle and causes it to open further to 160.68◦ . An angle in a
separate system was changed from 121.3◦ to 132.0◦ . When the electron was added
to this system, the O–Si–O angle opened to 164.5◦ . These results demonstrate that,
although a wide O–Si–O bond angle serves as an efficient precursor to spontaneous
electron trapping in a-SiO2 , thermally activated trapping can also take place at other
sites by opening an O–Si–O angle. These results also make apparent the link between
the geometric structure of a trap and its electronic properties and allows one to use
the criterion of wide O–Si–O angles as a fingerprint for identifying precursor sites
for spontaneous electron localization in initial a-SiO2 structures and estimating the
concentration of such sites.
6.3.3 Concentration of electron trapping sites in a-SiO2
As suggested above, by analysing the structure of an a-SiO2 sample for the presence of
O–Si–O angles exceeding 132◦ one can estimate the lower limit of the concentration of
precursor sites which can act as electron trapping centres. The results from the twenty
models of a-SiO2 samples indicate that the presence of an O–Si–O angle exceeding
132◦ always leads to spontaneous localization of extra electrons in a-SiO2 . This angle
is at the tail of the O–Si–O angle distribution in regular SiO2 structures constructed
using the ReaxFF potential and optimized using DFT.
Page 112
Atom Bond length/Å Values / mT

Si -50.98
-45.45
-45.23
O 1.82 -4.181
-2.660
-2.624
O 1.78 -5.714
-4.357
-4.298
O 1.70 -1.548
-1.216
-1.212
O 1.70 -1.581
-1.264
-1.259
Table 6.1: Geometrical parameters and the average principal values of the hyperfine
tensor of the electron trap in a-SiO2 from the four models. The bond lengths shown
are with respect to the Si atom on which the electron is trapped.
To test whether the existence of these precursor sites and their concentration de-
pends on the model of amorphous structure and to obtain better statistics, three
additional samples of amorphous SiO2 were created using the BKS interatomic po-
tentials 38 as described in section 5.2. These potentials are often used in studying the
structural properties of a-SiO2 and give good agreement with experimental data. 59,68,69
The three samples created have dimensions of 50 × 25 × 5 nm3 , 25 × 12.5 × 2.5 nm3 ,
and 12.5 × 7 × 1.5 nm3 and include 401,760; 55,296 and 8,640 atoms, respectively.
These models were searched for O–Si–O angles exceeding the fingerprint value of 132◦
to estimate the concentration of electron trapping precursor sites. The concentration
of such sites in all models proved to be very similar with a volume concentration of
4 × 1019 cm−3 in the sample containing 401,760 atoms. It is interesting to note that
in spite of the difference in cell sizes and force-fields used, this concentration agrees
well with the original observation of four trapping sites in twenty 216 atom samples.
Page 113
Figure 6.4: Atomic structure and spin density distribution of an intrinsic electron
trap in a-SiO2 . We highlight the SiO4 tetrahedron on which the electron traps and
show the spin density only on the nearest ions. The spin density isovalue is 0.02.
These calculations demonstrate that electrons can indeed trap at intrinsic sites in both
crystalline and amorphous SiO2 . The distribution of geometrical parameters in the
disordered a-SiO2 leads to the existence of precursor Si sites which can spontaneously
trap an electron in a state ≈ 3.2 eV below the bottom of the conduction band. The
estimated concentration of these precursor sites is 4 × 1019 cm−3 . The large average
distance between precursor sites suggests that diffusion of trapped electrons via a
thermally activated tunnelling mechanism should be quite inefficient and they are
more likely to move via thermal activation into the mobility edge states of amorphous
silica at high temperature.
These results differ from those previously reported by Bersuker 130 and Camel-
lone 131 which focussed on the effect of the Si–O bond length and its relation to intrinsic
electron trapping in SiO2 . The calculations presented above indicate that the O–Si–O
angle is a more efficient precursor for electron trapping in SiO2 . The differences in
these results from those presented by Camellone et al. could stem from our use of a
non-local functional as opposed to the generalized gradient approximation (GGA), 131
as GGA tends to underestimate the degree of electron or hole localization. 140,141
Page 114
The high volume concentration of large O–Si–O angle electron trapping precur-
sor sites suggests that the electron trapping can be abundant in a-SiO2 samples.
However, identifying these electron traps in relatively pure bulk samples may require
irradiating at liquid nitrogen temperatures, where both trapped electrons and holes
are immobile. 142 These results suggest that trapped electrons can be stable even at
room temperature. However, trapped holes become mobile below 200 K and can re-
combine with electrons. These results also support the common perception that the
abundance of impurities, such as Al, Ge, Li, Na and water in quartz as well as in silica
glass samples may lead to efficient electron trapping by impurity centres and further
hamper the identification of intrinsic electron traps in a silica network.
These intrinsic trapping states are correlated to electron traps identified exper-
imentally in MOS devices 143 at an energy of 2.8 eV below the conduction band of
a-SiO2 grown on Si and SiC crystals. 133–136 Experimentally, these traps were popu-
lated by illuminating the MOS structures by photons of energy sufficient to excite
electrons from the semiconductor valence or conduction band above the edge of the
SiO2 conduction band. They had initially been correlated with oxygen deficient cen-
tres at the near-interfacial oxide. 134,135,144 However, later experiments on nitridated
SiC/SiO2 samples questioned this attribution, particularly when taking into account
the fact that the density of known O-deficiency centres (E0γ and E0δ centres) rarely
approaches the density range of 1013 cm−2 found for the 2.8 eV deep electron traps.
Although these electron traps are especially pronounced in 4H-SiC/SiO2 devices, they
seem to play a role in all devices containing SiO2 as the insulating material, suggest-
ing that they may be intrinsic to the oxide. For instance, these traps are expected
to appear below the conduction band of Si nano-crystals in the case of quantum
confinement. 145,146
The intrinsic electron traps in a-SiO2 discussed in this chapter could be good
candidates for understanding these data. The calculated concentration of precur-
sor sites for intrinsic electron traps approaches the experimentally observed value for
the states filled by photo-stimulated tunnelling from SiO2 valence band. However,
Page 115
populating such a density of electron traps via electron injection from an electrode
through the SiO2 conduction band should be much less efficient because an electron
capture event requires dissipating about 1.5 eV of relaxation energy into phonons
during the trapping process. This is likely to be slower than fast electron transport in
the conduction band of a thin oxide towards an opposite electrode. In order to keep
the additional (unpaired) electrons on these centres, one must ensure sufficiently high
strength of electric field is externally applied to the interface. The latter can hardly be
realized under the conditions of an ESR experiment because the presence of conduct-
ing electrodes impairs the quality factor of the microwave resonator. Furthermore,
all available experimental evidence concerns interfaces of SiO2 with semiconducting
materials (Si, SiC), which might suggest that the experimentally observed high prob-
ability of electron trap occupation may be related to the strain in the SiO2 network
near the interface, while in the bulk of the film their concentration may be lower.
The observed trap photo-ionization energy at 2.8 eV is between the values calculated
for α-quartz (2.5 eV) and a-SiO2 (3.0 eV). The results in this chapter indicate that
the geometry of the oxide structure can significantly affect the position of the defect
level, and the discrepancy between the experimental value and our a-SiO2 value may
reflect the higher oxide density in thermally grown oxides 2,147 rather than the density
obtained in this work.
To summarize, these results demonstrate that, similar to holes, 142 electrons can
be trapped at structural precursor sites in an amorphous silica matrix, forming deep
electron states in the oxide band gap. The geometric structure of trapped electron
centres in a-SiO2 are qualitatively similar to the extrinsic electron trapping centres
studied in chapter 5. In a-SiO2 , these states may be responsible for the electron
trapping observed at interfaces of SiO2 -based MOS devices and should be present in
bulk SiO2 samples.
Page 116
Optical Absorption Spectra of Trapped
7
Electrons in SiO2
7.1 Introduction
The previous two chapters presented calculations which characterized intrinsic and
extrinsic electron traps in a-SiO2 . However, unambiguously identifying all sites re-
sponsible for electron trapping in a-SiO2 has proved particularly challenging because
of a large number of possible charge redistribution channels and the presence of water
and impurities in most samples. Thus, spectroscopic evidence of intrinsic electron
trapping in a-SiO2 is still missing. In this chapter we calculate the optical absorption
spectra of the intrinsic electron traps in a-SiO2 using the embedded cluster method
Chapter 7. Optical Absorption Spectra of Trapped Electrons in SiO2
and TDDFT described in chapter 2 and compare it to low temperature experimental

optical absorption spectra measured by Prof. Katsumi Tanimura.
In order to assess the accuracy of the method, we first calculate the optical ab-
sorption spectrum of the GEC in α-quartz, described in chapter 5, as an absorption
band at 4.2 eV has been universally attributed to it. We then calculate the spectra
for the intrinsic electron traps in a-SiO2 , and by analysing this data and comparing
it to experimental optical absorption spectra, an absorption band was identified that
can indeed serve as a signature of intrinsic electron traps.
7.1.1 Experimental Observations
In order to identify absorption bands associated with electron trapping centres in

SiO2 , optical absorption spectra of electron irradiated SiO2 were measured at 80 K by
Prof. Katsumi Tanimura of Osaka University. A detailed account of the experimental
observations is given below.
Three different, high-purity a-SiO2 samples were used in order to differentiate the
results that are specific to only a given sample from those that are common to all
samples. These include Suprasil W1 (SW1), Viosil-SCF (SCF), and Viosil-SZ (SZ);
all fabricated by Shin-Etsu Quartz, LTD. SW1 has less than a few ppm of OH ions
but includes a considerable amount of dissolved O2 . SZ contains almost no chemical
impurities - concentrations of both Cl and OH ions are less than 1 ppm - however,
it includes 500 ppm of oxygen deficiency (so called ODC(I)) centres that exhibit a
strong absorption peak at 7.6 eV. On the other hand, SCF has no Cl ions (less than
1 ppm), but includes about 1000 ppm of OH ions. All samples were cut and polished
to a size of 8 × 10 × 1 mm and irradiated with electron pulses (1 MeV, 20 ns duration)
generated from a Febetron (HP43710A) at 80 K. The optical absorption spectra were
measured using a spectrophotometer (Shimazu UV-3100). A Xe-arc lamp with optical
filters was used as the light source for optical excitation of the irradiated specimens.
Red light with a wavelength range from 600 to 780 nm, and UV light with a peak at
320 nm (3.87 eV) and a band pass of 10 nm was generated by combining appropriate
Page 118
Figure 7.1: a) Optical absorption spectra of electron-irradiated a-SiO2 (SZ). The black
curve was recorded at 80 K after electron irradiation. The red curve was recorded at
80 K after an anneal at 186 K. The blue curve is the difference between the black and
red curves, revealing the absorption components due to centres annihilated during
the anneal. b) and c) show the difference spectra between the absorption bands of
electron-irradiated a-SiO2 and those after selective excitation: b) by red light (600 -
780 nm) and c) by UV light (320 ± 10 nm). The difference spectra are de-convoluted
into the absorption bands due to STH2 (thin blue solid curve at 2.2 eV), a 3.7 eV
band (thick solid curve), a 4.7 eV band (thin red curve), the E0 band (broken curve
centred at 5.8 eV) and the rest of the spectrum (solid green curve).
filters.
The black curve in Fig. 7.1a shows the optical absorption spectrum induced by
electron-pulse irradiation at 80 K. Several peaks are clearly distinguishable. The
lowest-energy peak at 2.2 eV is the absorption band due to self-trapped holes with two-
centre type configurations (STH2). 148 The peaks at 5.1 eV and 5.7 eV are associated
with B2 and E0 centres, respectively. The strong absorption peaking at >6.5 eV
is usually attributed to ODC(II) or B2 centres. 149 We note that similar peaks with
almost the same heights have been generated in SW1 and SCF samples, which contain
148
almost no oxygen vacancies (see Fig. 2 in ref. ). The peak at 3.7 eV is an, as yet,
unidentified band which is commonly generated in a-SiO2 samples irradiated and
Page 119
measured at low temperatures and thus not related to specific pre-existing lattice
imperfections in a-SiO2 .
To reveal the origin of the 3.7 eV band, its growth kinetics was measured as a
function of the electron dose. In Fig. 7.2, the growth characteristics of the E0 centres
(5.7 eV), STH2 (2.2 eV), and the 3.7 eV absorption band are shown. The yield of E0
centres increases linearly with the dose, while those of the STH2 and the unidentified
centre are strongly correlated and show a tendency to saturate at doses greater than
≈ 1.5 × 105 Gy.
The thermal stability of the defects responsible for the optical absorption spectrum
in Fig. 7.1a was studied by a pulse annealing method. The red line in Fig. 7.1a shows
the absorption spectrum after warming the irradiated sample up to 186 K for 5 mins,
followed by absorption measurement at 80 K. The peaks at 2.2 eV and 3.7 eV are
completely annealed out, which is also demonstrated by the blue curve showing the
difference between the initial low temperature and the post-anneal spectra. It is also
evident that absorption bands with energies above 4.5 eV are also partially annihilated
after annealing at 186 K.
The results obtained during the growth and anneal processes indicate that the
3.7 eV band is either an absorption band of STH2 or is due to a centre formed as a
partner of a self-trapped hole. The former possibility can be ruled out based on the
Figure 7.2: Growth kinetics of the optical absorption bands of electron irradiated
a-SiO2 measured at 80 K. The 2.21 eV and 3.70 eV peaks show a trend to saturate
at a dose greater than 1.5 ×105 Gy.
Page 120
symmetry of the centre revealed by the polarized bleaching experiments. 148 The 2.2
eV band shows a clear dichroism by which the asymmetric spectral shape of the STH2
band could be determined. However, the 3.7 eV band does not show any dichroism,
indicating that the centre is not a two-centre type.
Alternatively, the 3.7 eV band could originate from an electron component of
electron-hole pairs created by the irradiation. Part of the holes self-trap, hence the
STH2 signal observed in our optical absorption spectra (see Fig. 7.1). Another part
could be trapped by oxygen vacancies forming E0 centres, which absorb at 5.7 eV, or
by other pre-existing defects. Electrons can also trap on pre-existing oxygen vacan-
cies, 150 impurities, or on intrinsic network sites. 151 The calculated optical absorption
of negatively charged oxygen vacancies peaks at about 3.3 eV 150 and they can con-
tribute to the 3.7 eV absorption in SZ samples which include considerable amount of
oxygen vacancies. However, as noted above, the 3.7 eV band is induced also in SW1
and SCF samples, which do not include pre-existing oxygen vacancies.
Optical bleaching experiments at 80 K can shed more light on the origins of ab-
sorption spectra by releasing trapped holes and electrons into valence and conduction
bands, respectively, as a result of selective excitation. Fig. 7.1b exhibits a difference
spectrum of the sample before and after bleaching by red light (600-780 nm), which
selectively excites only STH2. This excitation results in bleaching of the 3.7 eV band,
as well as other higher energy bands. Interestingly, the B2 band is little affected by
the excitation, although a substantial amount of the E0 band is bleached. Similarly,
Fig. 7.1c shows the spectral components after selective excitation by UV light (320
± 10 nm). The excitation in the 3.7 eV band results in bleaching of the 2.2 eV band
and some of the higher energy bands including that of the E0 centre.
The ratio between the 2.2 eV and 3.7 eV bands is not the same in the two different
photo-bleaching events, demonstrating that they are associated with different centres
and support the assertion that the centre giving rise to the 3.7 eV absorption band is
an electron partner associated with STH2. Figs. 7.1b and c also demonstrate that the
decay of 3.7 eV band is always associated with the decay of the spectral components
Page 121
around 4.7 eV and above 6 eV. This feature is not specific to the SZ sample, but is
common to SCF and SW1 and hence characteristic to localized electron centres in
a-SiO2 .
To summarize, optical absorption measurements made after electron irradiation
at low temperature exhibit an unidentified peak at 3.7 eV in addition to the peaks at
2.2, 4.7 and 5.7 eV associated with defects. The 3.7 eV peak is correlated with STHs
in a-SiO2 but is distinct from that centre. The remainder of the chapter is devoted
to calculating the optical absorption spectra of intrinsic electron traps in a-SiO2 and
comparing them to this experimental data.
To investigate the origins of the absorption bands of the intrinsic electron trap in
a-SiO2 , their optical transitions were calculated using an embedded cluster method
described in chapter 2.5 implemented in the GUESS code. 88 The inhomogeneous
broadening of the absorption spectra due to the structural disorder of the glass samples
was included by using seven independent models of a-SiO2 which were generated as
described in chapter 3. To assess the reliability of the method, the optical absorption
spectrum of the GEC in α-quartz was calculated initially and compared to its well
known experimental spectrum. To model the GEC, only one model of the α-quartz
lattice was used.
Each model was adapted for use in an embedded cluster model as follows. A
spherical nanocluster was generated from the unit cells obtained in chapter 3. Each
nanocluster was divided into two regions - I and II. Region I is located in the centre
of the nanocluster and is further subdivided into three concentric regions of different
physical descriptions. The centre of Region I consists of a cluster of atoms which are
treated using DFT implemented in the Gaussian09 code. We used a hybrid BxLYP
functional with x=32.5% exchange, which provides a better description of the band
gap than the B3LYP functional 35 and improves the functional’s linearity. 152 The Si
Page 122
and Ge atoms are described with a 6-31G* basis set while the O atoms with a 6-31G
basis set. The outer section of Region I is described by a core-shell variant of the
classical BKS potential. 87 In the core-shell model, an atom is divided into two point
charges connected by a spring, allowing the atoms to change their polarization. The
central quantum region is interfaced with the outer classical ions by the embedding
pseudopotentials, described in detail in ref. 87 The outer Region II is composed of
non-polarizable ions which are fixed in perfect lattice sites and interact with the ions
in Region I by the classical BKS potential without shells. 38 All atoms in Region I
are allowed to relax during geometry optimizations to provide the atomic structure
of the electron traps. The optical absorption spectra of the electron traps were then
calculated using the time-dependent DFT (TDDFT) method, 24 as implemented in
the Gaussian09 code.
7.3 Results of Calculations
7.3.1 Optical Absorption of the Ge Electron Trapping Centre
The optical absorption spectrum of the well-known Ge electron centre (GEC) in α-

quartz was calculated. The GEC was used to ensure that the embedding methods
can accurately calculate the optical absorption spectra of localized electron centres in
SiO2 , due to their similar structures and due to the fact that the GEC’s absorption
spectrum has been measured. The GEC in α-quartz is responsible for the Ge(I) ESR
signal. 109 This signal strongly correlates with an optical absorption band at 4.2 eV,
which is universally attributed to this centre. 125,148,153,154 Optical absorption spectra
of Ge-doped SiO2 include additional bands that peak at ≈ 5.8 eV and higher energies,
which have also been tentatively assigned to the GEC in a-SiO2 . 155
A quantum cluster containing 72 atoms and representing the α-quartz structure
was embedded into the rest of the crystal as described in the calculation details. A
single Si atom at the centre of the quantum cluster was substituted for a Ge atom and
the geometry of the system was optimized in the neutral charge state. The relaxation
Page 123
Figure 7.3: a) Atomic structure and spin density of the GEC in α-quartz. The GeO4
tetrahedron is highlighted, where the pink ball is Ge, the Si atoms are yellow balls,
the O atoms are red balls and the spin density is the blue polyhedron. b) Simulated
optical absorption spectrum of the Ge electron trap in α quartz. The black line
is a sum of all Gaussian broadened excitation energies weighted by their respective
oscillation strengths. The bar plot shows the excitation energies and their respective
oscillation strengths. Red bars correspond to electronic transitions in the α channel
while blue bars correspond to β channel.
resulted in an extension of the Ge–O bonds to 1.72 Å compared to the Si–O bonds
(1.64 Å) in α-quartz; however, the GeO4 tetrahedron remains intact. The calculated
band gap of the system is 7.2 eV with a Ge unoccupied state located just below the
conduction band and strongly localized on the Ge atom. These results are consistent
with the calculations on the GEC in a periodic model described in chapter 5.3.1.
An extra electron was then added to the system followed by a full geometry op-
timization of Region I. This resulted in a strong local relaxation, where the O–Ge–O
angle opened from 109◦ up to 148◦ and the extra electron localized on the Ge atom,
as shown in Fig. 7.3a. The Ge–O bond lengths extend asymmetrically, so that two
bonds measure 1.83 Å while the other two bonds measure 1.87 Å. This defect intro-
duces an occupied electron state which sits 5.0 eV below the α quartz conduction
band.
The optical transition energies and oscillator strengths of the GEC in α-quartz
were then calculated. The calculated optical absorption spectrum is plotted in Fig.
7.3b as a solid black line. The bar plot shows all one-electron transitions weighted by
Page 124
their oscillator strengths. Each absorption is Gaussian broadened by 0.3 eV in order

to simulate homogeneous broadening. The total spectrum shown by the black line is a
sum of broadened excitation energies, weighted by their respective oscillator strengths.
The optical transitions can be categorised into two types. The transitions of the first
type are highlighted as red bars while those of the second type are highlighted as
striped blue bars. The first type is caused by transitions of the unpaired electron
localized on the wide O–Ge–O angle into quasi-local states at the bottom of the a-
SiO2 conduction band (see also transition 1 in the inset of Fig. 7.4b). The calculated
spectrum exhibits a peak with a maximum at 4.2 eV which is a transition of this first
type, in excellent agreement with the experimental signal. The second type is due to
transitions of a β-spin electron from the occupied non-bonding oxygen ‘p’ orbitals,
which have broken away from the top of the a-SiO2 valence band due to the structural
distortion introduced by the defect, to the unoccupied state which is localized on the
wide O–Ge–O angle (see also transition 2 in the inset of Fig. 7.4b). In our calculated
spectrum, they comprise 2 peaks with maxima at around 5.8 and 7.1 eV. The O ‘p’
states, which split from the top of the valence band, are localized around the wide
O–Ge–O angle, similar to the O ‘p’ states of the nearest O neighbours highlighted in
Fig. 7.3a. The 5.8 eV peak agrees very well with the band recently assigned to the
GEC in a-SiO2 . 155
7.3.2 Optical Absorption of the Electron Trapping Centre
Intrinsic electron trapping sites were selected based on a geometrical criterion of an

O–Si–O angle greater than 132◦ . As has been shown in chapter 6, extra electrons can
trap at such precursor sites without barrier. Each quantum region was constructed
around such a wide O–Si–O angle, with the number of atoms ranging from 96 to
114. Due to the disordered nature of a-SiO2 and the rather stringent criterion of
O–Si–O angles being greater than 132◦ , it was not possible to standardize the numbers
of atoms described quantum mechanically across the seven different glass samples.
The quantum region was then interfaced with a classically polarizable region which
Page 125
contained between 910 and 933 atoms. The diameter of the entire system was ≈ 50
nm while the quantum regions are ≈ 1.5 nm in diameter. The positions of all atoms
in Region I were optimized in the neutral charge state. The band gap averaged over
the seven different samples in the embedded cluster model is 7.7 eV, ranging from 7.4
to 8.1 eV.
An extra electron was added into each of the seven samples and their respective
geometries were optimized, resulting in seven different intrinsic electron trap con-
figurations, where an electron has localized on an Si atom and an O–Si–O angle has
opened from 132◦ to 174◦ , ranging over 5◦ . These electron traps introduce an occupied
single-particle state located on average 4.3 eV below the bottom of the conduction
band in the a-SiO2 matrix. The single-particle wave function of this state is shown in
Fig. 7.4a.
The transition energies and oscillator strengths of all seven electron trapping con-
figurations were calculated using TDDFT and plotted together with equal weights in
Fig. 7.4b. The bar plot shows all one-electron transitions corresponding to the two
types explained in the inset of Fig. 7.4b with their respective oscillator strengths.
Each absorption line was then Gaussian broadened by 0.3 eV, the same as the GEC,
to simulate homogeneous broadening. The black line in Fig. 7.4b was produced as a
sum of all Gaussian broadened optical transitions in seven electron trapping centres
weighted by their respective oscillation strengths. The strongest transitions of the first
type are highlighted as red bars while the second type of transitions are highlighted
as striped blue bars. The first type is caused by transitions of the electron localized
on the wide O–Si–O angle (see Fig. 7.4a) into quasi-local states at the bottom of the
a-SiO2 conduction band composed of ‘d’ orbitals of nearby Si atoms and nearby O ‘s’
orbitals. It exhibits a peak with a maximum at 3.7 eV. The second type of excitations
is due to transitions of a β-spin electron from the occupied non-bonding oxygen ‘p’
orbitals, which have broken away from the top of the a-SiO2 valence band due to the
structural distortion introduced by the intrinsic electron trap, to the unoccupied state
in the band gap localized on the wide O–Si–O angle. They comprise a peak with a
Page 126
maximum at around 6.4 eV. The O ‘p’ states, which split from the top of the valence
band, are localized around the wide O–Si–O angle, similar to the O ‘p’ states of the
nearest O neighbours highlighted in Fig. 7.4a. The whole spectrum is much broader
than that for the GEC described in the previous section because transitions in each of
the seven configurations have different energies due to disorder in a-SiO2 and different
local environments.
Figure 7.4: a) The one-electron wave-function of the electron localized in a-SiO2 .

The tetrahedron with the wide O–Si–O angle is highlighted. b)The calculated optical
absorption spectrum of intrinsic electron traps in a-SiO2 . The inset shows a schematic
of the optical transitions that can occur in the intrinsic electron trap. Two types of
transitions are shown by two colours. See text for discussion.
We note that the calculated spectrum has transitions in the energy range of 3-7
eV. This is consistent with both the anneal and spectral bleaching data, which show
significant reduction of optical density at energies exceeding 4 eV (see Fig. 7.1).
The nature of both types of transitions is the same as those of the GEC in α-quartz
described in the previous section.
Calculated optical absorption spectra of intrinsic electron traps in a-SiO2 exhibit a

band centred at 3.7 eV. This compares very well with experimental optical absorption
spectra of electron-irradiated a-SiO2 recorded at 80 K, where an unidentified peak has
been measured at 3.7 eV and where self-trapped holes are immobile. These results
Page 127
provide the strongest evidence yet for the formation of intrinsic localized electron
centres in a-SiO2 as a counterpart to self-trapped holes and a detailed model of a
strongly localized electron in a disordered oxide. This localization is promoted by
strained O–Si bonds and occurs on distorted SiO4 tetrahedra. Therefore one can
expect that similar electron traps can exist in other types of silica glasses, 156 as well
as in feldspar and other silicates based on networks of SiO4 tetrahedra. Many of
these materials contain, as yet unexplained, deep electron centres responsible for
therm-luminescence employed in e.g. luminescence dating. 157 Therefore, the optical
signatures of intrinsic electron traps revealed here will be useful for understanding the
properties of other silicates as well as the properties and stability of electronic and
optical devices which utilize a-SiO2 .
Page 128
Hydrogen’s Interactions with Silicon
8
Dioxide
8.1 Introduction
Hydrogen plays an important role in a number of technological applications which

utilise a-SiO2 . In particular, hydrogen has long been used in both the microelec-
tronics and fiber optics industries to passivate defects which exhibit electronically or
optically detrimental properties (such as the Pb and E 0 centres 158 ). Despite its use in
passivating defects, it is also known for its various roles in the chemical and physical
processes at Si/SiO2 interfaces and in bulk amorphous SiO2 . 86,159–162 Hydrogen has
been found to be involved in some hydrogen-complexed defects, such as the hydrogen
Chapter 8. Hydrogen’s Interactions with Silicon Dioxide
bridge, 11,163 the 74 G-, and the 10.4 G-doublet centre. 164 These defects are the subject
of intense research efforts across a range of sciences; they may lead to severe electronic
device degradation 165 and they play a role in radiation induced attenuation (RIA) in
optical fibers. 154 The formation mechanisms of the hydrogen-complexed defects are
largely unknown. In this chapter, we investigate hydrogen complexed defects which
could be relevant to silica’s technological applications and propose mechanisms for
their formation and further interactions with hydrogen.
The E0 centre is a well known defect in SiO2 . It was first characterized using ESR
techniques 101 and a theoretical model was proposed by Rudra and Fowler whereby
a hole is trapped at an oxygen vacancy. 166 However, recent experiments using field-
dependent recombination of holes trapped in SiO2 with injected electrons 167 revealed
that the paramagnetic state of the E0 centre is not always correlated with the entity
bearing the positive charge. It has been suggested that the positive charge is protonic
in origin, a hypothesis later corroborated by a number of experimental results. 167,168
Consequently the O3 Si H entity in a-SiO2 , which is easily identified due to its
strong infra-red absorbance signal, has been suggested as a possible E0 precursor, 167
where upon hole trapping hydrogen dissociates from the Si–H bond in the form of a
proton leaving behind a neutral paramagnetic E0 centre. The question remains as to
whether the liberated proton is then trapped in SiO2 or diffuses through and escapes
the SiO2 layer. A section of this chapter is devoted to understanding how a neutral
E0 centre can be formed from Si–H bonds in a-SiO2 .
In addition, EPR analysis of OH-rich synthetic silica after irradiation suggests
that atomic hydrogen, which can diffuse rapidly through a-SiO2 , 169 can interact with
the a-SiO2 network to create a charge neutral defect. 170,171 A 0.08 mT doublet due
to proton hyperfine splitting has been assigned to a Si dangling bond coordinated by
two bridging oxygens and an OH group. This centre is thought to result from the
interaction of H0 with electronically excited, strained Si–O bonds. However, reactions
of atomic hydrogen with strained Si–O bonds have not been investigated theoretically
and the perception that it interacts weakly with the a-SiO2 network prevails in the
Page 130
literature. 172–174 We investigate whether neutral, atomic hydrogen can interact more
strongly with the a-SiO2 network.
Previous theoretical studies have shown that hydrogen, both atomic and molecular,
interacts with point defects in SiO2 . Using a cluster model and post Hartree-Fock
(HF) methods , Edwards et al. studied the interaction of hydrogen molecules with
E0 -type defects (i.e. a O3 Si moiety) in α-quartz. 175 The clusters studied were
terminated with O–H groups with the terminating H atoms fixed in space. They
showed that the H2 molecule dissociates at an E0 centre to leave behind a free H atom
and an Si–H bond, passivating the 3-coordinated Si with a barrier of 0.78 eV. DFT
calculations using hybrid functionals performed on small clusters by Lopez et al. also
indicate the H2 cracking reaction at E0 centres to be endothermic with a barrier of
less than 0.5 eV. 176
In summary, hydrogen atoms and molecules are known, both theoretically and
experimentally, to interact with 3-coordinated Si defects in a-SiO2 and can passivate
the defect. In contrast, hydrogen has been shown to interact negligibly with the a-SiO2
matrix. 86,172–174 However, recent experiments suggest that atomic hydrogen can in fact
create new defects, although the formation mechanisms are entirely unknown. 171 In
addition, experiments also show that the 3-coordinated Si defects may have a charge
neutral analogue in addition to the well known positively charged vacancies known as
E0 centres in a-SiO2 . 167 Furthermore, it has been suggested that the charge state of
the E0 centre may be neutral and that it may be formed from an Si–H precursor.
In this chapter, atomic and molecular hydrogen interactions with bulk a-SiO2 were
investigated and shown to create two types of charge neutral defects. The first defect
involves hole trapping at Si–H bonds in bulk a-SiO2 and is shown to create neutral,
paramagnetic E0 centres. The presence of a proton nearby is shown to significantly
affect the defect’s one-electron levels in the band gap. The interactions of atomic
hydrogen with pure bulk a-SiO2 have also been investivated. We find that hydrogen
can bind to a bridging oxygen in bulk, stoichiometric a-SiO2 to generate an under-
coordinated Si defect centre facing a hydroxyl [O–H] group which shall be referred to
Page 131
as a hydroxyl E0 centre herein. This hydrogen-induced defect is charge neutral and

is thermodynamically stable over a wide range of Fermi levels including across the Si
band gap, suggesting that it could be detrimental to electronic device applications of
a-SiO2 . The origins of these defects and how they may play a role in electronic device
and optical fiber reliability issues are discussed.
The calculations presented in this chapter make use of both classical force-fields and
ab initio theory to understand how H interacts with a-SiO2 . The ReaxFF 51 force-field
was used to generate 116 models of amorphous SiO2 each containing 216 atoms and
periodic boundary conditions were employed, as described in chapter 3. In addition,
twenty models of a-SiO2 were generated with an extra hydrogen molecule introduced
into the melting phase in order to study the formation of defects from Si–H bonds.
Thus, these twenty models each contained 218 atoms. In order to ensure that Si–H
bonds are present in the final structure, the distance between a single Si and H atom
was fixed 1.46 Å apart and their coordinates frozen for the entire melt-and-quench
procedure. After quenching these structures, they were further optimized using DFT
as described in chapter 3. After the optimization, each model contained one Si–H
bond and one O–H bond and there were no coordination defects in any of the models.
An example of one of the models containing an Si–H and O–H bond is shown in Fig.
8.1.
Density functional theory (DFT), implemented in the CP2K code, was then used
to further optimize geometries of all the structures. 27 The same technical details were
used as are described in the calculation details of chapter 5.
Page 132
Figure 8.1: a) A model of a-SiO2 containing an Si–H bond and an Si–O–H bond along
with the HOMO of this system. The Si atoms are yellow spheres, O atoms are the
darker red spheres H atoms are pale grey spheres. The distance between the Si–H
bond and Si–O–H bond is ≈ 6 Å. b) Electronic density of states of a-SiO2 containing
an Si–H and O–H bond. They are projected onto Si, O and H atoms as well as
showing the total density of states. A state associated with the Si–H bond sits just
above the valence band.
8.3 Results of calculations
8.3.1 Atomic Hydrogen in a-SiO2
Atomic hydrogen’s interactions with a-SiO2 have been studied in the literature. 172–174
Almost all theoretical studies confirm that atomic H occupies an interstitial position
in SiO2 and prefers to be in the positive or negative charge state depending on the
system’s Fermi level. That is, interstitial atomic H is not thermodynamically stable
in the neutral charge state. In this section, we investigate whether neutral, atomic
hydrogen can interact more strongly with the a-SiO2 network using the DFT methods
described in chapter 2.3. Atomic hydrogen was found to interact strongly with long
Si–O bonds to produce a defect that is referred to as the hydroxyl E0 centre. The
interaction of further atomic H with this centre is also discussed.
Interstitial Hydrogen Atom
As a reference calculation for hydrogen’s interactions with a-SiO2 , we calculated the

structure and energetics of an interstitial neutral H atom in 26 different a-SiO2 ma-
trices in order to build up the statistics required for a disordered sample. A single
H atom was placed in a random position in each a-SiO2 sample under the constraint
Page 133
Figure 8.2: Atomic configuration and spin density of the hydroxyl E0 centre and its
precursor. The Si atoms are the bigger yellow balls, the O atoms are the smaller red
balls and the H atom is the small white ball in each configuration. The spin density
are the blue, transparent polyhedra. a) This is an unperturbed SiO4 tetrahedron. The
circled bridging O center has a statistically long Si–O bond which acts as a precursor
for the following hydrogen defect configurations. b) The hydroxyl E0 center. a 3-
coordinated Si where the spin density is localised on the Si and its three O neighbors.
The 3-coordinated Si faces a hydroxyl group.
that it was further than 2 Å away from any nearest neighbor. The total energy of
each system was then minimized with respect to its atomic coordinates. Consistent
with the literature, we find that interstitial H atoms do not bind and barely interact
with a-SiO2 matrices, with the nearest neighbors found at ≈ 2.6 Å. The spin density
of the system is almost entirely localized on the H atom with an average Mulliken
spin moment of 0.96 from the 26 systems. The electronic structures of the systems
show that the H atom introduces a one-electron level which sits in the a-SiO2 band
gap; 0.7 eV on average and ranging from 0.1 to 1.2 eV above the valence band. The
intersitial H atom can sit in any of the voids in a-SiO2 . As the structure is amor-
phous, the total energy of the system will be different depending on where the H
atoms sits. To establish the range of interstitial H atom energies, we placed H in 10
different interstitial positions in the same a-SiO2 structure and find that the range of
total energies is spread over 0.2 eV. These results are in accord with what has been
previously reported in the literature. 177
Page 134
Figure 8.3: Histograms of the bond length distributions around the hydroxyl E0 center.
a) Distribution of the O–H bonds of the hydroxyl group belonging to the hydroxyl E0
center. b) Distribution of the three Si–O bond lengths around the 3-cooordinated Si
of the hydroxyl E0 center. c) Distribution of the Si- -O distance where the Si belongs
to the 3-coordinated Si and the O belongs to the hydroxyl group [see Fig. 8.2]b. Note
that ‘- -’ indicates a non-bonding interaction in order to distinguish from the bonding
interactions represented by ‘–’.
Neutral Atomic Hydrogen’s Interactions with Bridging Oxygen Sites
To investigate whether a neutral H atom interacts more strongly with O sites in the
a-SiO2 network, we placed it at a distance of 1 Å from bridging O atoms in 116
different a-SiO2 samples [see Fig. 8.2a]. These geometries were optimized resulting in
a new type of defect described below.
Hydroxyl E0 Centre
When atomic hydrogen was placed 1 Å away from long Si–O bonds, shown in Fig.
8.2a, a defect was formed which resembles an E0 centre, i.e. a 3-coordinated Si atom
with an unpaired electron, 166 facing a hydroxyl group, shown in Fig. 8.2b. It is
therefore referred to as the hydroxyl E0 centre. Preliminary calculations whereby H
atoms were placed at random Si–O bonds showed that this defect configuration formed
predominantly at long (> 1.65 Å) Si–O bonds. The concentration of such bonds in
the a-SiO2 models discussed in chapter 3 is 2.2%. To obtain a statistical distribution
of the hydroxyl E0 centre’s properties, H atoms were placed 1 Å away from O atoms
Page 135
Figure 8.4: Histogram of the one-electron level of 116 configurations of the hydroxyl E0
centre. The energy scale starts from 0.0 and finishes at 8.1 eV; this is the a-SiO2 band
gap as calculated using the PBE0 TC LRC functional. The area of the histograms
which are coloured dark red show occupied states while the area of the histograms
which are coloured blue show the unoccupied states.
associated with long Si–O bonds in 116 different a-SiO2 matrices and their geometries
were optimized. In all the systems, this resulted in the H atoms breaking an Si–O
bond, forming an O3 Si centre facing a hydroxyl group as illustrated in Fig. 8.2b.
The O–H bond averages at 0.99 Å and shows a very narrow distribution, as can be
seen in Fig. 8.3a. The Si–O bonds of the O3 Si moiety average at 1.65 Å and
their distribution can be seen in Fig. 8.3b, while the distance between the Si and
the O from which it dissociated averages at 2.63 Å, a very wide distribution as can
be seen in Fig. 8.3c. The average Mulliken spin moment of the 3-coordinated Si is
0.90, ranging from 0.84 to 0.98, indicating that the unpaired spin is highly localized
on the Si centre. Its average Mulliken charge is +1.08 |e|, relatively lower than the
average Mulliken charge of bulk, 4-coordinated Si sites from 116 models of defect-free
a-SiO2 of +1.42 |e|, indicating that this 3-coordinated Si centre is more negative than
the average Si site in a-SiO2 . The nascent Si dangling bond introduces a one-electron
state 3.12 eV above the a-SiO2 valence band on average, almost resonant with the
top of the Si valence band (c.f. Si/SiO2 valence band offset 107 ). Fig. 8.4 shows the
Page 136
Figure 8.5: Correlations between the one-electron levels and relative stabilities of the
hydroxyl E0 centre with the non-bonding Si– –O interaction. (a) A histogram of the
distribution of the Si– –O distances of the dissociated Si–O bond in the hydroxyl E0
centre [see Fig. 8.2b]; (b) energies of the one-electron defect levels with respect to the
top of the a-SiO2 valence band plotted against the Si– –O distance; (c) the relative
stability of the hydroxyl E0 with respect to an interstitial H atom, plotted against the
Si– –O distance.
position of the occupied and unoccupied one-electron levels of the hydroxyl E0 centre
within the a-SiO2 band gap. It is interesting to note that from the 116 configurations,
a normal distribution of electronic states emerges.
The non-bonding Si– –O distances correlate strongly with some of the defect’s
properties. Fig. 8.5a shows a histogram of the non-bonding Si– –O distance. The
position of the one-electron defect level with respect to the top of the valence band
is plotted against this distance and is shown in Fig. 8.5b. The average position of
the one-electron level is 3.1 eV above the a-SiO2 valence band. The further away the
Si and the negative O ion are from each other, the deeper (closer to the top of the
valence band) the defect level becomes. This may be due to the reduction in repulsion
between the unpaired electron and the O ion. This reduction in repulsion energy also
Page 137
increases the relative stability of the hydroxyl E0 centre with respect to an interstitial
H atom, as seen in Fig. 8.5c.
The hydroxyl E0 centre is by far more thermodynamically stable than an interstitial
H atom, being lower in energy by 0.9 eV on average and ranging from a minimum of 0.3
eV to a maximum of 2.3 eV lower than the interstitial H atom. We have calculated
the formation energies of 50 configurations each of the hydroxyl E0 centre and an
interstitial H atom with respect to the Fermi level of the system according to equation
B.1 in the appendix. The configurations utilized in these calculations are shown in Fig.
8.6a. Note that the term interstitial was chosen purely to maintain consistency with
previous studies of hydrogen in SiO2 , despite the fact that the charged variants are
interesting defects in their own right, but a detailed discussion on them is beyond the
scope of this thesis. Remarkably, we find a number of configurations of the hydroxyl
E0 centre are thermodynamically stable in the the neutral charge state over a range
of Fermi levels (an example of which is shown in Fig. 8.6b), in stark contrast to what
has been calculated for interstitial hydrogen in the literature. 172,174 We find that 15 of
the 50 configurations studied are more stable as the neutral hydroxyl E0 centre across
a range of Fermi levels.
Thermodynamic switching levels are the Fermi levels for which the formation
energies of two different charge states of a defect are equal. For example, the (+/0)
thermodynamic switching level is the Fermi level for which the formation energies
of the positive and neutral charge states are the same. A histogram of these levels
is shown in Fig. 8.7. For the configurations in which the (+/-) switching level was
lower, only this value was recorded; however, in the systems where (+/0) is lower,
the (+/0) and (0/-) levels were recorded. This histogram shows that the neutral
hydroxyl E0 centre can be stable over a wide range of Fermi levels, from ≈ 3.5 to just
under 6 eV above the a-SiO2 valence band in our calculations. Note that it is not
necessarily a single configuration that is stable over this range. This range overlaps
strongly with the position of the Si band gap, indicating that the hydroxyl E0 centre
will contribute to the more thermodynamically stable hydrogen interactions in the
Page 138
Figure 8.6: Formation energy of the hydroxyl E0 centre in a-SiO2 . a) Configurations

of the hydroxyl E0 centre and interstitial hydrogen atom used in the calculations of
the formation energies. The configurations are labelled as 0, +, and - for the neutral,
positive and negative charge states respectively, with D and I standing for defect and
interstitial, respectively. b) Plot of the formation energy versus the Fermi level (with
respect to the a-SiO2 valence band) of a single configuration of the hydroxyl E0 centre
and an interstitial H atom. Formation energies corresponding to the hydroxyl E0
centre are drawn as solid red lines while the interstitial H atom’s formation energies
are drawn as dashed black lines. This diagram clearly shows that the neutral hydroxyl
E0 centre (lower horizontal solid red line) is the most stable charge state across a range
of Fermi levels.
Page 139
Si/SiO2 systems widely used in microelectronic devices. We have also calculated

the hyperfine splitting due to the unpaired electron’s interaction with the matrix in
13 configurations. There is a strong isotropic hyperfine splitting resulting from the
interaction of the 3-coordinated 29 Si with the localized electron which averages at 48.4
mT and ranges from 44.1 to 52.2 mT. This value shows a similarity to the experimental
isotropic hyperfine splitting for the E0 centre in a-SiO2 which was measured as 42.0
mT. 101
Figure 8.7: Histogram of the thermodynamic switching levels of the hydroxyl E0

centre. This shows that the neutral hydroxyl E0 centre can be stable over an ≈ 2.5
eV range of the Fermi level. The position of the Si band gap is the shaded section.
To investigate the formation of the hydroxyl E0 centre, the barrier from an intersti-
tial H to the defect configuration was calculated using the CI-NEB method described
in chapter 2.6 in 13 different models. A schematic of this reaction along with the cal-
culated barriers is shown in table 8.1 labelled as Reaction 1. The barriers range over
1.0 eV due to the slight variations in geometry of each amorphous structure, clearly
emphasising the importance of obtaining statistics when modelling amorphous SiO2
and demonstrating that there are some configurations where it will be much more
likely for an atomic H to create this defect configuration. We note that the distri-
bution of calculated barriers is highly asymmetric and does not result in a normal
distribution. This is again due to the variations in the a-SiO2 structures, where each
local chemical environment modulates the calculated barrier.
Passivation of the E0 centre and 3-coordinated Si defects in a-SiO2 is experimen-
Page 140
tally known to occur, 165 with the electron spin resonance (ESR) signal of the E0 centre
significantly reduced after a soak in an ambient of H2 or forming gas. 178 We investi-
gated the passivation of the hydroxyl E0 centre in the presence of atomic hydrogen. To
study the passivation, an H atom was placed 1 Å away from the undercoordinated Si
defect centre in 26 different systems and their total energies were optimized. We find
that this defect is indeed passivated by H, forming a stable Si–H bond which averages
at 1.45 Å, ranging from 1.43 Å to 1.46 Å. The defect level that was ≈ 3 eV above the
valence band in the neutral configuration is now suppressed and sits instead at the
top of the a-SiO2 valence band. The binding energy of the Si–H bond was calculated
as:
T ot T ot
EBinding = EInterstitial − ESi–H . (8.1)
T ot
where EInterstitial is the total energy of a system with the 3-coordinated Si defect and
T ot
an interstitial H atom sat in an Si–O void and ESi–H is the total energy of the system
after the defect has been passivated. The calculated binding energy for the nascent
Si–H bond reveals that it is strong, averaging at 4.19 eV and ranging from 3.96 to
4.31 eV calculated from 13 systems. Reaction 2 in table 8.1 shows this reaction and
its binding energies. The Si–H bond in a-SiO2 seems to have a narrow range of bond
lengths and binding energies in contrast to the energies associated with the Si and
O atoms, such as the position of the defect level in the band gap and even the Si–O
bond length distributions. This suggests that the Si–H bond formed is very stable
and that the passivation is very effective.
8.3.2 Molecular hydrogen in a-SiO2
In this section we look at how molecular hydrogen interacts with a-SiO2 . In particular,
we focus on two specific cases. The first is the introduction of an H2 molecule into
the melting phase when generating the a-SiO2 models. Introduction of H2 generates
Si–H bonds which can lead to the formation of defects. The second case we look at
is the interaction of a hydrogen molecule with defect-free a-SiO2 lattices. Again, this
creates Si–H bonds which act as precursors to defect configurations. These two cases
Page 141
Figure 8.8: Atomic structure and spin density of the neutral E0 centre. The O–H
bond that remains in the a-SiO2 structure is highlighted on the right hand side of the
figure. It is located 6 Å away from the 3-coordinated Si.
are described in detail below.
Neutral E0 Centres from Si-H Bonds
As discussed in the calculation details, twenty models of a-SiO2 were generated with
an extra H2 molecule added during the melt procedure. This resulted in each model
containing a single Si–H and O–H bond, as can be seen in Fig. 8.1. The Si–H
bonds in these models introduce a state that sits near the top of the a-SiO2 valence
band, as can be seen in the density of states shown in Fig. 8.1. Starting from these
optimized models, a H atom was removed from the Si–H bond in each structure and
their geometries were reoptimized. This resulted in a single electron residing in a
dangling bond at what is now a 3-coordinated Si, as illustrated in Fig. 8.8. This
defect structure is very similar to what is known as the E0 centre 166 and the system
is charge neutral, therefore it is referred to as the neutral E0 centre herein. As can
be seen in Fig. 8.8, the spin density is highly localized in a dangling bond at this 3-
coordinated Si. The 3-coordinated Si introduces a state in the a-SiO2 band gap, which
sits 3.0 eV above the valence band on average from the 20 models. Its average Mulliken
charge from the 20 different models is +1.07 |e|, significantly lower than the average
of +1.42 |e| of bulk Si atoms in these models. Its average Mulliken spin moment taken
Page 142
over the 20 models is 0.89 and ranges from 0.88 to 0.92, further confirming that the
electron is indeed highly localized at the 3-coordinated Si. The isotropic hyperfine
splittings were calculated in 8 models. The most dominant contribution comes from
the 3-coordinated Si, with a splitting that averagest at 44.4 mT and ranges from 40.0
to 47.8 mT.
We now turn to how the aforementioned defect configuration can be formed start-
ing from an Si–H bond. The bond would need to break and the H atom would diffuse
away. This reaction can be described by the following equation.
O3 Si H O3 Si + H0 (8.2)
The formation energy of this reaction was calculated as the total energy difference
between the reactants and the products.
F or T ot T ot
EN eutralE 0 = ESi−H − E3cSi (8.3)
F or T ot
Where EN eutralE 0 is the formation energy, ESi−H is the total energy of the system
T ot
that contains an Si–H bond and E3cSi is the total energy of the system when an H
atom is moved away from the O3 Si moiety. The H atom that was moved sits in
an interstitial position as described in section 8.3.1. The average formation energy
obtained from the 20 a-SiO2 models is 4.2 eV, which indicates that the Si–H bond is
very stable and would require high temperature for thermally activated dissociation
to occur.
As the Si–H states sit near the top of the valence band, we then investigated
whether hole trapping would facilitate the Si–H dissociation in all 20 a-SiO2 models,
as in reaction:
O3 Si H + h+ O3 Si + H+ (8.4)
where h+ is a hole. In the neutral system the Si–H states are located close to the
Page 143
Figure 8.9: The neutral E0 centre formed by trapping a hole at an Si–H bond. An
unpaired electron localised on a 3-coordinated Si atom and a proton bound to a
bridging oxygen can be clearly seen in the figure. The H in this figure is a proton.
top of the SiO2 valence band. In two of the 20 models, addition of a hole results
in spontaneous dissociation of a Si–H bond releasing a proton and leaving behind a
neutral E0 centre. The proton binds to the nearest oxygen that is not bonded to the Si
atom from which the proton dissociated, as can be seen in Fig. 8.9. The Si–H states
of the two models in which the Si–H bond dissociated spontaneously are the highest
in the a-SiO2 band gap, sitting slightly above the top of the SiO2 valence band. These
are also the configurations that have the smallest H– –O distance with the O atom
which binds the proton, as can be seen in Fig. 8.10. In the remaining 18 models
there is a barrier to remove a proton, but the final defect state is always lower in
energy than the Si–H state just after the hole trapping. The barriers for removing the
H+ after hole trapping was calculated using CI-NEB. The barrier to proton removal
increases as the H. . .O distance increases. In our 20 models, this barrier does not
exceed 0.5 eV even when the proton has to cross the largest H. . .O distance of 3.2
Å.
In all cases the proton binds to an oxygen atom that sits across an Si–O void,
forming a hydronium-like configuration shown in Fig. 8.9 and leaving behind a 3-
Page 144
Figure 8.10: A graph showing how the energetic position of the Si–H state, mea-
sured with respect to the top of the a-SiO2 valence band, changes with the nearest
non-bonding oxygen distance. There is a clear anti-correlation, with the defect level
moving further into the valence band as the nearest non-bonding oxygen moves further
away.
coordinated Si atom with an electron localized on this Si atom. Moving the proton
to one of the O atoms belonging to the same tetrahedron as the 3-coordinated Si is
also energetically favourable, but requires overcoming a higher barrier of 0.6 eV, as
opposed to the 0.5 eV maximum required to cross a ring. The average relaxation
energy, from our 20 systems, for a proton dissociating and binding to an oxygen atom
across the ring is 0.66 eV. However, an electron can be trapped by an E0 centre with
a proton less than ≈ 3 Å away, leading to restoration of the Si–H bond. When the
proton is further away, the Si–H bond does not reform.
Although the majority of the hyperfine splitting comes from the interaction of the
unpaired electron with the nucleus of the Si atom on which it is localized (see Fig.
8.9), the hyperfine interaction with the proton is also possible if the Si. . .H distance
is small enough. If the Si. . .H distance is 2.3 Å or less, this hyperfine splitting is 1.6
mT on average.
Comparison of the one-electron states of the E0 centre with the proton less than
Page 145
Figure 8.11: The effect of a nearby proton on the defect levels of the neutral E0 centre.
a) The defect levels of the neutral E0 centre. b) The defect levels of the neutral E0
centre when a proton is placed less than 6 Å away. The defect’s levels are shifted
closer to the valence band when a proton is nearby.
6 Å away and a neutral E0 centre with no proton in the system reveals that the
proton can shift the defect states down by up to ≈ 0.5 eV. That is, the defect state
is highly dependent on the proton’s location. Fig. 8.11 shows the defect levels for
a system containing a O3 Si moiety and the same system with a proton located
less than 6 Å away. It can clearly be seen that the defect levels of the system with a
proton are shifted down. The position of the proton depends on how easily the proton
can diffuse through a-SiO2 . It has been suggested by Pasquarello et al. that proton
diffusion occurs predominantly via a ring-crossing mechanism. 179 This involves the
proton hopping from one oxygen to another across an Si-O ring. Barriers for this
hopping range from almost zero to ≈ 1.5 eV, dependent on the distance between
the initial and final O atoms. We also find that proton diffusion proceeds via ring-
crossing with the E0 centre making little difference to the energies previously calculated
for proton diffusion. 179 This means that the defect’s level will be modulated by the
protons that exist nearby.
Page 146
Hydroxyl E0 Centre from H2 Molecules
We now turn to the question of whether interaction of a H2 molecule with a pristine,

stoichiometric a-SiO2 matrix can result in a passivated configuration of the hydroxyl
E0 centre without initial creation of the active defect. In 10 different configurations
of bulk a-SiO2 , H2 molecules were placed in interstitial positions and their geometries
were optimized. Hydrogen molecules sit in the middle of an Si–O void and interact
negligibly with the a-SiO2 matrix, similar to the interstitial H atoms described earlier
in section 8.3.1 and in good agreement with what has been previously calculated in
the literature. 86,180 Passivated defect structures were then generated by forming an
Si–H bond and an O–H bond facing each other at a long Si–O bond and the geometry
of each configuration was then optimized, as described in section 8.3.1. We find a
variation in the qualitative stabilities of the passivated defect structures, being more
thermodynamically favourable than interstitial H2 molecules by 0.2 eV on average,
but ranging from 0.61 more to 0.29 eV less stable. The barrier to form the passivated
configuration from an interstitial H2 molecule was then calculated using the CI-NEB
method with the interstitial H2 and passivated configuration described above used as
the initial and final images respectively. A summary of the reaction and its barriers
are presented as Reaction 3 in table 8.1. The barriers for the forward reaction average
at 1.74 eV while they average at 1.94 eV for the reverse reaction. The forward barrier
is relatively high and would require a high temperature to overcome and form a
passivated hydroxyl E0 centre, but the reverse barrier indicates that once it is formed
it would be kinetically very stable.
H2 Cracking at Hydroxyl E0 Centres
Previous theoretical studies have investigated the reaction kinetics of a H2 molecule

cracking at an E0 centre, passivating the defect and leaving behind an interstitial H
atom. 175,176,178 This reaction has been shown to be endothermic at multiple levels
of theory, i.e., the passivated configuration is thermodynamically unfavourable with
respect to the active defect and an interstitial H2 molecule. We have performed sim-
Page 147
Forward Reactions Reverse Reactions ∆E

Reaction Min. Max. Avg. Min. Max. Avg.
Hydrogen reactions in defect free a-SiO2
O HO
1) O3 Si Si O3 + H0 O3 Si Si O3 0.50 1.71 0.91 1.25 2.40 1.66 -0.75
HO HO
2) O3 Si Si O3 + H0 O3 Si H Si O3 0.00 0.00 0.00 3.96 4.31 4.19 -4.19
O HO
3) O3 Si Si O3 + H2 O3 Si H Si O3 1.07 2.15 1.74 1.57 2.45 1.94 -0.20
HO HO
4) O3 Si Si O3 + H2 O3 Si H Si O3 + H0 0.46 0.78 0.65 0.10 0.24 0.20 0.45
Table 8.1: Hydrogen’s reactions with a-SiO2 . All results are in eV.
ilar calculations at the hydroxyl E0 centre and a schematic is shown as Reaction 4 in

table 8.1. Starting from the neutral configuration of the hydroxyl E0 centre, an H2
molecule was placed within 2 Å of the defect and the geometry of the system was
optimized. The H2 remains in an interstitial position and interacts with neither the
defect nor the a-SiO2 matrix. Using this configuration and the passivated configu-
ration with an extra interstitial H atom, we used CI-NEB to calculate the barrier
for the passivation of the undercoordinated Si by an H2 molecule in 13 models. The
calculated barriers are qualitatively similar to previous calculations of H2 dissociation
at E0 centres, 175,176,181 with the defect configuration in the presence of a H2 molecule
being thermodynamically favourable by 0.45 eV on average, ranging from 0.22 to 0.68
eV more stable. The barrier for passivation from an interstitial H2 molecule averages
at 0.65 eV, ranging from 0.46 eV to 0.78 eV. However, the reverse barrier is much
lower and averages at 0.20 eV, ranging from 0.10 eV to 0.24 eV. Similar to the Si–H
binding energies, these barriers exhibit a much narrower range and seem to be more
resistant to fluctuations in the local environment of the a-SiO2 matrix.
Our modelling confirms that Si–H bonds in a-SiO2 can act as precursors to formation
of neutral E0 centres and that the presence of strained Si–O bonds in a-SiO2 gives rise
to an additional channel of interaction of H atoms with a-SiO2 networks, predicting
the formation of a hydroxyl E0 centre.
The Si–H dissociation is facilitated by hole injection and the barrier to Si–H dis-
Page 148
sociation is at most 0.5 eV in our models. Proton diffusion away from the neutral E0
centre occurs via a ring crossing mechanism with barriers less than 1.5 eV. Electron
injection can restore the Si–H bond if the proton is less than 3 Å away. The hyperfine
interaction of the neutral E0 centre is in good agreement with the experimental values.
This defect introduces a one-electron level that sits 3.0 eV on average above the SiO2
valence band. This position makes this defect a candidate for the defects that could
be involved in electronic device reliability issues. An interesting aspect of this model is
the 1.6 mT signal associated with the proton. Experimentally there is a weak satellite
signal of 1.3 mT associated with the E0 centre controversially attributed either to a
29
Si atom 182 or hydrogen atoms. 183 Our calculations suggest that this satellite signal
may be related to hydrogen present in the system after Si–H dissociation. Experi-
mentally the 1.3 mT signal has been seen to disappear at temperatures greater than
100◦ C while the E0 signal is increasing. 184 In this model the proton can diffuse away
after overcoming some barrier, presumably leading to a disappearance of the 1.6 mT
signal, even while neutral E0 centres are being generated.
The hydroxyl E0 was shown to have deep levels in the a-SiO2 band gap, which are
nearly resonant with the top of the Si valence band. Some configurations are stable in
the neutral charge state over a range of Fermi level positions in the a-SiO2 band gap.
These results demonstrate that the hydroxyl E0 centre could have an effect on the
technological applications of a-SiO2 . Its defect levels are located close to the top of
the Si valence band in a Si/SiO2 system (see Fig. 8.4), typically used in metal-oxide-
semiconductor (MOS) electronic devices. Holes at the top of the Si valence band,
typical of a pMOS device, can tunnel into the hydroxyl E0 centre so that it can act as
a source of trapped charge in such a device. Hence H0 is not always a benign agent
in defect-free a-SiO2 networks and can produce thermodynamically stable neutral
defects in a-SiO2 , adding to the density of dangling bond defects, such as E0 centres,
which are implicated in reliability issues of devices which utilize a-SiO2 . The results
presented for atomic hydrogen’s interactions with a-SiO2 may also shed new light on
the behaviour of atomic hydrogen in other amorphous solids, in which H0 is thought
Page 149
to interact negligibly. 172

Both the neutral E0 and hydroxyl E0 centres have defect levels which sit in similar
positions in the a-SiO2 band gap, ≈ 3.5 eV above the valence band. In an Si/SiO2
system, this state would sit near the top of the Si valence band and holes could
potentially tunnel into these states. Thus they would be detrimental to the operation
of an electronic device made using this system. Furthermore, the position of the
defect levels suggests that they could absorb radiation in the ultra-violet range of
the spectrum. Further calculations would need to be performed to calculate their
optical absorption spectra and assess how detrimental these defects would be to optical
devices.
The reactions of hydrogen species with defect-free a-SiO2 are summarized in table
8.1. Our results for H2 dissociation at the hydroxyl E0 centre can explain the slow
passivation kinetics of the E0 centre signals under an H2 atmosphere. The barrier to
H2 cracking and passivation is quite high while the final configuration is thermody-
namically unfavourable in the presence of an H atom, indicating that the passivation
reaction would proceed rather slowly. It is also important to note the consequence
of the backwards reactions, which essentially re-activates the passivated defect. The
barriers for the de-passivation reactions are rather low and, if a source of atomic H
exists in the system, it can de-passivate Si–H bonds and activate a hydroxyl E0 centre.
Notably, the de-passivation reactions studied are also similar to the experimentally ob-
served de-passivation of the Pb centre, a defect consisting of a 3-coordinated Si with a
trapped electron, which exists at the Si/SiO2 interface. 185 The de-passivation reaction
is insensitive to how the passivated configuration was generated in the first place and
is only limited by the concentration and diffusion of atomic H in the system. Atomic
H can be released into the a-SiO2 layer in metal-oxide-semiconductor (MOS) devices
in a number of conditions. 186–188 Coupled with the results for the de-passivation reac-
tions and the positions of the hydroxyl E0 defect levels, these results strongly suggest
that the hydroxyl E0 can be a potentially detrimental defect for electronic and optical
technologies.
Page 150
Concluding Remarks
9
9.1 Summary
A number of novel charge trapping defects in SiO2 were investigated in this thesis.
They were shown to strongly affect their local atomistic surroundings and significantly
change the electronic structure of the oxide. Defects such as this can significantly
affect the operation of an electronic device and their electronic structures show that
they could be involved in device reliability issues.
Extrinsic sites in SiO2 were shown to trap electrons. Substitution of an Si atom
with Ge in α-quartz results in a GeO4 tetrahedral structure, similar to the SiO4
tetrahedral unit but with longer Ge–O bonds. This Ge atom is associated with a state
which sits just below the α-quartz conduction band in the neutral system. Addition
Chapter 9. Concluding Remarks
of an electron causes a significant structual rearrangement around the Ge impurity.

Two of the Ge–O bonds extend significantly and the angle between these two bonds is
severely distorted (150◦ ) from the ideal tetrahedral angle (109.5◦ ). The extra electron
was found to localize on this wide angle. The state on which the electron localises
lies deep in the α-quartz band gap, 4.2 eV below the conduction band, making it a
deep electron trap. An interstitial Li atom in α-quartz was also investigated. In the
neutral charge state, it was found to sit equidistantly between two SiO4 tetrahedra.
The Li atom introduces an occupied state which sits just below the bottom of the
α-quartz conduction band. However, opening a nearby O–Si–O angle resulted in an
electron trapped at that site, similar to the that of Ge-doped α-quartz. A defect state
sits 2.5 eV below the conduction band and an electron was shown to have transferred
onto an Si site. The calculated hyperfine splittings of the trapped electron are in good
agreement with the experimental signals of the [SiO4 /Li]0 centre in α-quartz. 118 This
was the first time that the atomistic structure of the [SiO4 /Li]0 was calculated.
Electrons were also shown to spontaneously trap at intrinsic sites in a-SiO2 . These
traps create a defect state which sits 3.2 eV below the conduction band on average.
The precursor sites at which electrons could trap were shown to be wide O–Si–O
angles. By constructing a model of a-SiO2 which contains almost half a million atoms,
the concentration of these precursor sites was estimated as 4 ×1019 cm−3 . Although
this concentration is rather high when compared to that of a typical point defect in
a-SiO2 , for example, the oxygen vacancy, experimental results on SiC/SiO2 electronic
devices show a similar concentration of electron traps which are not attributed to any
existing defect. The optical absorption spectra of these electron traps were calculated
in an embedded cluster model using TD-DFT for the first time. The calculated spectra
show a strong absorption at 3.7 and 5.8 eV. These peaks are in good agreement with
experimental measurements at low temperatures, providing the strongest evidence yet
for intrinsic electron trapping in a-SiO2 . The position of the defect levels suggest that
they could indeed play a role in electronic device reliability issues.
Finally, Hydrogen’s interactions with SiO2 were investigated as it is widely used in
Page 152
electronic device fabrication processes. It was shown that both atomic and molecular
hydrogen can be involved in defect generation processes. Introduction of molecular
hydrogen into the melted phase of SiO2 followed by quenching results in Si–H bonds.
They were shown to interact with holes to generate a neutral 3-coordinated Si defect
which is referred to as the neutral E0 centre. During this process, a proton is liberated
and binds to a nearby bridging O. The distance between this liberated proton and
the 3-coordinated Si has a strong effect on the defect’s level in the SiO2 band gap.
Atomic hydrogen was also shown to generate defects. In particular, it can interact with
bridging O at long Si–O bonds to generate a defect which has not been previously
characterized in the literature. It consists of a 3-coordinated Si facing a hydroxyl
group and is referred to as the hydroxyl E0 centre. Both the neutral and hydroxyl
E0 centres have defect levels which suggest that they can both be involved in device
reliability issues.
9.2 Future Work
This thesis investigated a number of defects in SiO2 using hybrid functionals and
DFT. Despite being characterised extensively, this work can be extended in a number
of directions.
A common theme of the intrinsic electron traps (see chapter 6) and the hydroxyl E0
centre (see chapter 8) is that the disorder of a-SiO2 allows for creation of novel defects.
Many other amorphous materials are now being considered for use in electronic and
other technological applications. For instance, hafnium dioxide is now regularly used
in electronic devices as a replacement for SiO2 , and studies show that it is deposited
in its amorphous state. 189,190 Understanding the role of disorder in defect creation
in other materials used in technological applications would help in producing more
reliable devices.
It is important to consider what nuclear quantum effects may have on the results
presented for the hydrogen related defects in chapter 8. These effects are known to
Page 153
play a large role in systems such as water and ice due to the small mass of the hydrogen
nuclei. They have been shown to have an effect on calculated barriers for hydrogenic
processes in SiO2 175 and are bound to have an effect on the results presented in chapter
8. A standard approach to investigating these effects would be to use a path integral
formulation of quantum mechanics, as has been reported in the literature. 191
Optical absorption spectra of intrinsic electron trapping in SiO2 was calculated
using TD-DFT in chapter 6. Although TD-DFT has had numerous successes in
describing excitations in a range of materials, it has trouble describing anything higher
than single excitations. This may have an effect on the calculated absorption spectra
presented in that chapter.
The defects investigated in this thesis could be involved in electronic device reli-
ability issues. However, further investigation is required to assess exactly what kind
of role they may play. A more direct method of analysing their role would be to
calculate charge transfer rates between the defects and a charge reservoir in the ma-
terial used to fabricate the device. As experimental techniques improve, more and
more data exists for these charge transfer rates. Being able to calculate these rates
and compare directly to experiment would be a very useful tool in implicating such
defects in electronic device reliability issues.
Furthermore, these defects are implicated in reliability issues where SiO2 is inter-
faced with another material. The most widely used MOS system centres around the
Si/SiO2 interface, which has been modelled in chapter 3. The point defects in the
oxide will be affected by the substrate, therefore studying how these defects would
change in an interface system would provide valuable information to the device relia-
bility community.
Page 154
BKS Parameterisation for Si and O
A
i−j Aij / eV bij / Å−1 cij / eVÅ−6 Atomic Charges

O–O 1388.7730 2.76000 175.0000 qO =-1.2
Si–O 18003.7572 4.87318 133.5381 qSi = +2.4
Table A.1: BKS Parameterisation for Si and O atoms used in this work.
Formation Energy Calculation
B
To assess the thermodynamic stability of the defects studied, their formation energies
were calculated as:
Ef or (F ) = Edef ect − (Ebulk + EH 0 ) + q(F + ∆V ) + Ecorr , (B.1)
where Edef ect is the total energy of the defective system, Ebulk is the energy of the
defect-free system, EH 0 is the energy of a H atom calculated using the same method,
q is the charge state of the defect, F is the Fermi level referenced to the top of the
a-SiO2 valence band, ∆V is a potential alignment term, and Ecorr is a correction term
for the periodic interaction between localized charges in the charged systems. The
∆V term was found to be negligible (< 0.05 eV) and was therefore ignored.
Chapter B. Formation Energy Calculation
B.1 Charge Correction
The Lany and Zunger method for charge correction was chosen for its ability to de-
scribe the interaction between a localized charge and extended delocalized screening
charge density, which comes out of DFT calculations of charged cells. 123,124 The ana-
lytic form of the charge correction is the same for all the defects, irrespective of the
character of localization, and is calculated as:
2
π 1 q α
Ecorr = 1− 1− , (B.2)
3α ε 2εL
where ε is the macroscopic dielectric constant of SiO2 (3.9 192 ), q is the charge of the
cell, α is the Madelung constant for a single charge in a periodic array and L is the
supercell length.
Page 157
Bibliography
[1] C. E. Association, Smartphones, HDTVs Are the Most Planned CE Pur-

chases This Year, According to CEA Study, http://www.ce.org/News/News-
Releases/Press-Releases/2012-Press-Releases/Smartphones,-HDTVs-Are-the-
Most-Planned-CE-Purchas.aspx, Accessed June 24, 2015.
[2] A. C. Diebold, D. Venables, Y. Chabal, D. Muller, M. Weldon and E. Garfunkel,

Mater. Sci. Semicond. Process., 1999, 2, 103–147.
[3] C. Hamaguchi, in Basic semiconductor physics, Springer, Berlin, 2010, pp. 334–
335.
[4] S. M. Sze and K. K. NG, in Physics of semiconductor devices Third edition,

Wiley, New Jersey, 2007, pp. 293–319.
[5] A. Stesmans and V. V. Afanas’ev, J. Phys.: Condens. Matter, 1998, 10, L19.
[6] A. Kawamoto, Ph.D. thesis, Stanford Universtity, 2001.
[7] J. H. Stathis and D. J. DiMaria, IEDM, 1998, 167–170.
[8] R. Degraeve, G. Groeseneken, R. Bellens, M. Depas and H. E. Maes, IEDM,

1995, 863–866.
[9] J. H. Stathis and S. Zafar, Microelectron. Reliab., 2005, 46, 270–286.
[10] Y. Miura and Y. Matukura, Jpn. J. Appl. Phys., 1966, 5, 180.

Chapter Bibliography
[11] P. Blöchl and J. Stathis, Phys.Rev.Lett., 1999, 83, 372–375.
[12] J. H. Stathis and E. Cartier, Phys. Rev. Lett., 1994, 72, 2745.
[13] E. Cartier, J. H. Stathis and D. Buchanan, Appl. Phys. Lett., 1995, 63, 1510.
[14] R. E. Stahlbush, E. Cartier and D. Buchanan, Microelectron. Eng., 1995, 28, 3.
[15] F. Schanovsky, W. Gös and T. Grasser, J. Vac. Sci. Technol. B, 2011, 29,
01A201.
[16] G. Pacchioni and C. Mazzeo, Phys. Rev. B, 2000, 62, 5452.
[17] L. Verlet, Phys. Rev., 1967, 159, 98–103.
[18] P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
[19] W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
[20] A. Dreuw, J. L. Weisman and M. Head-Gordon, J. Chem. Phys., 2003, 119,

2943–2946.
[21] S. Hirata and M. Head-Gordon, Chem. Phys. Lett., 1999, 302, 375 – 382.
[22] F. Furche and R. Ahlrichs, J. Chem. Phys., 2002, 117, 7433–7447.
[23] E. Runge and E. K. U. Gross, Phys. Rev. Lett., 1984, 52, 997–1000.
[24] M. Marques and E. Gross, Annu. Rev. Phys. Chem., 2004, 55, 427–455.
[25] G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–487.
[26] N. Troullier and J. L. Martins, Phys. Rev. B, 1991, 43, 1993.
[27] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinelo, T. Chassaing and

J. Hutter, Comp. Phys. Comm., 2005, 167, 103.
[28] J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip.

Rev. Comput. Mol Sci., 2014, 4, 15–25.
Page 159
[29] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 1980, 45, 566–569.
[30] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
[31] J. P. Perdew and M. Levy, Phys. Rev. Lett., 1983, 51, 1884.
[32] J. L. Gavartin and A. L. Shluger, Rad. Eff. Defects Solids, 2001, 155, 311–315.
[33] R. Dovesi, R. Orlando, C. Roetti, C. Pisani and V. Saunders, Phys. Status Solidi
B, 2000, 217, 63–88.
[34] M. Guidon, J. Hutter and J. Vandevondele, J. Chem. Theory Comput., 2009,

5, 3013–3021.
[35] C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785.
[36] A. Nakano, R. K. Kalia, K. ichi Nomura, A. Sharma, P. Vashishta, F. Shimojo,

A. C. van Duin, W. A. Goddard, R. Biswas and D. Srivastava, Comp. Mater.
Sci., 2007, 38, 642 – 652.
[37] F. F. Abraham, R. Walkup, H. Gao, M. Duchaineau, T. D. De La Rubia and

M. Seager, Proc. Natl. Acad. Sci. U.S.A, 2002, 99, 5777–5782.
[38] B. W. H. van Beest, G. J. Kramer and R. A. van Santen, Phys. Rev. Lett., 1990,
64, 1955–1958.
[39] D. Herzbach, K. Binder and M. H. Müser, J. Chem. Phys., 2005, 123, 124711.
[40] T. Watanabe, H. Fujiwara, H. Noguchi, T. Hoshino and I. Ohdomari, Jap. J.

Appl. Phys., 1999, 38, L366.
[41] A. K. Rappé and W. A. Goddard, J. Phys. Chem., 1991, 95, 3358–3363.
[42] R. T. Sanderson, Chemical bonds and bond energy, Academic press, 1976.
[43] T. R. Lucas, B. A. Bauer and S. Patel, Biochim. Biophys. Acta, 2012, 1818,
318 – 329.
Page 160
[44] S. Patel, A. D. Mackerell and C. L. Brooks, J. Comp. Chem., 2004, 25, 1504–
1514.
[45] T. Shan, B. D. Devine, J. M. Hawkins, A. Asthagiri, S. R. Phillpot and S. B.

Sinnott, Phys. Rev. B, 2010, 82, 235302.
[46] T. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott and S. R. Phillpot, Phys.

Rev. B, 2010, 81, 125328.
[47] J. Ludwig, D. G. Vlachos, A. C. Van Duin and W. A. Goddard, J. Phys. Chem.

B, 2006, 110, 4274–4282.
[48] A. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu and W. Goddard, J.

Phys. Chem. A, 2001, 105, 9396–9409.
[49] A. Yasukawa, JSME It J., 1996, 39, 313.
[50] P. Richet, Y. Bottinga, L. Denielou, J. Petitet and C. Tequi, Geochim. Cos-

mochim. Acta, 1982, 46, 2639 – 2658.
[51] A. C. T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu and W. Goddard,

J. Phys. Chem. A, 2003, 107, 3803–3811.
[52] J. C. Fogarty, H. M. Aktulga, A. Y. Grama, A. C. T. van Duin and S. A. Pandit,

J. Chem. Phys., 2010, 132, 174704.
[53] A. C. T. van Duin, J. M. A. Baas and B. van de Graaf, J. Chem. Soc., Faraday
Trans., 1994, 90, 2881.
[54] K. J. Laidler and M. C. King, J. Phys. Chem., 1983, 87, 2657–2664.
[55] N. Govind, M. Petersen, G. Fitzgerald, D. King-Smith and J. Andzelm, Comp.

Mater. Sci., 2003, 28, 250 – 258.
[56] D. J. Wales, J. Chem. Phys., 1994, 101, 3750.
[57] J. Zhao, S. L. Simon and G. B. McKenna, Nat.Commun., 2013, 4, .
Page 161
[58] W. H. Zachariasen, J. Am. Chem. Soc., 1932, 54, 3841–3851.
[59] S. Mukhopadhyay, P. V. Sushko, A. M. Stoneham and A. L. Shluger, Phys. Rev.

B, 2004, 70, 195203.
[60] L. S. Wang, J. B. Nicholas, M. Dupuis, H. B. Wu and S. D. Colson, Phys. Rev.

Lett., 1997, 78, 4450–4453.
[61] S. Plimpton, J. Comp. Phys., 1995, 117, 1–19.
[62] J. R. Chelikowsky, Phys. Rev. B, 1998, 57, 3803–3811.
[63] J. R. Chelikowsky and X. Jing, Phys. Rev. B, 1994, 50, 3803–3811.
[64] M. Litniewski, Mol. Sim., 2003, 29, 223–229.
[65] A. R. Leach, in Molecular modelling: Principles and applications, Pearson Pren-

tice Hall, 2001, pp. 359–362.
[66] T. Gerber and B. Himmel, J. Non-Cryst. Solids, 1985, 83, 324–334.
[67] J. Neuefeind, K. D. Liss and B. Bunsenges, J. Phys. Chem., 1996, 100, 1341.
[68] K. Vollmayr, W. Kob and K. Binder, Phys. Rev. B, 1996, 54, 15808–15827.
[69] A. Roder, W. Kob and K. Binder, J. Chem. Phys., 2001, 114, 7602–7614.
[70] A. Navrotsky, Diffus. Defect Data, 1987, 61, 53–54.
[71] S. R. Elliott, in Physics of Amorphous Materials, Longman Scientific & Tech-

nical, New York, 1990, ch. 3.1, pp. 71–132.
[72] D. Price and J. Carpenter, J. Non-Cryst. Solids, 1987, 92, 153 – 174.
[73] M. Guidon, J. Hutter and J. VandeVondele, J. Chem. Theory Comput., 2010,

8, 2348–2364.
[74] R. A. B. Devine and J. Arndt, Phys. Rev. B, 1987, 35, 9376–9379.
[75] Y.-N. Xu and W. Y. Ching, Phys. Rev. B, 1995, 51, 17379–17389.
Page 162
[76] F. Urbach, Phys. Rev., 1953, 92, 1324–1324.
[77] P. W. Anderson, Phys. Rev., 1958, 109, 1492–1505.
[78] H. Akatsu, Y. Sumi, T. Iwabuchi, S. Kaneko, T. Ueno and I. Ohdomari, J.

Appl. Phys., 1992, 72, 1906.
[79] N. Miyata, H. Watanabe and M. Ichikawa, Phy. Rev. B, 1998, 58, 13670.
[80] S. Miyazaki, H. Nishimura, M. Fukuda, L. Ley and J. Ristein, Appl. Surf. Sci.,
1997, 113/114, 585.
[81] A. C. Diebold, D. Venables, Y. Chabal, D. Muller, M. Weldon and E. Garfunkel,

Mat. Sci. Semicon. Proc., 1999, 2, 103–147.
[82] K. Hirose, H. Nohira, T. Koike, K. Sakano and T. Hattori, Phys. Rev. B, 1999,
59, 5617.
[83] J. H. Oh, H. W. Yeom, Y. Hagimoto, K. Ono and M. Oshima, Phys. Rev. B,

2001, 63, 205310.
[84] D. M. Fleetwood, S. T. Pantelides and R. D. Schrimpf, Defects in microelecronic

materials and devices, CRC Press, 2009.
[85] D. L. Griscom, J. Non-Cryst. Solids, 2011, 357, 1945–1962.
[86] P. Blochl, Phys. Rev. B, 2000, 62, 6158.
[87] V. B. Sulimov, P. V. Sushko, A. H. Edwards, A. L. Shluger and A. M. Stoneham,

Phys. Rev. B, 2002, 66, 024108.
[88] P. V. Sushko, A. L. Shluger and C. R. A. Catlow, Surf. Sci., 2000, 450, 153–170.
[89] L. Levien, C. T. Prewitt and D. Weidner, Am. Mineral., 1980, 65, 920–930.
[90] S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B, 1996, 54, 1703–1710.
[91] J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
Page 163
[92] Z.-Y. Lu, C. J. Nicklaw, D. M. Fleetwood, R. D. Schrimpf and S. T. Pantelides,

Phys. Rev. Lett., 2002, 89, 285505–1–4.
[93] C. J. Nicklaw, Z. Y. Lu, D. M. Fleetwood, R. D. Schrimpf and S. T. Pantelides,

IEEE Trans. Nucl. Sci., 2002, 49, 2667 – 2673.
[94] M. Walters and A. Reisman, J. Electrochem. Soc., 1991, 138, 2756–2762.
[95] T. Grasser, B. Kaczer, W. Goes, H. Reisinger, T. Aichinger, P. Hehenberger,

P.-J. Wagner, F. Schanovsky, J. Franco, M. T. Luque and M. Nelhiebel, IEEE
Trans. Electron Devices, 2011, 58, 3652–3666.
[96] D. Donaldio, M. Bernasconi and M. Boero, Phys. Rev. Lett., 2001, 87, 195504–
1–4.
[97] P. V. Sushko, S. Mukhopadhyay, A. M. Stoneham and A. L. Shluger, Micro-

electron. Eng., 2005, 80, 292–295.
[98] A. S. Mysovsky, P. V. Sushko, S. Mukhopadhyay, A. H. Edwards and A. L.

Shluger, Phys. Rev. B, 2004, 69, 085202.
[99] R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980,

72, 650.
[100] A. D. McLean and G. S. Chandler, J. Chem. Phys., 1980, 72, 5639.
[101] M. G. Jani, R. B. Bossoli and L. E. Halliburton, Phys. Rev. B, 1983, 27, 2285.
[102] A. Stoneham, M. Szymanski and A. Shluger, Phys.Rev.B, 2001, 63, 241304.
[103] A. Bongiorno and A. Pasquarello, Phys.Rev.Lett., 2002, 88, 125901.
[104] Z. Celik-Butler, P. Vasina and N. Vibhavie Amarasinghe, IEEE Trans. Electron

Devices, 2000, 47, 646–648.
[105] F. Schanovsky, W. Gös and T. Grasser, J. Comput. Electron., 2009, 9, 135–140.
Page 164
[106] S. Mukhopadhyay, P. V. Sushko, A. M. Stoneham and A. L. Shluger, Phys. Rev.

B, 2005, 71, 235204–1–9.
[107] V. V. Afanasev, M. Bassler, G. Pensl, M. J. Schulz and E. Stein von Kamienski,

J. Appl. Phys., 1996, 79, 3108–3114.
[108] J. A. Weil, Phys. Chem. Miner., 1984, 10, 149–165.
[109] J. Isoya, J. A. Weil and R. F. C. Claridge, J. Chem. Phys., 1978, 69, 4876.
[110] D. L. Griscom, Opt. Mater. Express, 2011, 1, 400–412.
[111] A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
[112] J. P. Douglas, Semicond. Sci. Technol., 2004, 19, R75.
[113] J. H. Jang, L. Wantae, M. S. Phen, K. Siebein, S. J. Pearton and V. Craciun,

Appl. Phys. Lett., 2009, 94, 202104.
[114] I. M. Lee and C. G. Takoudis, J. Vac. Sci. Technol. A, 1997, 15, 3154.
[115] E. Long, A. Azarov, F. Klow, A. Galeckas, K. A. G. and S. Diplas, J. Appl.

Phys., 2012, 111, 024308.
[116] S. J. Kilpatrick, R. J. Jaccodine and P. Thompson, J. Appl. Phys., 2003, 93,

4896.
[117] D. C. Paine, C. Caragianis and A. F. Schwartzman, J. Appl. Phys., 1991, 70,

5076.
[118] M. G. Jani, L. E. Halliburton and A. Halperin, Phys. Rev. Lett., 1986, 56,
1392–1395.
[119] L. Levien, C. T. Prewitt and D. J. Weidner, Am. Mineral., 1980, 65, 920–930.
[120] B. Civalleri and P. Ugliengo, J. Phys. Chem. B, 2000, 104, 519–532.
[121] M. D. Towler, N. L. Allan, N. M. Harrison, V. R. Saunders, W. C. Mackrodt

and E. Apra, Phys. Rev. B, 1994, 50, 5041–5054.
Page 165
[122] R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980,

72, 650.
[123] S. Lany and A. Zunger, Modelling Simul. Mater. Sci. Eng., 2009, 17, 084002.
[124] H.-P. Komsa, T. T. Rantala and A. Pasquarello, Phys. Rev. B, 2012, 86, 045112.
[125] W. Hayes and T. J. L. Jenkin, J. Phys. C: Solid State Phys., 1986, 19, 6211–
6219.
[126] M. E. Markes and L. E. Halliburton, J. Appl. Phys., 1979, 50, 8172.
[127] T. M. Wilson, J. A. Weil and P. S. Rao, Phys. Rev. B, 1986, 34, 6053–6055.
[128] A. Sartbaeva, S. A. Wells and S. A. T. Redfern, J. Phys.: Condens. Matter,

2004, 16, 8173–8189.
[129] T. Y. Chan, K. K. Young and C. Hu, IEEE Electron Device Lett., 1987, 8,
93–95.
[130] G. Bersuker, A. Korkin, Y. Jeon and H. Huff, Appl. Phys. Lett., 2002, 80, 832.
[131] M. Farnesi Camellone, J. C. Reiner, U. Sennhauser and L. Schlapbach, Phys.

Rev. B, 2007, 76, 125205.
[132] L. Giordano, P. V. Sushko, G. Pacchioni and A. L. Shluger, Phys. Rev. Lett.,

2007, 99, 136801.
[133] V. V. Afanas’ev and A. Stesmans, Phys. Rev. Lett., 1997, 78, 2437.
[134] V. V. Afanas’ev and A. Stesmans, Appl. Phys. Lett., 1997, 70, 1260.
[135] V. V. Afanas’ev and A. Stesmans, Microelectron. Eng., 1997, 36, 149.
[136] V. V. Afanas’ev and A. Stesmans, J. Phys.: Condens. Matter, 1997, 9, L55.
[137] V. V. Afanas’ev, A. Stesmans, M. Bassler, G. Pensl and M. J. Schulz, Appl.

Phys. Lett., 2000, 76, 336–338.
Page 166
[138] N. S. Saks and A. K. Agarwal, Appl. Phys. Lett., 2000, 77, 3281–3283.
[139] N. S. Saks, S. S. Mani and A. K. Agarwal, Appl. Phys. Lett., 2000, 76, 2250–
2252.
[140] M. Städele, M. Moukara, J. A. Majewski and P. Vogl, Phys. Rev. B, 1999, 59,
10031.
[141] A. D. Becke, J. Chem. Phys. C, 1993, 98, 5648.
[142] D. L. Griscom, J. Non-Cryst. Solids, 2006, 352, 2601.
[143] I. Pintilie, C. M. Teodorescu, F. Moscatelli, R. Nipoti, A. Poggi, S. Solmi,

L. S. L. vlie and B. G. Svensson, J. Appl. Phys., 2010, 108, 024503.
[144] V. V. Afanasev, F. Ciobanu, S. Dimitrijev, G. Pensl and A. Stesmans, J. Phys.:

Condens. Matt., 2004, 16, S1839.
[145] V. V. Afanas’ev and A. Stesmans, Phys. Rev. B, 1999, 59, 2025–2034.
[146] G. Seguini, S. Schamm-Chardon, P. Pellegrino and M. Perego, Appl. Phys. Lett.,

2011, 99, 082107.
[147] A. Waseda and K. Fujii, IEEE Trans. Instrum. Meas., 2007, 56, 628–631.
[148] Y. Sasajima and K. Tanimura, Phys. Rev. B, 2003, 68, 014204.
[149] L. Skuja, J. Non-Cryst. Solids, 1998, 239, 16–48.
[150] A. V. Kimmel, P. V. Sushko, A. L. Shluger and G. Bersuker, ECS Transactions,

2009, 19, 3–17.
[151] A.-M. El-Sayed, M. B. Watkins, V. V. Afanas’ev and A. L. Shluger, Phys. Rev.

B, 2014, 89, 125201.
[152] S. Lany and A. Zunger, Phys. Rev. B, 2009, 80, 085202.
[153] T. J. L. Jenkin, J. Koppitz and W. Hayes, J. Phys. C: Solid State Phys., 1987,
20, L367–L371.
Page 167
[154] S. Girard, J. Kuhnhenn, A. Gusarov, B. Brichard, M. Van Uffelen, Y. Ouerdane,

A. Boukenter and C. Marcandella, IEEE Trans. Nucl. Sci., 2013, 60, 2015–2036.
[155] A. Alessi, S. Agnello, S. Grandi, A. Parlato and F. M. Gelardi, Phys. Rev. B,

2009, 80, 014103.
[156] A. N. Trukhin, M. N. Tolstoi, L. B. Glebov and V. L. Savelev, Phys. Stat. Sol.

B, 1980, 99, 155–162.
[157] R. H. Kars, N. R. J. Poolton, M. Jain, C. Ankjærgaard, P. Dorenbos and

J. Wallinga, Radiat. Meas., 2013, 59, 103–13.
[158] P. Lenahan, in Defects in Microelectronic Materials and Devices, ed. D. Fleet-

wood, R. Schrimpf and S. Pantelides, Taylor and Francis/CRC Press, 2008,
ch. 6, pp. 163–214.
[159] R. Van Ginhoven, H. Hjalmarson, A. Edwards and B. Tuttle, Nucl. Instr. &
Meth. Phys. Res. Sect. B, 2006, 250, 274–278.
[160] S. Pantelides, S. Rashkeev, R. Buczko, D. Fleetwood and R. Schrimpf, IEEE

Trans. Nucl. Sci., 2000, 47, 2262–2268.
[161] S. Rashkeev, D. Fleetwood, R. Schrimpf and S. Pantelides, IEEE Trans. Nucl.

Sci., 2001, 48, 2086 –2092.
[162] P. Bunson, M. D. Ventra, S. Pantelides, D. Fleetwood and R. Schrimpf, IEEE

Trans. Nucl. Sci., 2000, 47, 2289–2296.
[163] A. Alkauskas and A. Pasquarello, Phys. B Condens. Matter, 2007, 401-402,

546–549.
[164] J. Conley and P. Lenahan, IEEE Trans. Nucl. Sci., 1992, 39, 2186–2191.
[165] S. T. Pantelides, L. Tsetseris, S. N. Rashkeev, X. J. Zhou, D. M. Fleetwood and

R. D. Schrimpf, Microelectron. Reliab., 2007, 47, 903–911.
[166] J. K. Rudra and W. B. Fowler, Phys.Rev.B, 1987, 35, 8223–8230.
Page 168
[167] V. V. Afanas’ev and A. Stesmans, J. Phys.: Condens. Matter, 2000, 12, 2285.
[168] V. V. Afanas’ev and A. Stesmans, Europhys. Lett., 2001, 53, 233–239.
[169] K. Kajihara, L. Skuja, M. Hirano and H. Hosono, Phys. Rev. Lett., 2002, 89,
135507.
[170] L. Skuja, K. Kajihara, M. Hirano, A. Saitoh and H. Hosono, J. Non-cryst. Sol.,

2006, 352, 2297–2302.
[171] L. Skuja, K. Kajihara and H. Hirano, M. andHosono, Nucl. Instr. and Meth. in
Phys. Res. B, 2008, 266, 2971–2975.
[172] H. Li and J. Robertson, J. Appl. Phys., 2014, 115, 203708.
[173] C. Van de Walle and J. Neugebauer, Nature, 2003, 423, 626–628.
[174] J. Godet and A. Pasquarello, Microelectron. Engineering, 2005, 80, 288–291.
[175] A. H. Edwards, J. A. Pickard and J. E. Stahlbrush, J. Non-cryst. Solids, 1994,

179, 148–161.
[176] M. Vitiello, N. Lopez, F. Illas and G. Pacchioni, J. Phys. Cem. A, 2000, 104,
4674–4684.
[177] A. Yokozawa and Y. Miyamoto, Phys. Rev. B, 1997, 55, 13783–13788.
[178] Z. Li, S. J. Fonash, E. H. Poindexter, F. Harmatz, F. Rong and B. W. R., J.

Non-Cryst. Solids, 1990, 126, 173.
[179] J. Godet and A. Pasquarello, Phys. Rev. Lett., 2006, 97, 155901.
[180] P. Bunson, M. D. Ventra, S. Pantelides, R. Schrimpf and K. Galloway, IEEE

Trans. Nucl. Sci., 1999, 46, 1568–1573.
[181] H. A. Kurtz and S. P. Karna, J. Phys. Chem A, 2000, 104, 4780–4784.
[182] G. Vaccaro, S. Agnello, G. Buscarino, L. Nuccio and S. Grandi, J. Appl. Phys.,

2009, 105, 093514.
Page 169
[183] L. Skuja, K. Kajihara, M. Hirano and H. Hosono, Nucl. Instrum. Methods Phys.
Res., 2008, 266, 2971–2975.
[184] J. Li, S. Kannan, R. L. Lehman and G. H. S. Jr., Appl. Phys. Lett., 1995, 66,
2816.
[185] E. Cartier, J. H. Stathis and D. A. Buchanan, Appl. Phys. Lett., 1993, 63,
1510–1512.
[186] J. Suñé and E. Y. Wu, Phys. Rev. Lett., 2004, 92, 087601.
[187] P. Nicollian, A. Krishnan, C. Bowen, S. Chakravarthi, C. Chancellor and

R. Khamankar, IEEE International Electron Devices Meeting 2005, Techincal
Digest, 2005, pp. 403–406.
[188] E. H. Poindexter, J. Non-Cryst. Solids, 1995, 187, 257 – 263.
[189] Q. Fang, J.-Y. Zhang, Z. Wang, M. Modreanu, B. O’Sullivan, P. Hurley, T. Leed-

ham, D. Hywel, M. Audier, C. Jimenez, J.-P. Senateur and I. W. Boyd, Thin
Solid Films, 2004, 453454, 203 – 207.
[190] V. Singh, S. K. Sharma, D. Kumar and R. K. Nahar, Microelectron. Engineering,

2012, 91, 137 – 143.
[191] X.-Z. Li, B. Walker and A. Michaelides, Proc. Natl. Acad. Sci., 2011, 108,
6369–6373.
[192] S. Muller and T. I. Kamins, Device Electronics for Integrated Circuits, Wiley,
2003.
Page 170

Thesis

Uploaded by

Copyright:

Available Formats

Thesis

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis

Uploaded by

Copyright:

Available Formats

Atomistic Modelling of Charge Trapping Defects

Al-Moatasem Bellah El-Sayed

Department of Physics and Astronomy

University Collge London

A thesis submitted in partial fulfilment of the requirements for the

degree of Doctor of Philosophy (Ph.D.) of the University College

Supervisor: Prof. Alexander L. Shluger

2 Theoretical Background and Methods 11

3 Modelling Amorphous Silicon Dioxide and the Si/SiO2 interface 37

4 Defects in Bulk SiO2 66

6 Intrinsic Electron Trapping Sites in Silicon Dioxide 102

7 Optical Absorption Spectra of Trapped Electrons in SiO2 117

8 Hydrogen’s Interactions with Silicon Dioxide 129

9 Concluding Remarks 151

A BKS Parameterisation for Si and O 155

B Formation Energy Calculation 156

1.1 A transmission electron micrograph (TEM) image of a silicon MOSFET. 3

2.1 Silicon-silicon bond-order dependence on interatomic distance. . . . . 29

4.1 Crystal structure of α-quartz . . . . . . . . . . . . . . . . . . . . . . . 68

5.1 Ge impurity in α-quartz . . . . . . . . . . . . . . . . . . . . . . . . . 93

7.1 Experimental optical absorption spectra of electron-irradiated a-SiO2 . 119

3.1 Bond lengths of the Si3 Oy clusters . . . . . . . . . . . . . . . . . . . . 41

4.1 Comparison of neutral oxygen vacancy geometries in different models 69

5.1 Hyperfine splittings of the Li electron trap in α-quartz. . . . . . . . . 99

8.1 Hydrogen’s reactions with a-SiO2 . . . . . . . . . . . . . . . . . . . . . 148

1. A. -M. El-Sayed, Y. Wimmer, W. Goes, T. Grasser, A. L. Shluger, “Theoretical

2. A. -M. El-Sayed, K. Tanimura, A. L. Shluger “Electron Localization in Amor-

3. A. -M. El-Sayed, M. B. Watkins, T. Grasser, V. V. Afanasev, A. L. Shluger,

4. A. -M. El-Sayed, T. Grasser, V. V. Afanas’ev, A. L. Shluger, “Hydrogen Induced

5. A. -M. El-Sayed, M. B. Watkings, A. L. Shluger, V. V. Afanas’ev, “Nature of

6. A. -M. El-Sayed, M. B. Watkins, A. L. Shluger, V. V. Afanas’ev, “Identification

7. S. Ling, A. -M. El-Sayed, M. B. Watkins, V. V. Afanas’ev, A. L. Shluger, “A

8. A. -M. El-Sayed, A. L. Shluger, “Atomistic Modelling of Defects Implicated in

I would like to take this opportunity to offer my heartfelt appreciation to a number of

1.2 Electronic Device Reliability Issues

Figure 1.1: A transmission electron micrograph (TEM) image of a silicon MOSFET. 2

Figure 1.3: a) Schematic band diagram of tunnelling of electrons through defects in

cause of the negative bias temperature instability (NBTI). 9 NBTI is characterised by

1.3 Atomistic Nature of Defects Implicated in Re-

To aid in the understanding of reliability issues, computational methods are used

1.4 Scope of this Work

This thesis presents atomistic calculations of a number of defects in silicon dioxide

an energetic barrier. The character of localization at intrinsic sites was shown to

2.2 Molecular Dynamics

Molecular dynamics is a computer simulation method in which the force on each

Adding equations (2.3) and (2.4) gives:

r(t + δt) = 2r(t) − r(t − δt) + δt2 a(t). (2.5)

r(t + δt) − r(t)

2.3 Density Functional Theory

DFT is a theory providing information on the electronic properties and chemical

E[ρ(r)] = E T [ρ(r)] + E V [ρ(r)] + E H [ρ(r)] + EXC [ρ(r)], (2.7)

a landmark paper, Kohn and Sham introduced an auxiliary system of non-interacting

ψ1 (r1 ) ψ2 (r1 ) . . . ψN (r1 )

1 ψ1 (r2 ) ψ2 (r2 ) . . . ψN (r2 )

Egs = hΨgs [ρ]|Ĥ|Ψgs [ρ]i, (2.11)

vef f = vext (t) + vHartree (t) + vxc (t), (2.17)

where ρ(t) is the time-dependent density. By self-consistently calculating a stationary

Aia,jb = δij (a − b ) + (ij|jb) + (ia|fxc |jb), (2.20)

2.3.1 Gaussian Plane Waves Method

Gaussian Basis Sets

Aia,jb = δij (a − b ) + (ij|jb) + (ia|fxc |jb), (2.20)