Bond Insertion at Distorted Si(001) Subsurface Atoms

Abstract

Using density functional theory (DFT) methods, we analyze the adsorption of acetylene and ethylene on the Si(001) surface in an unusual bond insertion mode. The insertion takes place at a saturated tetravalent silicon atom and the insight gained can thus be transferred to other saturated silicon compounds in molecular and surface chemistry. Molecular orbital analysis reveals that the distorted and symmetry-reduced coordination of the silicon atoms involved due to surface reconstruction raises the electrophilicity and, additionally, makes certain σ bond orbitals more accessible. The affinity towards bond insertion is, therefore, caused by the structural constraints of the surface. Additionally, periodic energy decomposition analysis (pEDA) is used to explain why the bond insertion structure is much more stable for acetylene than for ethylene. The increased acceptor abilities of acetylene due to the presence of two π*-orbitals (instead of one π*-orbital and a set of σ*(C–H) orbitals for ethylene), as well as the lower number of hydrogen atoms, which leads to reduced Pauli repulsion with the surface, are identified as the main causes. While our findings imply that this structure might be an intermediate in the adsorption of acetylene on Si(001), the predicted product distributions are in contradiction to the experimental findings. This is critically discussed and suggestions to resolve this issue are given.

1. Introduction

The discovery of stable molecules containing silicon–silicon double bonds four decades ago [1] has sparked the research interest in the chemistry of disilenes [2,3,4,5]. Unlike in alkenes, the substituents at these double bonds are not arranged in a planar fashion, but tilted, a result of the bonding situation which has been described as a double dative bond [6,7]. The most stable conformation is a trans arrangement of the four substituents [6,8,9,10], which is favored not only from orbital overlap aspects at the Si–Si bond, but also Pauli repulsion within the molecule [11]. While these unusual double bonds can undergo similar reactions as alkenes, e.g., 1,2-addition, cycloaddition, or coordination in metal complexes [3], the reactivity is increased and sterically demanding substituents are needed to stabilize the molecules [5]. Additionally, the tilted geometry leads to changes in reaction mechanisms and stereochemistry, as a recent study on the addition of molecular hydrogen to a disilene showed [12]: here, the system favors anti-addition to syn-addition or a stepwise mechanism.

A cis arrangement of substituents in disilenes is rarely observed. Two examples include silicon atoms as part of a borane cluster [13] and a base-stabilized adduct [14]. In both cases, the coordination of the two silicon atoms in the most stable conformation is not symmetric (formally, C_2v symmetry, in the case of four identical substituents), but asymmetric (formally, C_s symmetry). This distortion changes the electronic structure and, as a consequence, the reactivity: in the base-stabilized adduct [14], one silicon atom acts as an acceptor in a dative bond, a situation that would not be possible in carbon–carbon double bonds. Further investigation of cis-substituted disilenes could, therefore, reveal new aspects of silicon chemistry.

The Si(001) surface, a widely-used substrate in surface science with a high relevance for application [15], features a structure similar to cis-substituted disilenes (Figure 1): in the surface reconstruction process of the bulk (a), which crystalizes in a diamond structure, two adjacent atoms form a covalent bond and tilt, yielding the characteristic buckled dimers of this surface (b). Although the bond could be, in a local picture, formally described as a double bond, the buckling leads to a localization of the occupied π-type orbital at the upper Si_up atom, while the empty π*-type orbital is mostly localized at the lower Si_down atom (c) [16]. Therefore, the surface possesses both nucleophilic (at Si_up) and electrophilic (at Si_down) characters. This enhances the reactivity with organic molecules: In addition to typical double bond reactions like 1,2-addition and cycloaddition, the surface can form dative bonds and act as a reagent in nucleophilic substitution [17,18]. This underlines that, in many cases, it can be treated as a molecular reagent [19]. The study of reactions on Si(001) can, therefore, aid in the understanding of molecular disilenes.

