The a-SiO₂/Si (100) and (111) interface models used in this work are derived from our previous works (Li et al., 2019a; Li et al., 2019b; Hong and Zuo, 2020). From a-SiO₂ and crystalline Si, the models are constructed through thermally oxidization of Si by a-SiO₂ using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code (Plimpton, 1995; van Duin et al., 2001; van Duin et al., 2003; Fogarty et al., 2010). Passivation can improve the stability of the structure and properties, so it has important scientific significance and broad application prospects (Ya-Hui et al., 2020). In this manuscript, hydrogen is used to passivate dangling bond defects produced during interface thermal oxidation. Dimensions of the final a-SiO₂/Si (100) and (111) interface models are 15.5 Å × 15.5 Å × 35.0 Å and 13.3 Å × 15.4 Å × 49.4 Å, respectively. In the a-SiO₂/Si (100) interface model, the thicknesses of the seven-layer Si (100) moiety and a-SiO₂ moiety are 8.1 Å and 12.7 Å, respectively, with a vacuum layer of 11.8 Å and a hydrogen passivation layer of 2.4 Å. In the a-SiO₂/Si (111) interface model, the thicknesses of the eight-layer Si (111) moiety and a-SiO₂ moiety are 10.2 Å and 15.1 Å, respectively, with a vacuum layer of 20.9 Å and a hydrogen passivation layer of 3.2 Å.

The first-principle calculations based on the Density Functional Theory are performed with the plane-wave basis, as implemented in the Vienna ab initio simulation package (VASP) (Kresse and Furthmüller, 1996). The electron core interaction is included by using the frozen-core all-electron projector augmented wave (PAW) method (Blöuml, 1994). The exchange correlation functional is treated by the generalized gradient approximation (GGA) parameterized by PERDEW, BURKE and ERNZERHOF (PBE) (Perdew et al., 1996). The plane-wave basis set cut-off energy is set as 520 eV for all calculations. Since the a-SiO₂/Si models are large in size, the Brillouin-zone integration is limited to the Γ point only. The structure is relaxed until the total energy is minimized. The convergence criterion of structural optimization is that the total energy difference is less than 10⁻⁵ eV, and the force on each atom is smaller than a tolerance of 0.05 eV/Å.

Harmonic Transition State Theory (HTST) (Vineyard, 1957) is derived from a well-known theory—Transition State Theory (TST) (Henry, 1935; Voter and Doll, 1984). For many solid-state reactions, HTST is a good approximation as the energy barrier is usually larger than the average kinetic energy k_B T in system. Another reason for using HTST as a reasonable approximation is that the harmonic approximation is often effective. Near the saddle point, the potential energy increases rapidly and forms a small bottleneck area. When using HTST, the first step is to find the saddle point corresponding to the reaction. There are many ways to find saddle points, and in this study we will use Climbing Image–Nudged Elastic Band (CI-NEB) method.

G. Henkelman developed the CI-NEB method to search for a transition state which is the first-order saddle point (Henkelman et al., 2000). The CI-NEB method specifically considers the problem of locating the transition state. The key difference between CI-NEB and Nudged Elastic Band (NEB) is the definition of the force on the point with the highest energy. In CI-NEB, this point will not be subjected to the spring force of the adjacent points, which can prevent its position from being pulled away from the transition state, therefore making this point climb to the transition state along the path in the direction of energy increase. The CI-NEB method requires only a few points, such as a total of five or even three points including the initial and final states to accurately locate the transition state. It is one of the most efficient methods to find a transition state. Between the initial and final structures of a reaction, appropriate number of images are inserted. In this work, the transition state and activation energy of the P _b-type defect depassivation reaction are obtained by the CI-NEB method.

