In order to calculate the contributions of each mechanism to HCD described Chapter 1, the Si–H bond and its characteristics need to be explored and analyzed by means of ab initio methods. The bond breakage trajectory of how H would eventually dissociate away from Si at the Si/SiO\(_2\) environment is largely unknown and no consensus has been reached on a microscopic level. Furthermore, the theoretical treatment of resonance induced excitations requires the knowledge of the associated excited states. An electron scattering into a resonance forms a temporary anionic system1, which evolves within its lifetime on the excited potential energy profile. Therefore, a major ingredient in the above proposed descriptions are the energetic positions of the resonances, their lifetimes and the excited potential energy surfaces (PESs). Additionally, the Si-H bonds’ response upon an applied electric field, i.e. its dipole moment, needs to be investigated. All of the aforementioned interactions trigger vibrational excitations and eventually dissociation of the Si–H bond. Therefore, the vibrational coupling with its environment, described as a phonon bath of oscillator, also needs to be accounted for. The coupling strength determines the vibrational relaxation, and thus the lifetime of vibrationally excited states, which influences the dynamics of the bond.
1 Conversely, a hole temporarily forms an excited cationic system.
The following part of this thesis focuses on the theoretical description of Si–H bonds and Si dangling bonds (DBs) at the Si/SiO\(_2\) interface. In order to be confident that the results are truly representative of the properties of Si/SiO\(_2\) structures and the corresponding defects, the initial atomistic models must be consistent with experimental perceptions of the atomic configurations and characteristics. This Section discusses the techniques used to create realistic Si/SiO\(_2\) models and the methods applied to subsequently model the behavior and properties of interfacial Si–H bonds.
Models of Si/SiO\(_2\)/Si interfaces were created using classical potentials by adapting the well known melt–quench technique applied to construct amorphous bulk models, similar to recent approaches reported in the literature [176–181]. The obtained models were further optimized using the DFT setup described below. The electronic structure as well as geometric properties have been extracted and validated against experimental data from various results reported in the literature. Additionally, the Si dangling bond defects, which natively formed during the melt and quench procedure, have been carefully investigated and compared against the stringent limitations of \(P_\mathrm {b}\) center characteristics found experimentally. Their structural properties as well as their behavior upon charge capture agree well with the measured specifications.
All simulations were carried out using the Cp2k package, a DFT code which uses a mixed Gaussian and plane waves approach to represent the electrons in the system. 3D periodic boundary conditions were used for the Si/SiO\(_2\) system throughout this work. A double–\(\zeta \) Gaussian basis set optimized for condensed phase systems for Si, O and H was employed in conjunction with the Goedecker–Teter–Hutter (GTH) pseudopotentials. The plane wave cutoff used in these calculations was set to \(\SI {650}{Ry}\)2. In order to minimize errors related to energies, i.e. band gaps and barrier calculations, the non–local, hybrid functional Pbe0_Tc_LRC was used for most of the work. It contains a \(20\%\) contribution of Hartree–Fock exchange and correlation in addition to a truncated Coulomb operator with a cutoff radius set to \(\SI {2}{\angstrom }\). To mitigate the computational cost of calculating the Hartree–Fock integrals within the non–local functional, the auxiliary density matrix method (ADMM) was employed. In addition to the main basis set, a more sparse auxiliary basis set (pFIT basis set for Si, O and H) was used to calculate the exchange terms, which allows to further reduce the computational expenses. All geometry and cell optimizations were performed using the BFGS algorithm to minimize forces on atoms and the cell to within \(\SI {1.5e-2}{eV/\angstrom }\) (\(\SI {5e-3}{GPa}\)). Energy barriers between different atomistic configurations were calculated using the climbing–image nudged elastic band (CI–NEB) method. The algorithm linearly interpolates between the initial and final configuration to generate a trajectory which is connected by springs. Subsequently, the whole band is then simultaneously relaxed with the spring’s force constant set to \(\SI {19.5}{eV/\angstrom ^2}\).
2 The cutoff was increased to \(\SI {1200}{Ry}\) for calculations including an electric field.
Constructing credible interface structures between amorphous oxides and a crystalline substrate has proven to be an extremely challenging problem [176–181], mainly due to the lack of detailed information regarding the microscopic interfacial region which, however, also depends on the exact process conditions. It is known from measurements performed in thin SiO\(_2\) films that thermal oxidation of silicon results in an intrinsic compressive stress within the oxide [182–186]. Depending on the oxide thickness and the processing temperature, the reported values vary between \(\SI {0.5}{}\) and \(\SI {2.0}{GPa}\)3. X–ray reflectivity experiments [187] confirm this trend and, furthermore, found an increased SiO\(_2\) density in an Si/SiO\(_2\) system compared to bulk amorphous silica together with a compresed Si–O–Si angle distribution [188, 189]. Moreover, transmission electron microscopy (TEM) images [190, 191] and atomic–scale electron–energy–loss spectroscopy (EELS) [192, 193] indicate a transition region of approximately \(\SI {5}{}-\SI {7}{\angstrom }\) between crystalline silicon and the amorphous oxide structure. Within this interfacial region the band gap as well as the dielectric constant continuously change between Si and SiO\(_2\), as was also proposed in recent DFT studies [194, 195]. Together with photoemission measurements [190, 196], which investigated the chemical structure of ultrathin interfaces, the data suggest three transition layers with increasing concentration of partially oxidized silicon. In summary, all experimental investigations provide a number of constraints an interface model must satisfy to reliably reflect realistic conditions.
The utilized methodology to create a 3D periodic Si/SiO\(_2\)/Si model with two interface regions adapts the melt and quench procedure previously used to create \(a\)–SiO\(_2\) cells [62, 174] and is schematically shown in Fig. 3.1. Due to the excessive number of calculation steps required within the molecular dynamics (MD) simulations, the ReaxFF force field [197] implemented in the Lammps code [198] was used. ReaxFF is a sophisticated potential, describing covalent interactions in terms of a bond order which allows for the dynamic formation and annihilation of bonds. In a recent study the applicability of ReaxFF to model glassy silica has been investigated [199]. The authors concluded that ReaxFF is a promising potential for creating and modeling amorphous systems which would even allow to study the surface reactivity thanks to its dynamical ability. Starting from a \(3\times 3 \times 3\) \(\beta \)–cristobalite cell, whose lateral cell dimensions were fixed to that of silicon, the energy has been minimized by allowing the cell to relax in the \(c\) direction. It was found that this procedure gives an oxide density of \(\sim \SI {2.3}{g/cm^3}\), in agreement with measured values [190]4. Subsequently, the cell dimensions were fixed during the melt and quench procedure. Crystalline silicon was added at the top and bottom and the respective atoms on either side were frozen while the remaining atoms were given random velocities from a Gaussian distribution. Furthermore, reflective walls have been introduced along both interfaces which prevent the oxygen atoms to penetrate into the \(c\)–Si lattice. In a first step the system was equilibrated at \(\SI {300}{K}\) for \(\SI {10}{ps}\) using a Berendsen thermostat. Afterwards the temperature was linearly increased to \(\SI {5000}{K}\) over \(\SI {60}{ps}\) to melt SiO\(_2\) and the system was further equilibrated for an additional \(\SI {100}{ps}\)5. Finally, the system was quenched to \(\SI {0}{K}\) at a rate of \(\SI {1.6}{K/ps}\) and the whole structure, including all atoms, was allowed to relax at \(\SI {300}{K}\). Selected snapshots of the different phases are shown in Fig. 3.1. In total over 100 Si/SiO\(_2\) models were created, which resulted in various defect configurations at the interface as well as in the oxide. The three most promising models, in terms of their geometrical properties such as coordination number, bond angles and lengths, were chosen and further examined as well as optimized using the DFT setup. Note that all undercoordinated atoms, at the interface as well as within the oxide, have been passivated with H for subsequent calculations.
3 A recent study which investigates low–temperature processing of SiO\(_2\) yields comparable results. A lower processing temperature effectively increases the interfacial strain and hence the defect density [MJC1].
4 The density of \(\beta \)–cristobalite is closer to the one of \(a\)–SiO\(_2\) compared to crystalline \(\alpha \)–quartz.
5 Within the simulations the melting point is at \(\sim \SI {4500}{K}\) and by \(\SI {5000}{K}\) the system is in its liquid phase.
Simulations of the nonoptimized interfaces, i.e. structures obtained from ReaxFF, indeed gave rise to a rather high intrinsic compressive stress of \(~\sim \SI {3.5}{GPa}\). Thus, in order to obtain a more realistic model structure which is compatible with experimental values mentioned above for the subsequent simulations and investigations, full cell optimizations including the ionic positions of the models have been performed with an external pressure of \(\SI {1}{GPa}\) parallel to the interface. Since the initial simulations did not show large off diagonal elements of the stress tensor, i.e. a shear stress, the original tetragonal cell symmetry was kept fixed during the optimizations. All three structures converged to a very similar cell size, \(a=b=\SI {16.203}{\angstrom }\), \(c=\SI {32.965}{\angstrom }\), with a variance of \(\SI {0.007}{\angstrom }\) in the lateral dimensions and \(\SI {0.014}{\angstrom }\) in the \(c\) direction.
The extracted geometrical properties of the final interface models are in good agreement with experimental data. The Si–O bond length ranges from \(\SI {1.46}{}\) to \(\SI {1.84}{\angstrom }\) and peaks at around \(\SI {1.64}{\angstrom }\), which is consistent with infrared absorption measurements [190]. The Si–O–Si angle averages at \(\SI {133.5}{\degree }\), which agrees well with the conclusions reached in [189], that the angle is reduced to \(\SI {135}{\degree }\) compared to \(\SI {148}{\degree }\) in amorphous bulk structures. On the other hand, the O–Si–O angles remain almost unchanged compared to bulk SiO\(_2\), around \(\SI {109}{\degree }\). In order to further quantify the interface quality, the deformation of the \(c\)–Si lattice was calculated. This deformation is defined as the deviation of the atomic positions from their respective positions in crystalline silicon, see Fig. 3.2 (middle panel). Already three layers away from the Si/SiO\(_2\) interfacial region the crystalline structure is almost completely restored, with displacements less than \(\SI {0.2}{\angstrom }\). However, directly at the transition region, distortions are stronger due to the formation of Si–O bonds, giving rise to a higher strain directly at the interface between Si and SiO\(_2\). A detailed discussion on the geometrical properties as well as the influence of the cell size and the oxide thickness is given in Appendix A.
In addition to the structural quality, the electronic structure of the interface systems have also been investigated. Mulliken population analysis was used to determine the oxidation state of silicon across the Si/SiO\(_2\) model, see Fig. 3.3. The profile clearly shows a gradual change from Si\(^0\) in the crystalline Si part to its fully oxidized form Si\(^{+4}\) (which corresponds to a net charge of \(\sim 1.4e\) in our DFT setup) in the SiO\(_2\) region over a \(\sim \SI {5}{\angstrom }\) transition region. Such a characteristic is qualitatively very similar to the measurement results in [196]. Furthermore, this transition region is also reflected in the calculated DOS across the model. While in the bulk regions the corresponding band gaps compare well to the known values, one can see a continuous change within the interfacial regions, see Fig. 3.2 (left and right panels). However, additional electronic states within the band gap for the interface as well as in bulk SiO\(_2\) are visible. These induced states can be identified as silicon band wave functions extending into the near oxide region, see Fig. 3.4. Such features have been also reported in recent experimental studies [192, 193], where the authors concluded that this actually limits the practical oxide thickness to \(\SI {1.2}{nm}\). It may be inferred that within the models used here, having an oxide thickness of \(\SI {1.1}{nm}\), even region 3 does not exhibit actual full bulk SiO\(_2\) properties.
The preselected (using ReaxFF) and optimized (using DFT) interface models possess three different silicon dangling bond configurations (which have been passivated by H) at the interface as shown in Fig. 3.5: a) a Si–DB back–bonded to two Si’s and one O atom, b) a Si–DB directly at the interface with an hydroxyl group \(\SI {1.9}{\angstrom }\) away and c) a subinterfacial (one layer away from the interface) Si trivalently back–bonded to three other Si atoms with no other O or H atoms in the direct vicinity. All defect structures possess an unpaired electron, as indicated by the spin density, with no other defect being present in the model. Compared to ESR studies [83, 88, 91–93], which suggest that a \(P_\mathrm {b}\) center consists of a Si vacancy bonded to three Si atoms with an \(sp^3\) hybridized dangling bond orbital, configuration (c) fully meets these requirements, and is therefore chosen as a prototypical defect within this study. Subsequently, the different charge states, neutral, negative and positive, of different defect configurations have also been calculated. The structural reconfigurations with respect to their neutral configuration are shown in Fig 3.6. One can see that due to the rigid crystalline silicon structure the Si–DB does not undergo large atomic relaxations, neither in the negative (blue) nor positive (red) charge state. However, the flexibility of the \(a\)–SiO\(_2\) network allows close–by oxygen atoms to reconfigure according to the Si–DB charge state, particularly when positively charged. This, however, is to be expected, since oxygen is slightly negatively charged in the SiO\(_2\) system. In order to reliably extract the formation energies and the respective charge transition levels, the electrostatic correction methodology presented in [200] and implemented in sxdefectalignment2d has been applied to calculate the correction energy. All three defects exhibit two charge transition levels (\(+/0\) and \(0/-\)) within the Si band gap, see Fig 3.6, as is also shown in [84–86, 88, 93, MJC3]. Even the defect structure (c) in Fig. 3.6, although potentially resembling a \(P_\mathrm {b1}\) configuration, exhibits two transition levels within the Si band gap6. A detailed defect analysis can be found in Appendix B.