Previously, a new reaction type on Si(001) was reported in a theoretical study [20]: the bond insertion of acetylene into a Si_down–Si_sub bond (compare Figure 1b), called the sublayer mode from now on. This behavior is highly unusual, since, until now, the Si_up and Si_down atoms were considered reactive centers in most cases. Si_sub atoms were disregarded because they are tetravalent, therefore fully saturated and missing the “dangling bond” orbitals of Si_up and Si_down (Figure 1c). The observation is also in stark contrast to the adsorption of alkenes like ethylene, which undergo a cycloaddition with one or two surface dimers mediated by a coordination of the molecular π system to a Si_down atom (π complex) [21,22,23,24,25,26,27]. For the adsorption of acetylene, cycloadducts are the most stable adsorption modes as well, as confirmed by theoretical calculations and experimental measurements at room temperature [20,21,28,29,30,31,32,33,34,35,36]. While there is experimental evidence for the presence of an intermediate state (called precursor in surface science) [37], there is no consensus on the nature of this state. Statements that have been made include a mobile precursor [37], a π complex [28,30], no stable precursor at all [21,35], and the bond insertion structure [20].

Using molecular orbital analysis and bonding analysis, we will now provide a description of acetylene adsorption on Si(001) that advocates the insertion mechanism (Scheme 1). Even though this mechanism has not been reported for other systems so far, we will explain how the electronic structure of acetylene benefits this reaction by comparing it to the adsorption of ethylene. Additionally, we will provide an explanation for the reactivity of Si_sub atoms by molecular orbital analysis of distorted silane. This will show that the calculated reactivity is a reasonable alternative to existing hypotheses.

2. Computational Details

2.1. Molecules

Structural optimization was performed using the PBE functional [38,39], the D3(BJ) dispersion correction [40,41] and the def2-TZVPP basis set [42] in ORCA 3.0.2 [43] (standard convergence criteria and integration grids). Molecular orbitals were calculated and visualized using PBE-D3(BJ)/TZ2P in ADF 2016 [44,45,46] (standard convergence criteria and integration grids).

2.2. Surfaces

All energies and structures were calculated with the Vienna Ab Initio Simulation Package (VASP) [47,48,49,50] version 5.3.5 using the PBE and HSE06 [51] functionals (structures and frequencies: PBE only), the D3(BJ) dispersion correction and the PAW formalism [52,53] with a basis set cutoff of 400 eV. Electronic k space was sampled using a Γ(221) grid. Self-consistent field (SCF) calculations were converged to an accuracy of 10⁻⁶ eV, and structural optimizations to 10⁻² eV·Å⁻¹. Reaction paths were calculated using the Climbing-Image Nudged Elastic Band method [54], while the transition state structures were refined using the Dimer method [55]. Harmonic vibrational frequencies were calculated from structures converged to 10⁻³ eV·Å⁻¹ by a finite differences approach and construction of the Hessian using Cartesian displacements of 0.01 Å (SCF convergence: 10⁻⁸ eV). Gibbs energies were calculated at T = 300 K, p = 1 bar in an approach described elsewhere [27]. All minimum and transition state geometries and their respective total energies are given in the Supplementary Materials.

The Si(001) surface was modeled as a six-layer slab in c(4 × 2) reconstruction and frozen double layer approximation. The bottom was saturated with hydrogen atoms in a tetrahedral arrangement at a distance of d(Si–H) = 1.480 Å, the experimental equilibrium distance in silane [56]. Cell parameters a and b were set to 15.324 Å, corresponding to a 4 × 4 cell with an optimized bulk lattice parameter of 5.418 Å. In the c direction, a vacuum layer of at least 10 Å was ensured. Convergence studies on these parameters can be found in a previous work [57]. The bonding energy E_bond was defined as the difference between the energy E_tot of the relaxed total system and the energies E_surf and E_mol of the relaxed and isolated surface and molecule:

E_bond = E_tot − E_surf − E_mol

(1)

Please note that surface science convention is the use of the adsorption energy E_ads with inverse sign convention (E_ads = −E_bond). Scanning tunneling microscopy (STM) topographies were calculated in the Tersoff-Hamann approximation [58,59] using bSKAN [60,61]. The approach outlined has delivered accurate and reliable results for organic/semiconductor systems in the past [18,26,27,57,62,63].