In order to simulate the stress induced by the final passivation layer, a biaxial strain is applied at the a-SiO₂/Si (111) and (100) interfaces (Zhou et al., 2014; Grote and Berger, 2015; Al-Shami et al., 2018). The biaxial strain is imposed on the relaxed unit cell along x-axis and y-axis ([100] and [010] direction) according to the equations (Al-Shami et al., 2018): ε a ( % ) = ( a − a 0 ) / a 0 ;   ε b ( % ) = ( b − b 0 ) / b 0 , where ε _a and ε _b represent strain degrees of the model in x and y direction, a and a ₀ are the strained and original lattice constants in x direction, and b and b ₀ are the strained and original lattice constants in y direction. We successively change the lattice constants a and b and relax models to simulate the effect of stress on the interfaces.

Structures of the P _b defect are shown in Figures 1C,D. Si(1) is the defective atom, and the three silicon atoms bonded with Si(1) are directly connected to the crystalline silicon moiety through Si-Si bonds. In order to simulate the depassivation reaction of the P _b defect at the a-SiO₂/Si (111) interface, the initial and final states of the reaction need to be constructed first. In the initial state of the reaction, the defective atom Si(1) is passivated with the atom H(1), and the proton H(2) is stably connected to the Si-O(1)-Si bridge (Yue et al., 2018). In the final state of the reaction, a hydrogen molecule is formed in the void near the defect. It is worth mentioning that before carrying out further calculations, the stability of the initial and final state structures needs to be confirmed (Pei et al., 2021; Ying-Ying et al., 2021). The atomic vibrational frequencies of initial and final state structures must be calculated first to make sure that there is no imaginary frequency. It is based on this operation that the reaction processes can be simulated correctly in this manuscript. The initial and final states of the P _b defect depassivation reaction are shown in Figure 1.

Since the P _b defect generated by the depassivation is positively charged, the defective Si atom and the surrounding three silicon atoms form a nearly planar structure in the final state, as shown in Figure 1D. The tetrahedral distortion parameter t is defined as t = ( 360 − ∑ i = 1,3 α i ) / 31.5 where ∑ i = 1,3 α i is the sum of three bond angles (Blöchl, 2000). In this manuscript t is used to measure the hybridization of the defective silicon atom. For an sp ³ hybridized silicon atom with a perfect tetrahedral structure, t is one, whereas for an sp ² hybridized silicon atom with a planar trigonal structure, t is zero (Blöchl, 2000). The comparisons between the P _b defect and the P _b ⁺ defect on three Si-Si-Si bond angles and the parameter t are listed in Table 1. It can be seen from Table 1 that the t parameter of the P _b defect is 1.15 when it is neutral, and the defect atom is close to an ideal sp ³ hybridization. When the P _b defect is positively charged as P _b ⁺, the t parameter is 0.09, and the defect atom is close to an ideal sp ² hybridization.

In the initial state 00 of the depassivation reaction, the bond lengths of Si(1)-H(1) and O(1)-H(2) are 1.508 and 0.984 Å, respectively, and the distance between H(1) and H(2) is 2.857 Å. From the initial state 00 to the state 02, H(1) and H(2) tend to move towards each other, and the lengths of Si(1)-H(1) and O(1)-H(2) bonds keep increasing until the bonds are broken. When the reaction proceeds to the transition state 02, the energy of the system reaches a highest point of 1.398 eV. In the 02 state, the distances between Si(1) and H(1), O(1) and H(2), and H(1) and H(2) are 1.889, 1.323, and 1.100 Å, respectively. As the reaction proceeds from state 02 to state 04, the distance between H(1) and H(2) is decreasing, while the distances between Si(1) and H(1), and O(1) and H(2) keep increasing. In the final state 04, the H(1)-H(2) bond length stabilizes at a value of 0.755 Å. The three intermediate states and the energy profile curve are shown in Figure 2. The bond lengths and energies of the states during the depassivation reaction are listed in Table 2.

It can be extracted from Figure 2D that the 02 state is the transition state and the activation energy is 1.398 eV. VASP is hired to calculate the vibration frequencies. There is no imaginary frequency in the initial and final states, and there is one and only one in the transition state. The calculated pre-exponential factor is 1.44 × 10^–7 cm³s⁻¹.

