The search for the detailed kinetics responsible for the creation of silicon dangling bonds at the Si/SiO\(_2\) interface has proven to be a challenging task with no consensus being reached within the past decades. While
experimental investigations for thermal dissociation in conjunction with ESR studies by Brower [73] and Stesmans [201] converged to an activation energy between \(\SI {2.5}{}\) and \(\SI {2.8}{eV}\), the theoretical
picture behind this phenomenon is not well understood. Pioneering theoretical calculations have been conducted by Tuttleet al. [202–206] and Van de Walleet
al. [207–209] investigating the structural and energetic properties of various
configurations of hydrogen in periodic models of crystalline and amorphous Si. Both authors concluded that H dissociation along the Si–H bond direction (stretching mode) into SiO\(_2\) is not the main dissociation path, due to
its high reaction barrier of \(\SI {3.6}{eV}\) , basically the Si–H binding energy7. The authors proposed that the Si–Si bond–center (BC) in bulk Si provides a stable configuration for the
hydrogen atom, particularly in its neutral charge state [48, 210, 211], which is confirmed by experimental data [212–214]. By investigating different configurations it was shown that placing H between strained Si–Si bonds is between \(\SI {1.7}{eV}\) and \(\SI {2.5}{eV}\)
higher in energy than the Si–H equilibrium configuration, depending on the actual position of the BC site. However, it was argued that the reaction barrier would be essentially the same as the difference in energy, indicating that
the H would spontaneously reform the Si–H bond as there is no reverse barrier. In [203] it was then speculated that by bending the H across the next
BC site into the Si–H antibonding configuration (AB), which is a \(\SI {180}{\degree }\) flipped position of the H, emerging energy levels in the Si band gap may aid the dissociation process by allowing the Si–H complex to
become charged. However, as we have shown recently [MJJ3], the AB site is actually only a metastable configuration without active trapping levels
in bulk Si and as such a dead–end in the reaction dynamics. Instead, a new dissociation trajectory was proposed where in an initial step the hydrogen bends towards an adjacent silicon atom, and subsequently relaxes into the
next but one Si–Si bond on the silicon side of the Si/SiO\(_2\) interface. The details of these calculations and the derived conclusions will be discussed in the following.
7 Assuming a void above the Si–H bond which renders it compatible to dissociation into free space.
3.2.1 Reaction Paths via Metadynamics
Finding the reaction dynamics and the corresponding pathway, i.e. locating a saddle point on a multidimensional PES \(V(\mathbf {q})\), is a highly challenging issue in modern computational chemistry methods. Among the
various methods, the NEB method [215, 216] is a popular way to find reaction
pathways and minimum energy paths (MEPs) when both, the initial and final states, are known. However, for the problem at hand, the reaction kinetics for Si–H, only the equilibrium position, i.e. the initial
state, is known, which renders the NEB method inappropriate. A complimentary method to NEB is the so–called dimer method [217], which only
requires an initial state as input. Without any assumption regarding the unstable normal mode associated with the transition state, the initial dimer vector is randomly chosen. Test calculations, however, showed that this approach
fails to converge to any reasonable results for finding the transition state for a single bond configuration in a solid–state system. A more rigorous approach would require a reasonable guess for the transition state together with the
respective imaginary phonon mode, which, however, is the great unknown and a major subject of the current investigations.
Figure 3.7: Representation of the free energy landscape \(V(r,\phi ,\theta )\) in 3D space obtained from the well–tempered metadynamics (WTMD) simulations parametrized with respect to the
three collective variables. Left: Isosurface for \(E=\SI {0.8}{eV}\); one can see the vicinity of the Si–H equilibrium positions and the AB site. Middle: Isosurface for \(E=\SI {1.8}{eV}\); the trajectory between
the initial and AB configuration is already explored. Furthermore, the newly found next but one BC site is visible. Right: Isosurface for \(E=\SI {2.4}{eV}\); an emerging transition path between the equilibrium
and the next but one BC configuration.
Figure 3.8: The free energy landscape together with a schematic of the different atomistic configurations in the vicinity of the interfacial Si–H bond obtained from WTMD with the ReaxFF
force field. Three distinct minima are visible: The equilibrium position 1 at \(r=\SI {1.37}{\angstrom }\) and \(\phi =\SI {0}{\degree }\), the hydrogen atom in the AB site labeled with 3, as well as the BC\(_{2,3}\) configuration 5 where the H moved in between the next but one Si\(_2\)–Si\(_3\) bond. The crosses 2 and 4 mark the transition barriers for the extracted minimum energy paths (dashed and solid line) connecting the different minima.
A very recently developed enhanced–sampling method is WTMD [218, 219].
This method allows one to efficiently explore the free energy landscape (FES) by adding an artificial Gaussian potential to the real energy landscape which continuously expands the
sampling of the FES and simultaneously discourages the system of returning to previously sampled configurations. Biasing certain collective variables (CVs), therefore, drives the Si–H complex out of its equilibrium and triggers
possible reactions along the stretching and bending mode, while at the same time the environment can be kept at a constant temperature, e.g. \(T=\SI {300}{K}\). Thus, reaction barriers can be overcome which are not accessible
in conventional molecular dynamics (MD) simulations. Having obtained the free energy surface, parametrized in terms of the CVs, the MEP connecting two configurations can be extracted.
The simulations conducted here are performed using the rather large Si/SiO\(_2\) model system (\(\sim 500\) atoms) described above. Together with the considerable amount of simulation steps necessary to converge such
simulations, only classical force–fields, such as ReaxFF [197], can be applied. All calculations again use the MD
engine Lammps [198] in conjunction with the library Plumed [220]. The starting point is again an equilibration phase for \(\SI {100}{ps}\) at \(T=\SI {300}{K}\) by assigning random, normally distributed velocities to
the whole system. Subsequently, several WTMD simulations with different bias factors and bias heights have been conducted for a maximum of \(\SI {50e6}{}\) timesteps with a stepsize of \(\SI {0.5}{fs}\). Detailed
informations about the convergence is given in Appendix F. Within the simulations two CVs were biased and monitored to describe the
system. The Si–H bond distance \(r\) as well as the polar angle \(\phi \) with respect to its equilibrium position. The additional azimuthal angle, however, was only monitored and not explicitly biased during the WTMD
simulations. The accessible region was limited to within \(\SI {4}{\angstrom }\) to ensure sampling of the free energy landscape only in the direct vicinity of the Si–H bond’s equilibrium position.
The raw data, \(V(r,\phi ,\theta )\), are shown in 3D space as isosurfaces for selected values of the potential energy, see Fig 3.7. In order
to give a more intuitive picture of the free energy landscape, the potential was mapped into 2D space by only considering the CV \(r\) and the polar angle \(\phi \). Three distinct minima can be identified: 1 the intact Si–H bond in its equilibrium configuration, 3 corresponds to the H being in the AB site and 5 represents a newly identified minimum formed by the H being in the next but one BC configuration between Si\(_2\) and Si\(_3\). The respective atomistic configurations are schematically illustrated in
Fig. 3.8 together with the extracted FES. The marked minimum energy paths connecting the three configurations have been verified in separate
subsequent simulations using more efficient CVs, such as \(r^\prime =1/\sqrt {2}(r_\mathrm {Si_3,H}-r_\mathrm {Si_1,H})\) and simply the angle \(\phi \).
The free energy map shows that inverting the H around Si\(_1\) into the AB site 3, which is \(\SI {0.8}{eV}\) higher in energy, proceeds via the dashed path. The transition barrier is \(\SI {1.75}{eV}\) and its configuration corresponds to the H atom in the BC\(_\mathrm {1,2}\) site 2. Interestingly, although the H forms a bond–center configuration between Si\(_1\) and Si\(_2\), it is not a preferred metastable configuration for H along the given path, contrary to it being a stable site for H in bulk
\(c\)–Si. In addition to the configuration 3 the hydrogen can also be moved into the BC\(_{2,3}\) site 5, shown as the solid black line. The extracted MEP yields a transition state 4 where the H is stretched away from its initial silicon, Si\(_1\), and attached to the adjacent Si\(_2\). With a forward barrier of \(\SI {2.25}{eV}\) and a reverse activation energy of \(\SI {1.05}{eV}\) this
trajectory could potentially explain the measurement data by Stesmans [201], reporting a barrier height of \(\SI
{2.83}{eV}\) for the Si–H depassivation process. The reverse passivation barrier \(\Rightarrow \)Si–H, on the other hand, is assumed to be mainly driven by the cracking of molecular H\(_2\), +H\(_2\Rightarrow \)Si–H+, where the additional (neutral) atomic H potentially is released into the SiO\(_2\) side. Therefore, the backward barrier can not be directly compared to the measured values in [221, 222]. Nevertheless, a barrier of around \(\SI {1}{eV}\) prevents the H from
directly going back to the created Si–DB and potentially enables new insight into device degradation due to interface defects. A further discussion will be given in Section 3.6.
3.2.2 Detailed Dynamics via DFT
Building on top of the results derived above, the Si–H kinetics have been further investigated using density functional theory (DFT). The initial 1 and final 5 configurations obtained from the classical force–field calculations were used to construct a trajectory with 13 frames in total by using linear interpolation8. Subsequently, CI–NEB simulations optimized the
whole band, including its endpoints. Fig. 3.9 shows the resulting energy along the converged trajectory, i.e. the one–dimensional PEC for the introduced
RC which represents a collective motion of atoms. The resulting path possesses very similar features as already predicted by ReaxFF observed in Fig. 3.8. First, the hydrogen moves towards Si\(_2\), marking the transition state, and eventually moves in between the Si\(_2\)–Si\(_3\) bond which was stretched
from \(\SI {2.34}{\angstrom }\) to \(\SI {3.13}{\angstrom }\), thereby forming a BC configuration. The total reaction barrier along the trajectory separating the intact Si–H configuration and the BC\(_{2,3}\) site is
\(\SI {2.77}{eV}\) for the forward reaction and \(\SI {1.30}{eV}\) for repassivating the Si dangling bond. Compared to the calculations with the classical force–field ReaxFF, the activation
energy for creating a Si–DB is even closer to the experimentally extracted barrier of \(\SI {2.83}{eV}\) [201]. Both, the transition state 4 as well as the final state 5, introduce two localized electronic levels in the Si band gap; a filled state close to the valence band edge and an empty level in the upper half of the band gap, see Fig. 3.9. Furthermore, the molecular orbitals (MOs) associated with these states, the highest (lowest) occupied (unoccupied) MO, are fully localized around the
unpassivated Si, see Fig 3.10, indicating the creation of an electrically active interface defect [88, 93, MJC3]. In contrast the next lower (higher) MOs clearly show the delocalized characteristics expected for band states. Compared to the Si–DB characteristics
shown above, the spin density, however, suggests a slightly distorted \(sp^3\) hybridized DB orbital with a significant spread onto the back–bonded Si atoms.
Figure 3.9: The results of a climbing–image nudged elastic band (CI–NEB) simulation investigating the direct MEP connecting the intact Si–H configuration 1 and the next but one BC site between Si\(_2\) and Si\(_3\) 5 via the transition state 4 where the H is attached to the adjacent Si\(_2\). Top: Simulation snapshots showing the atomistic configurations of the initial, final and transition state. Left: The total energy required for the forward
reaction is \(\SI {2.77}{eV}\), whereas the backward barrier is \(\SI {1.30}{eV}\) for moving the H back to its initial configuration. Right: Difference of the projected DOS onto the Si\(_1\) atom along the path
compared to the intact Si–H configuration. It is clearly visible that the transition and also the final state introduce localized electronic states in the Si band gap.
Figure 3.10: Left: The HOMO-1 and LUMO+1 are delocalized Si band states. However, they partially penetrate into the oxide region as was already discussed in Sec. 3.1.
Right: The HOMO and the lowest unoccupied molecular orbital (LUMO) are clearly localized around the Si–DB. Comparing the wavefunctions to the prototypical defect configuration in Fig. 3.5 suggests the dangling bond orbital to be constituted by a slightly distorted unpaired Si \(sp^3\) hybrid.
Additionally, the charges associated with the atoms along the bond breakage path and the various atomic configurations have been analyzed using Mulliken and Bader charge analysis [223–225]. Both methods show that the hydrogen actually dissociates in its neutral charge
states, leaving one remaining electron on the created Si–DB, in accordance with the occupation of the MOs. Details of the Bader charge analysis can be seen in Fig. 3.11, which shows the change of the electronic density at the transition and final state. Integrating over the associated Bader volumes of Si\(_1\) shows
that the respective charge changes by \(0.78\,e\) as the H migrates to Si\(_2\), and by \(0.95\,e\) for the H being in its final position. This clearly indicates that one electron remains on the Si–DB. Analyzing the H along the
trajectory shows that its charge changes by \(-0.18\,e\), becoming slightly negatively charged in the final BC configuration. However, in summary one can conclude that the H would dissociate in its neutral charge state. Another
important detail can be extracted using the Bader charge analysis. Along the path the H first attaches to Si\(_2\) until it relaxes into the final position between Si\(_2\) and Si\(_3\). Note that such a configuration was already
confirmed in the literature to be a stable position for H\(^0\) [210–214]. A
closer look reveals that the H would actually be bound to Si\(_3\), as suggested by the electronic density as well as the distance between the H and the respective Si atoms. While at the transition point the H–Si\(_2\) distance is
\(\SI {1.65}{\angstrom }\), it increases to \(\SI {1.68}{\angstrom }\) for the final position, compared to \(\SI {1.55}{\angstrom }\) between the H and Si\(_3\). Thus, the Si\(_2\) also possesses the character of a
Si–DB with one unpaired electron. However, no additional states in the Si band gap are created due to the interaction with the nearby H atom which moves its energy levels inside the valence and conduction band.
Figure 3.11: The charge density is represented by the blue and red translucent profiles along the dissociation trajectory (blue: increase of the electron density, red: decrease of the electron density w.r.t. the initial configuration).
Left: The equilibrium Si–H configuration and the corresponding distribution of charge. Middle: In the initial step the H attaches to an adjacent Si atom (Si\(_2\)) with a distance of \(\SI {1.57}{\angstrom }\).
Thereby, one electron remains on the initial Si\(_1\) as indicated by the Bader charge analysis. Right: In the final configuration the H forms a Si–H–Si complex between the next but one Si–Si bond. Being in this position,
the H is slightly negatively charged and almost fully bonded to Si\(_3\) as suggested by the charge density difference as well as the distances to Si\(_2\) and Si\(_3\).
8 Note that this corresponds to a direct connection between the minima in Fig. 3.8.
3.2.3 Statistical Analysis
The inherent structural disorder at the Si/SiO\(_2\) interface gives rise to a distribution of Si–Si and Si–O bond lengths and angles. Linking theoretical data and experimental results is, therefore, only possible at the statistical
level. To expand the simulation results, three different Si/\(a\)–SiO\(_2\) models have been used with a total number of 13 variations of pristine Si–H bond configurations, consistent with the perceptions of a \(P_\mathrm {b}\)
center. The different Si–H bonds were created by breaking and passivating selected Si–Si or Si–O bonds at the interface. Thus, the statistics include \(P_\mathrm {b0}\) (back–bonded to three Si’s) and \(P_\mathrm {b1}\)
(back–bonded to two Si’s and one O)–like types as well as configurations with a nearby H atom and hydroxyl groups, respectively, see Fig. 3.5. In order to facilitate comparison with experimental results, all initial Si–H configurations were placed within \(\SI
{4}{\angstrom }\) of the subinterfacial Si side. For all 13 starting configurations the three next but one BC sites were chosen and the direct trajectories were calculated using the CI–NEB method. Overall, this procedure gave 38
dissociation paths together with their respective energetics9. Owing to computational limitations, the CI–NEB simulations have been performed using the Pbe functional.
Subsequently, single point calculations based on the hybrid Pbe0 functional were carried out on the optimized structures along the trajectory.
Figure 3.12: Statistical analysis of 38 calculated dissociation trajectories and barriers for various initial Si–H configurations using the CI–NEB method. Left: While the colored bars show the full results using all defect
configurations, the gray bars present a limited statistics where the initial Si–H bonds have been selected to be at least one layer away from the interface. Right: The extracted potential energy curves together with the mean
PEC. One can clearly see the variation of the reaction coordinate (RC) due to the amorphous interface region.
The simulation results show a broad distribution of forward (backward) barriers of the investigated reaction, see Fig. 3.12, ranging from \(\SI
{2.07}{}\) (\(\SI {0.95}{}\)) to \(\SI {2.95}{eV}\) (\(\SI {1.94}{eV}\)) with a standard deviation of \(\SI {0.20}{}\) and \(\SI {0.23}{eV}\), respectively. The mean value for breaking the bond averages at \(\SI
{2.57}{eV}\) and nicely matches the experimental value of \(\SI {2.83}{eV}\); however, the extracted standard deviation \(\sigma \) is actually much larger than the value reported by Stesmans and Brower (\(\sigma =\SI {0.06}{eV}\)). Two main factors possibly influence the results leading to the observed discrepancy. First, the created and utilized atomistic
models possibly introduce an artificial local strain at the transition region due to their finite dimensions. This would result in artificially distorted bond angles and lengths, see Appendix A, which broadens the calculated statistics. However, further increasing the model sizes challenges the computational limits and prohibits accurate calculations using hybrid
functionals. Secondly, the spatial placement of the initial Si–H bond could bias the resulting statistics. Within the presented statistics the pristine configurations are also placed directly at the Si/SiO\(_2\) interface, where the
atoms experience the maximum distortion due to the amorphous oxide, see Sec. 3.1 and Fig. 3.2. Limiting the statistics to initial configurations residing in a subinterfacial Si environment, at least one layer away from the
interface, results in a much narrower distribution, see Fig 3.12. The concept of interface defects not being directly at the transition region, was
already proposed in recent publications [92, 94], motivated by the small
\(\sigma \) measured for the activation as well as the narrow ESR spectra, and thereby is consistent this assumption.
9 One calculation did not converge to the defined criteria and was therefore discarded.
\(\Rightarrow \)Si–H, on the other hand, is assumed to be mainly driven by the cracking of molecular H\(_2\), 
, where the additional (neutral) atomic H potentially is released into the SiO\(_2\) side. Therefore, the backward barrier can not be directly compared to the measured values in [221, 222]. Nevertheless, a barrier of around \(\SI {1}{eV}\) prevents the H from
directly going back to the created Si–DB and potentially enables new insight into device degradation due to interface defects. A further discussion will be given in Section 3.6.