2.3. Periodic Energy Decomposition Analysis

Bonding analysis was performed at PBE-D3(BJ)/TZ2P, Γ only k sampling, using pEDA [64] in ADF-BAND 2016 [65,66]. The pEDA method allows to dissect E_bond into well-defined quantities that allow interpreting the bonding between two fragments in a system (here: molecule and surface) in a chemically meaningful way. In the first step, E_bond is partitioned into the intrinsic interaction energy ∆E_int and the respective preparation energies ∆E_prep of the molecule (M) and surface (S):

E_bond = ∆E_int + ∆E_prep(M) + ∆E_prep(S)

(2)

Since a dispersion correction is applied to the DFT calculations, ∆E_int can be divided into a dispersion term (disp) and an electronic term (elec):

∆E_int = ∆E_int(disp) + ∆E_int(elec)

(3)

The actual pEDA procedure then decomposes ∆E_int(elec) into contributions from Pauli repulsion (∆E_Pauli), electrostatics (∆E_elstat), and orbital interaction (∆E_orb):

∆E_int(elec) = ∆E_Pauli + ∆E_elstat + ∆E_orb

(4)

Additionally, the NOCV extension [67,68] was applied to the pEDA calculation in this work. This allows expressing ∆E_orb as the sum of individual energy contributions ∆E_orb(i) of different character (e.g., σ/π bonding, donation and back donation). The assignment to a character was done by visual inspection of the deformation densities ∆ρ(i) which show the corresponding charge transfer.

3. Results

3.1. Distorted Coordination of Tetravalent Silicon

The enhanced reactivity of Si_sub atoms can be understood from the structural distortion in the reconstruction process (Figure 1a): since the position of the surface atoms Si_up and Si_down is changed with respect to the bulk, the coordination of the Si_sub atoms changes from tetrahedral to a symmetry-reduced arrangement (Figure 1b). This affects the electronic structure, as can be seen from a molecular orbital analysis of a similarly distorted silane molecule: Figure 2a shows a Walsh diagram [69] depicting the Kohn-Sham energies of selected molecular orbitals as a function of the distortion angle α (Figure 2b), i.e., the tilting of two Si–H bonds from tetrahedral position without a change in bond length. The distortion reduces the symmetry of the molecule from point group T_d to C_s and causes degenerate orbitals to split in energy. This is most evident in the triply-degenerate 1t₂ orbital, the HOMO of the tetrahedral molecule: Orbital 4a′ rises approximately linearly in energy with increasing α, while at the same time, orbital 1a″ is unaffected and 3a′ decreases in energy. This is caused by changes in overlap of the involved 3p(Si) and 1s(H) atomic orbitals. The unoccupied and anti-bonding 2t₂ orbitals are affected in a similar manner and the LUMO of the distorted molecule (5a′) gets significantly lower in energy with increasing α. At 50° distortion, it features a large lobe at the “empty coordination site”, i.e., the location of a fifth bonding partner in a five-fold coordinated structure. Since this orbital becomes spatially and energetically more accessible, the distortion makes the system more electrophilic. The ability of silicon to accommodate more than four substituents is well-known in molecular silicon chemistry.

Still, this does not explain the activation of the Si_down–Si_sub bonds, since the HOMO (4a′), while raised in energy, is located at three Si–H bonds. However, on Si(001) the structure is even further reduced in symmetry and the two distorted bonds are not tilted by the same amount. This is because every Si_sub atom is connected to one Si_down and one Si_up in the minimum configuration (Figure 1b) and the corresponding bonds are tilted by a different α: Approximately 30° (towards Si_up) and 50° (towards Si_down). When this geometry is reproduced in silane, an additional symmetry reduction takes place (C_s to C₁), the HOMO (Figure 3a) becomes a mixture of orbitals 4a′ and 1a″ of the C_s-symmetric molecule and features a large lobe at the Si–H bond that is tilted by 50°, i.e., the one representing the Si_down–Si_sub bond on Si(001). At the same time, the LUMO (Figure 3b) is essentially unchanged with respect to the C_s-symmetric molecule. This asymmetric distortion provides a much better explanation of the bond activation, since the HOMO could interact with an empty orbital of another molecule (adsorbate) and form a covalent bond by breaking the Si–H (Si_down–Si_sub) bond, whereas the LUMO could interact with an occupied orbital and form a second covalent bond, reinforcing the bond insertion reactivity.