Structures of the P _b1(AOD) defect are shown in Figures 3C,D. Si(1) is the defective atom, and Si(4) as one of the three silicon atoms connected to Si(1), is linked to the a-SiO₂ moiety through one Si-O-Si bridge and to the crystalline silicon moiety through two Si-O-Si bridges. In the initial state of the P _b1(AOD) defect depassivation, the proton H(1) is stably connected to the Si-O(1)-Si bridge to form the O(1)-H(1) bond (Yue et al., 2018), and the atom H(2) is chosen to passivate the defective atom Si(1). In the final state, a hydrogen molecule is formed in the void near the P _b1 ⁺ defect. The initial and final state diagrams are shown in Figure 3.

It can be seen from the final state partial Figure 3D that the defective silicon atom Si(1) and the surrounding three atoms form a plane-like structure. The comparisons between the P _b1 defect and the P _b1 ⁺ defect on three Si-Si-Si bond angles and the parameter t are listed in Table 3. When the P _b1 defect is neutral, the t parameter is 1.25 close to the value of one for a perfectly sp ³ hybridized atom. For the positively charged P _b1 ⁺ defect, the t value is 0.41 close to the perfect sp ² hybridization value zero.

In the initial state 00, the bond lengths of Si(1)-H(2) and O(1)-H(1) are 1.501 Å and 0.989 Å, respectively, and the distance between H(1) and H(2) is 2.665 Å. From the initial state 00 to the intermediate state 01, H(1) gradually separates from O(1), and H(2) gradually separates from Si(1). When the reaction proceeds to the 02 state, the system energy reaches the maximum value of 0.931 eV. In the 02 state, the distances between Si(1) and H(2), O(1) and H(1), and H(1) and H(2) are 1.758, 1.353, and 1.035 Å, respectively. As the reaction progresses, the distance between Si(1) and H(2) keeps increasing, and the distance between O(1) and H(1) is also gradually increasing while the distance between H(1) and H(2) is gradually decreasing. Finally the reaction proceeds to the final state 05, where the H(1)-H(2) bond length is 0.740 Å. The four intermediate states diagrams and the energy profile curve are shown in Figure 4. The bond lengths and energies of the states during the depassivation reaction are listed in Table 4.

In Figure 4E, the 02 state is the transition state, and the activation energy is 0.931 eV. VASP is used to calculate the vibration frequencies and there is no imaginary frequency in the initial and final states, and one and only one in the transition state. The calculated pre-exponential factor is 2.71 × 10^–8 cm³s⁻¹.

In order to simulate the stress induced by the final passivation layer, a biaxial strain is applied to the a-SiO₂/Si (111) and (100) interfaces, and the strain degree ε ranges from −7% to 7% with a step of 1%. The stress-strain scatter diagrams are shown in Figure 5. The external pressure is the equivalent pressure on the interface after stress is applied. A positive value indicates that the interface is compressed, and a negative value indicates that the interface is stretched. It is shown by the figure that the external pressure (or its absolute value) is almost linearly proportional to the strain degree.

In addition to the stress-strain curves, the change of average bond length and average bond angle also reflects the degree of deformation of models to a certain extent. In order to facilitate a more detailed analysis of the structural properties of models, and to further verify the rationality of the stress applied, a statistical analysis of the bond lengths and bond angles is conducted, including the average values and the increase rates. The increase rate is defined as the ratio of the increase value of bond length (or bond angle) in a certain stage to the initial value. Tables 5, 6, respectively show the average lengths of Si-O bond in a-SiO₂ part and Si-Si bond in crystalline silicon part, and the average angles of ∠O-Si-O and ∠Si-O-Si in a-SiO₂ part with strains of −7%, 0% and 7%.