Two selected crystal orbitals (COs) of a Si(001) slab at the Γ point in k space show that the electronic structure at Si_sub atoms is similar: the HOCO (Figure 3c, −4.85 eV) is localized primarily at the Si_down–Si_sub bonds, while the LUCO+5 (Figure 3d, −3.47 eV) features large lobes at the Si_sub atoms. It is surprising to see the Si_down–Si_sub bond orbital as the HOCO, since the non-bonding electron pairs at Si_up (Figure 1c) are usually considered as frontier orbitals in the surface reactivity of Si(001). However, at the Γ point, the highest energy orbital representing the Si_up states (HOCO–1) is very close in energy (−5.02 eV) to the HOCO. Additionally, this orbital becomes the HOCO in other parts of k space, as a recent theoretical determination of the band structure using hybrid functionals has shown (Figure 5c in [70]). Hence, both the non-bonding electron pairs at Si_up and the activated Si_down–Si_sub bonds are high enough in energy to be considered reactive towards adsorbates. The energetic location of the LUCO+5, the other crystal orbital involved in bond insertion reactivity, is more in line with previous considerations, since the orbitals representing the empty states at Si_down (Figure 1c) are still lower in energy (−4.43 to −3.61 eV) and, therefore, more accessible. All of this shows that the effects described for silane also apply to Si_sub atoms on the Si(001) surface and that the enhanced reactivity of these saturated atoms is caused by the distorted geometric arrangement.

3.2. Bond Insertion of Acetylene and Ethylene on Si(001)

The optimized structures and bonding energies of acetylene and ethylene on Si(001) in the π complex and sublayer adsorption modes are given in Figure 4. As previously reported [20], the π complex mode of acetylene is not a minimum for the computational approach chosen here, but a first-order saddle point (ν_imag = 71i). This is also evident in the energy profile connecting the two stationary points (c) which shows no local energy maximum. The acetylene π complex structure can therefore be seen as a transition state (TS) between two symmetry-equivalent sublayer minima. In contrast to this, the ethylene π complex structure is a minimum, the sublayer mode is higher in energy, and the energy profile features a maximum. Furthermore, a local minimum is apparent: this represents the lowest energy rotamer of the ethylene π complex, in which the C–C bond is rotated by 30–45° approximately in the xy plane [27,71]. Since the structure shown in Figure 4b is only 3 kJ·mol⁻¹ higher in energy and allows a better comparability with the corresponding acetylene TS structure due to a similar molecular orientation, it is analyzed here instead of the minimum rotamer.

The two π complex structures show several similarities, beginning with a C–C bond elongation by 3–4% compared to the gas phase value (acetylene: 1.208 → 1.239 Å, ethylene: 1.333 → 1.381 Å). This is already an indication of electron donation from the π system to the surface. An indicator of back donation is the Si_down–Si_sub bond length, which is elongated by 2% compared to its minimum value in the isolated surface slab (2.340 → 2.396 (acetylene), 2.340 → 2.389 Å (ethylene)). The nature of this interaction would be a donation of electron density from the σ(Si_down–Si_sub) crystal orbital (Figure 3c) to the molecular π* orbital. Carbon-silicon distances d(C–Si_down) are in a similar range for both systems as well but slightly shorter for acetylene (2.175 Å, ethylene: 2.234 Å). While this could indicate a stronger interaction between molecule and surface in comparison to ethylene, the E_bond value of acetylene (−64 kJ·mol⁻¹) implies a weaker interaction (ethylene: −71 kJ·mol⁻¹). Bonding analysis will provide an explanation for this apparent contradiction later on. In comparison with [20], where this reaction was reported first, our E_bond values are lower in energy by ~25 kJ·mol⁻¹. This difference can be mainly attributed to the dispersion correction to the DFT energies, which was not applied in the literature study, and other differences in the computational setup. Our G_bond values (Figure 4) are higher in energy by 47 (acetylene) and 58 kJ·mol⁻¹ (ethylene) than E_bond values, highlighting the importance of molecular entropy loss upon adsorption described previously [27,72,73]. However, it should be noted that, in the case of the acetylene TS, the partition function used to calculate G_bond has one fewer degree of freedom, the imaginary mode. Hence, this G_bond value should not be taken as a bonding energy of this state in thermodynamic equilibrium, which is not possible for a saddle point, but rather the effective energy in a reaction where the structure appears as a TS, i.e., conversion between two sublayer minima.

The sublayer mode is easily characterized by a drastically increased Si_down–Si_sub bond length (acetylene: 3.128 Å, ethylene: 3.127 Å) with respect to the isolated slab value (2.340 Å), a strong indication that this covalent bond has been broken. At the same time, two C–Si covalent bonds are formed (d(C–Si_down/C′–Si_sub) = 1.848/1.910 (acetylene), 1.882/1.944 Å (ethylene)) and the C–C′ bond elongates to 1.348 (acetylene) and 1.529 Å (ethylene), typical values for double and single bonds, respectively. This implies that one π bond is broken in the insertion process and the C–C′ bond order is reduced by one. The slightly larger d(C–Si) values for ethylene are caused by Pauli repulsion between the CH₂ groups and the surface, as bonding analysis will show later. While the acetylene sublayer mode is markedly lower in energy than the π complex TS at values of −120 (E_bond) and −65 kJ·mol⁻¹ (G_bond), the ethylene structure is higher in energy than the π complex at −65 (E_bond) and −2 kJ·mol⁻¹ (G_bond). This was reported in the previous study [20] as well (E_bond = −97 (acetylene) and −22 kJ·mol⁻¹ (ethylene)) and used as an argument on why this mode was not observed in any adsorption experiments of ethylene. The differences of the literature energies to our values can again be attributed to the missing dispersion correction and the different computational setup. In particular, the smaller cell size of 2 × 2 atoms per layer (this study: 4 × 4) appears to destabilize the sublayer modes. A similar effect was previously found for the bridge mode of ethylene [27]. This destabilization is also apparent at the TS for the π complex → sublayer conversion in the ethylene system, since our values for the energy barrier (E_a: 14, G_a: 13 kJ·mol⁻¹) are considerably lower than the literature value of E_a = 41 kJ·mol⁻¹.

Considering that both this study and [20] were performed using GGA functionals, one could argue that the stability of the acetylene sublayer mode might be an artifact of the methodology. To address this issue, we performed single point calculations at the optimized structures using the HSE06 hybrid functional. The resulting energies of the stationary points (Figure 4c) show only a minor decrease in E_bond values of up to 15 kJ·mol⁻¹ and barely any changes in energy differences (HSE06-D3 energies of the full reaction paths in Figure 4c are given in the Supplementary Materials). This confirms that the PBE functional is an appropriate choice for the description of these systems. An independent test of the DFT energy values, e.g., by wavefunction-based methods, would nevertheless be desirable.

Bonding analysis results obtained from pEDA calculations (see Section 2.3) are given in

Table 1. pEDA bonding analysis of acetylene and ethylene adsorbed on Si(001) in the π complex and sublayer modes ¹.

	π Complex				Sublayer
	Acetylene (TS)		Ethylene		Acetylene		Ethylene