It is shown by the above two tables that, as the strain degree ε of the a-SiO₂/Si interface increases from −7% to 7%, the average values of bond lengths and bond angles also increase. However, the increase rate γ is different. When ε increases in a certain extent, i.e., from −7% to 0% or from 0% to 7%, the Si-Si bond exhibits a higher average increase rate of 2.65%, while that of the Si-O bond is 0.45%. It can be concluded that the Si-Si bond is relatively more sensitive to stress. For the bond angles, ∠Si-O-Si has an average increase rate of 2.68% while ∠O-Si-O has an average increase rate of 0.10%. It can be seen that compared to ∠O-Si-O, ∠Si-O-Si is more sensitive to stress.

After applying stress to the interface models, the depassivation reactions of the interface defects under different stresses are investigated. Tensile and compressive strains are applied to the a-SiO₂/Si (111) and (100) interfaces and the activation energies of P _b and P _b1 defects depassivations under different strains are listed in Tables 7, 8. The energy curves during depassivation reactions and activation energy profile curves (Figures 6–9) under different stresses are shown below.

From Figure 6 and Table 7, it can be concluded that when the tensile strain degree ε is between 0% and 4%, the activation energy of P _b defect depassivation reaction is stable at about 1.35 eV, and the fluctuation range is less than 0.07 eV. When the strain degree ε increases from 4% to 5%, the activation energy surges by 0.76–2.05 eV. Afterwards as the strain degree ε increases, the activation energy stabilizes at a higher value of about 2.00 eV.

In Figure 7 and Table 7, as the strain degree ε decreases from 0% to −3%, the P _b depassivation reaction activation energy stabilizes at about 1.40 eV, and the fluctuation is less than 0.07 eV. Within this strain range, the forward reaction barrier of the depassivation reaction is always higher than the reverse one. When the strain degree ε decreases from −3% to −4%, the activation energy drops drastically to 1.17 eV, and the forward reaction barrier becomes smaller than the reverse reaction barrier. Then, as the strain degree ε decreases from −4% to −7%, the reaction activation energy also gradually decreases and is less than 1 eV, and the forward reaction barrier of the reaction is always smaller than the reverse one.

In Figure 8 and Table 8, as the tensile strain degree ε increases, the activation energy of the P _b1 defect depassivation reaction does not change much. It is basically stable at about 0.90 eV and the fluctuation is less than 0.08 eV. From the depassivation reaction curves, when ε increases to 6%, the reverse reaction barrier of the depassivation reaction becomes higher than the forward reaction barrier.

In Figure 9 and Table 8, the activation energy of the depassivation reaction of P _b1 defects is stable at about 0.92 eV under compressive strain, and the fluctuation is less than 0.07 eV. As the degree of compressive strain increases, the reverse reaction barrier becomes higher than the forward one.

According to Eq. 8, as long as the activation energy and the pre-exponential factor are obtained, the reaction rate k _d can be solved for a chemical reaction. The pre-exponential factor is the ratio of the product of the vibration frequencies of the initial state and the product of all real vibration frequencies of the transition state of the reaction. In this study, VASP code is used to calculate pre-exponential factors of different depassivation reactions. Combined with the reaction activation energy shown in the previous part, the depassivation reaction rates of P _b and P _b1 defects under different strain degrees at room temperature (300K) are calculated. The results are exhibited in Tables 9, 10 and curves are plotted in Figure 10.

It can be seen that at room temperature the depassivation reaction rate of P _b1 defects is generally faster than that of P _b defects. For P _b defects, when the strain degree ε is 5%, the depassivation reaction rate is the slowest; when ε is −7%, the depassivation reaction rate is the fastest. In addition, when the strain degree ε is in the range of −7%∼−5%, the depassivation reaction rate is relatively fast, which is 8–11 orders of magnitude faster than that without stress; when ε is in the range of 5%–7%, the depassivation reaction rate is relatively slow, which is 10–12 orders of magnitude slower than that in the case of no stress. For P _b1 defects, the depassivation reaction rate changes generally steady, and there is no significant increase or decrease as stress is applied. When the strain degree ε = -2%, the depassivation reaction rate is the fastest, which is two orders of magnitude faster than that without stress. The increase or decrease of depassivation rate under other strain degrees is insignificant. It can be inferred that the depassivation of P _b1 defects on the a-SiO₂/Si (100) interface is insusceptible of stress, while the depassivation of P _b defects on the a-SiO₂/Si (111) interface can be greatly affected by stress.