∆Eint	−105		−106		−656		−583
∆Eint(disp) 2	−23	(22%)	−27	(25%)	−22	(3%)	−30	(5%)
∆Eint(elec) 2	−82	(78%)	−79	(75%)	−634	(97%)	−553	(95%)
∆EPauli	701		653		1587		1628
∆Eelstat 3	−367	(47%)	−342	(47%)	−912	(41%)	−942	(43%)
∆Eorb 3	−416	(53%)	−390	(53%)	−1310	(59%)	−1240	(57%)
∆Eorb(M → S) 4	−292	(70%)	−289	(74%)	−551	(42%)	−536	(43%)
∆Eorb(S → M) 4	−88	(21%)	−76	(19%)	−714	(55%)	−641	(52%)
∆Eprep(M)	15		10		364		341
∆Eprep(S)	17		16		156		161
Ebond 5	−73	(−64)	−80	(−71)	−136	(−120)	−81	(−65)

¹ All values in kJ·mol⁻¹, calculated using PBE-D3/TZ2P. Fragments: Molecule and surface. Fragmentation: Closed-shell singlet (π complex), triplet (sublayer). M: Molecule, S: Surface. ² Percentage values give relative contributions of dispersion and electronic effects to the interaction energy ∆E_int. ³ Percentage values give relative contributions of the attractive pEDA terms ∆E_elstat and ∆E_orb. ⁴ Percentage values give relative contributions to the total orbital interaction energy ∆E_orb. ⁵ PAW values (in parentheses) given for comparison.

. By comparing the π complex values, it becomes apparent that the lower E_bond value of ethylene is caused by a reduced molecular preparation energy (∆E_prep(M) = 10 vs. acetylene: 15 kJ·mol⁻¹) and an increased dispersion interaction (∆E_int(disp) = −27 vs. acetylene: −23 kJ·mol⁻¹). The electronic interaction energy ∆E_int(elec), as well as the stabilizing pEDA terms electrostatic (∆E_elstat) and orbital interaction (∆E_orb) are all weaker in case of ethylene. These results resolve the apparent contradiction of shorter C^1/2–Si_down bonds, but less negative E_bond values in the case of acetylene: The shorter distance is caused by a larger stabilization through electrostatic and orbital interaction, which are overcompensated by the reduced dispersion interaction and higher preparation energy. However, the character of the chemical bond as given by pEDA is comparable in both structures: a dative bond between the molecule and surface, apparent from the ∆E_elstat/∆E_orb ratio (47:53), a typical value in covalent bonds [74], and a dominating molecule-to-surface orbital contribution ∆E_orb(M → S) (acetylene: 70% of ∆E_orb, ethylene: 74%).

In the sublayer mode of acetylene, ∆E_int is drastically lower than in the π complex TS, at −656 kJ·mol⁻¹. Since ∆E_int(disp) barely changes, this is an electronic effect that can be attributed to the formation of the two shared-electron bonds between molecule and surface. ∆E_Pauli, ∆E_elstat, and ∆E_orb rise significantly in absolute value, and orbital interaction is now even more predominant in stabilization at 59% vs. 41% for ∆E_elstat. This is mainly caused by a large increase in back bonding (surface to molecule) contributions ∆E_orb(S → M), now the dominating orbital term (55% of ∆E_orb). The increase occurs because the nature and number of the covalent bonds has changed from one dative to two shared-electron bonds. These major contributions to stabilization, however, are partially compensated by large preparation energies in both molecule and surface (364 and 156 kJ·mol⁻¹, respectively), so that the resulting E_bond value is only ~60 kJ·mol⁻¹ lower than the energy of the π complex TS. The ethylene sublayer analysis shows that the destabilization with respect to acetylene is not a result of increased molecular deformation, as previously suggested [20], since the ∆E_prep(M) value for ethylene is lower than the corresponding value for acetylene. Furthermore, ∆E_int(disp) is lower for ethylene as well, so the effect is entirely contained in ∆E_int(elec), which is 81 kJ·mol⁻¹ less stabilizing than for acetylene. The pEDA decomposition reveals that the destabilization is contained in ∆E_Pauli and ∆E_orb, but not ∆E_elstat, which is more stabilizing for ethylene. This is in contrast to the π complex structures, where the weaker ∆E_int was reflected in a lower absolute value of all three pEDA terms. Here, the ∆E_orb destabilization is mainly caused by a reduction of back bonding contributions ∆E_orb(S → M), which are weaker by 73 kJ·mol⁻¹ for ethylene, whereas ∆E_orb(M → S) contributions are only reduced by 15 kJ·mol⁻¹. This can be understood from the electronic structure of the molecules: Whereas acetylene has two low-lying π*-orbitals that can act as acceptors for electron density from the surface, in ethylene, only one π*-orbital is available. The back donation into σ*(C–H) bonds (negative hyperconjugation) that is possible in ethylene yields much less stabilization since the orbitals are higher in energy. Therefore, the acceptor ability of ethylene is reduced (lower ∆E_orb(S → M) value) and the amount of ∆E_orb stabilization decreases, while, at the same time, Pauli repulsion increases since the two additional C–H bonds are interacting with the surface. The only stabilizing effects that arise from the additional atoms and bonds in ethylene are in ∆E_elstat and ∆E_int(disp), but these are not able to compensate the energy loss from the two other terms. It can therefore be concluded that the enhanced acceptor abilities of acetylene combined with the lower number of C–H bonds leading to the reduced Pauli repulsion are the determining factors that benefit a stabilization of the sublayer mode for this molecule.

3.3. Comparison with Experiment

As already mentioned in the Introduction, these findings imply that the sublayer mode could be a precursor intermediate in the adsorption mechanism of acetylene on Si(001). Up to now, scanning tunneling microscopy (STM) and spectroscopic measurements were performed at room temperature or higher [29,31,32,33], where the more stable on-top and bridge modes (Scheme 1) dominate and all molecules that could have possibly been trapped in a precursor have already reacted to these states. However, if one assumes that these reactions are under kinetic control, as observed for many reactions of organic molecules on Si(001) [19], the distribution of on-top and bridge structures at room temperature can be estimated from calculated energy barriers G_a (Equation (5)) and compared with experimental findings.

$$\frac{N_{on-top}}{N_{bridge}}=\exp \left(-\frac{G_{\mathrm{a},on-top}-G_{\mathrm{a},bridge}}{k_{\mathrm{B}}T}\right)$$

(5)

For ethylene, this approach showed excellent agreement with experiment [27]. Here, the calculated barriers G_a,on-top/G_a,bridge, 55/35 (PBE-D3) and 62/48 kJ·mol⁻¹ (HSE06-D3), yield distributions of N_on-top/N_bridge = 0.0003 (PBE-D3) and 0.004 (HSE06-D3) at 300 K. Unfortunately, an experimental STM measurement at the same temperature showed a much higher abundance of on-top structures (28%) at low coverage and the corresponding N_on-top/N_bridge value of 0.39 would indicate a difference in G_a values of less than 3 kJ·mol⁻¹. Hence, our proposed model for the adsorption mechanism (Scheme 1) is not able to reproduce the experimental findings with the calculated energies and is therefore either incorrect or incomplete. The alternative of formulating a π complex as a precursor from which the system converts to on-top and bridge is not supported by DFT calculations as well: previous studies were unable to determine a transition state on the π complex → on-top path, whereas for the path towards bridge, a transition state could be found [21,35]. Hence, conversion to on-top should be highly preferred, which is again not in line with the experimental observations. To minimize the possibility that we missed other precursor structures by our structural optimization approach, we carried out ten ab initio molecular dynamics simulations of the gas-surface dynamics. Eight trajectories ended up in the sublayer mode while the other two ended up directly in on-top and bridge, respectively (see the Supplementary Materials). Thermal effects are also an unlikely source of the disagreement, since molecule and surface were initially put in random states of thermal excitation corresponding to T = 300 K in these simulations. Additionally, our DFT approach provided reliable results in the description of adsorption kinetics for ethylene and other organic molecules on Si(001) [18,27,62]. Thus, we assume that the disagreement between theory and experiment is not due to a failure of the approximations in our density functional approach, but rather a result of an incomplete model of the adsorption dynamics. An alternative idea would be to assume that the reaction is under thermodynamic control. This also does not deliver a coherent picture, since all theoretical studies predict the on-top mode to be more stable than bridge at low coverage [20,21,30,35,57]. Therefore, a new, presumably more complex model of the adsorption dynamics must be devised, which we cannot provide at this moment.