After obtaining the depassivation reaction rate constant k _d, the curve of the relative concentration value (x) of the defect with time (t) can be calculated according to Eq. 5. Researchers have studied the hydrogen/deuterium diffusion in solid-state crystallized (SSC) poly-Si and low-pressure chemical-vapor-deposited (LPCVD)-grown polycrystalline silicon. The study pointed out that the H surface concentration amounts to ≈10¹⁹ and ≈3×10¹⁹ cm⁻³ in LPCVD-grown and SSC poly-Si, respectively (Nickel et al., 1996). Based on the previous research [H⁺] in Eq. 5 is set to 10¹⁹ cm⁻³ in our work, and the relative concentration x _Pb and x _Pb1 of P _b and P _b1 defects as a function of time are plotted in Figure 11.

It can be seen from Figure 11 that compared with the P _b defect, the depassivation of the P _b1 defect is relatively insensitive to stress, which can be explained through their local structures. For the P _b1 defect, one of the three silicon atoms bonded to the defective atom is also bonded to three oxygen atoms, so the defective silicon atom is connected to the silicon substrate through a [·Si-Si-(O)₂-Si] structure, whereas the defective silicon atom of the P _b defect is directly connected to the silicon substrate through a [·Si-Si] structure. In a nutshell, the Si-O-Si bridge exists in the P _b1 defect, but not in the P _b defect. The different depassivation performances for the two defects under stress are owing to the softness of the Si-O-Si bridge. Compared to the Si-Si-Si bond angle, the energy required to change the Si-O-Si bond angle is much lower, so the P _b1 defect can release local stress through a torsion of its structure, therefore making itself insensitive to stress (Li et al., 2019b). However, the P _b defect is short of a structure like the Si-O-Si bridge, so it is more sensitive to the applied stress.

For P _b1H defects at a-SiO₂/Si(100) interface, it takes up to 10³–10⁶ s to depassivate them with protons to P _b1 ⁺ at room temperature, which is consistent with the experimental results summarized by Schwank et al. (2008) on the accumulation rate of interface traps. When the strain degree ε of the a-SiO₂/Si(100) interface is −2%, 1% or 2% (equivalent to the interface pressure of about 1.36 GPa, −0.30 GPa or −0.80 GPa), the depassivation reaction rate of P _b1H defects is the fastest. The total time for protons to depassivate all P _b1H into P _b1 ⁺ is about 10³–10⁴s. When ε is −3%, 6% or 7% (equivalent to the interface pressure of about 1.78 GPa, −2.30 GPa or −2.60 GPa), the depassivation reaction rate of P _b1H defects is the slowest. At this time, it takes about 10⁵–10⁶ s for the P _b1H depassivation reaction to reach saturation.

For P _bH defects at a-SiO₂/Si (111) interface, stress greatly impacts the time taken for the depassivation reaction. When the strain degree ε of a-SiO₂/Si (111) interface is −5%, −6% or −7% (equivalent to the interface pressure of about 3.20, 3.50 or 3.90 GPa), the P _bH depassivation reaction rate is relatively high, and all P _bH depassivation reactions can reach saturation within 10⁵ s, especially when ε is −6% or −7%, it takes only 10² s to depassivate all P _bH into P _b ⁺ defects. When the strain degree ε of a-SiO₂/Si (111) interface is 5%, 6%, 7% (equivalent to the interface pressure of about −1.90, −2.30 GPa or −2.60 GPa), the P _bH depassivation reaction rate is very slow. In these cases, the time needed to depassivate all P _bH at the interface to P _b ⁺ is about 10²³∼10²⁴ s.