As an aid to identify different possible precursor species in experimental measurements, Figure 5 shows simulated STM images of the adsorption modes sublayer, on-top, and bridge in comparison with the clean Si(001) surface. The bias voltages were chosen as −2.0, +0.8 and +1.9 V, typical experimental values. While at negative voltage, occupied orbitals are probed and the non-bonding electron pairs at Si_up atoms are visible on the clean surface (compare Figure 1), unoccupied orbitals are probed at a positive voltage and the empty orbitals mostly localized at Si_down atoms are visible. At a voltage of +1.9 V, the nodal planes of the π* type orbitals become apparent. Adsorption modes on-top and bridge are easily identifiable at all voltages by the quenched signals of surface atoms, thereby appearing as dark spots of different orientation [31]. In contrast, the sublayer mode is invisible at negative bias voltage, as is the ethylene π complex structure [26]. However, in contrast to the π complex, which is evident at positive bias voltage from a quenched surface atom signal, this mode might be identified from a displaced Si_down signal, which is due to the change in coordination and electronic structure of this surface atom. However, if the STM resolution does not allow the identification of the displacement of the Si_down signal, the structure would appear invisible at all commonly-used voltages. If the π complex of acetylene is, in contrast to all theoretical investigations, the precursor, it should appear in STM measurements with the same features as the equivalent ethylene structure described above. It should be noted that in order to identify this mode, experimentally demanding low-temperature measurements at both positive and negative bias voltages would be ideal, since the signals at positive voltages could also stem from on-top structures, whereas the molecule is invisible at negative voltages. On the basis of our results, signals that are visible at positive and invisible at negative voltages would be strong indications for the π complex.

Alternatively, vibrational spectroscopy measurements at low temperatures could provide indications about the nature of the precursor. The π complex of ethylene is evident from hindered vibrations at ~200 and 600 cm⁻¹ [22,26] and it can be expected that an equivalent acetylene structure would show a similar vibrational fingerprint. The sublayer mode, on the contrary, could be identified from a large ~100 cm⁻¹ splitting of the C–H stretching vibration band whereas on-top and bridge show a more narrow splitting of ~20–30 cm⁻¹.

An alternative to the precursor-based models discussed up to now could be the existence of direct or pseudo-direct adsorption pathways as present in the adsorption of cylooctyne on Si(001) [62,75]. Future experiments or more accurate computations might resolve the remaining open questions and help to address the difficulties of established DFT methods in describing the adsorption behavior of acetylene on Si(001).

4. Conclusions

Using molecular orbital analysis of silane, we have shown that the geometric distortion of saturated silicon atoms leads to enhanced reactivity. In particular, the electrophilic character of the silicon atom is increased while the distorted σ bonds are activated. On Si(001), this distortion is enforced by the surface reconstruction and leads to an unusually high reactivity of the saturated subsurface atoms. Similar behavior might be observed in molecular systems if steric or electronic effects of substituents could be designed to lead to similar structures. The results therefore present a type of bond activation not driven by orbital interaction, as is often the case, but rather the constrained geometry of the system.

Furthermore, pEDA analysis has shown that the increased acceptor ability and smaller number of C–H bonds of acetylene lead to a stabilization of the sublayer bond insertion mode on Si(001), whereas in the case of ethylene adsorption, this bonding pattern leads to a rather unstable structure. While this mode is a possible candidate for a precursor intermediate in the adsorption of acetylene, the results are not in agreement with the experimental findings. In particular, the ratio of product states on-top and bridge at low coverage deviates by two to three orders of magnitude. However, alternative pathways also do not lead to a better agreement and further experiments or more accurate computations are needed to settle the adsorption mechanism. We provided STM simulations of the sublayer adsorption mode in comparison with the product states which might help in identifying molecules that might be trapped in this intermediate state at low temperatures. In any case, this system provides an opportunity to improve and revise either established adsorption models and/or the theoretical methods used to model them.