Molecular Sciences Model, First-principle Calculation of Ammonia Dissociation on Si(100) Surface. Importance of Proton Tunneling

The dissociation of an ammonia molecule on a cluster of Si atoms simulating the 100 silicon crystal structure with two Si dimers has been investigated by means of the DFT and an approximate instanton methods. The model corresponds to the low coverage limit of the surface. Absolute rate constants of two different dissociation paths are evaluated together with deuterium isotope effects. It is demonstrated that, even at room temperatures, the process is dominated by tunneling and that dissociation to a silicon atom of the adjacent dimer, rather than a silicon within the same dimer, is the prevailing mechanism. This leads to creation of a metastable structure which will slowly decay through a two-step hydrogen atom migration towards the absolute minimum on the potential energy surface corresponding to the NH 2 group and the hydrogen atom residing in the same dimer.


Introduction
Silicon (100)-(2 x 1) surface is known for its high chemical reactivity resulting from the presence of dangling bonds that leads to formation of Si dimers on the surface [1].The interaction of ammonia with silicon surfaces has been studied experimentally at various temperatures [2][3][4][5][6][7][8][9].In particular, Avouris and coworkers studied extensively the (100) surface [3] using scanning tunneling microscopy (STM), X-ray (XPS) and ultraviolet photoemission (UPS) spectroscopies.They found that at temperatures as low as 90 K, dissociation of ammonia occurs and that it leaves the Si dimers essentially intact.Hill et al. [4] observed Si-NH x species on the surface by UPS and XPS spectroscopies.Dresser et al. [5] and Johnson et al. [6] have demonstrated the presence of NH 2 fragments on the Si surface and dissociative adsorption of NH 3 with the formation of Si-NH 2 and Si-H bonds.Dissociation of ammonia in the NH 3 → NH 2 + H reaction is thus the primary process to be studied.It involves transfer of a hydrogen atom from the chemisorbed ammonia to a silicon atom located presumably within the same Si dimer that adsorbed the ammonia molecule.These findings were supported by other research groups with a variety of experiments [7][8][9].Finally, some vibrational frequencies of modes involving adsorbed NH 3 , NH 2 and H fragments on the (100) surface were obtained experimentally by means of high resolution vibrational electron energy loss spectroscopy (HREELS) [7,11,12].
Recently, Srivastava and co-workers [10] investigated the dissociative adsorption of ammonia, phosphine and arsine molecules on the (100) surface by means of density functional theory (DFT).
They found that this process is strongly exothermic for all three molecules and that there is a substantial barrier (10-15 kcal/mol) to be overcome for the hydrogen atom transfer.On the other hand, the process of binding of an ammonia molecule to Si(100) was found to be barrierless.The authors studied only the case when the Si-NH 2 and Si-H bonds are within the same Si dimer.
Since the dissociation of ammonia on Si(100) surface is observed at very low temperatures (90 K), hydrogen atom tunneling must be the predominant mechanism under such conditions and, consequently, large deuterium effect in the low temperature dissociation must be displayed.Despite this, the tunneling mechanism was never addressed and no theoretical or experimental studies of such effect have been undertaken so far.The present paper fills this gap, in an attempt to assess the importance of hydrogen tunneling and the effect of deuterium substitution on the dissociation process.
We demonstrate that tunneling dominates the dissociation process even at room temperatures and that it leads to marked deuterium effects.DFT generalized gradient approximation (GGA) method, together with a proton-dymamics method based on the instanton formalism are used in our model calculations.

Methods
We study a cluster of 15 Si atoms representing an equivalent of two Si dimers on its '100' surface.
All other Si atoms of the cluster are saturated with hydrogen atoms fixed at their positions in a Si crystal.The structure of this cluster is depicted in Fig. 1a.An ammonia molecule binds to one of the four 'surface' silicon atoms (Si 2 ).Subsequently, it dissociates by donating one of the hydrogens to another silicon atom in the same (Si 1 ) or in a neighbouring dimer (Si 4 ).This model excludes multiple adsorptions of ammonia molecules in neighbouring dimers and thus is applicable only for low coverage of the silicon surface.We use Becke's hybrid three-parameter exchange functional [13] and Lee-Yang-Parr correlation functional [14,15] (B3LYP), together with the 6-31G* basis set.Gaussian 98 suite of programs was used in the calculations [16].The electronic-structure and force field output of these calculations were used to generate the multimensional potential energy surface that governs the dissociation dynamics.The latter was analyzed by the Approximate Instanton Method (AIM) [17], as implemented in the program DOIT 1.2 [18] and used to calculate proton-transfer rate constants.
Below we give a brief summary of the method; detailed description can be found in the recent reviews [17] and references therein.
AIM is a multidimensional semiclassical approach developed to describe tunneling phenomena.It has been applied successfully to a number of proton transfer processes in the ground and excited states of large polyatomic molecules, molecular clusters and biological models [19][20][21][22][23][24][25][26][27][28].The method is based on the concept of least action path, called instanton, that represents a compromise between energetically favorable but long paths, like the minimum energy path, and short but energetically costly paths.It avoids an explicit evaluation of the multidimensional instanton trajectory by extrapolating exact solutions for the instanton action based on two-and three-dimensional models [29].
The Hamiltonian is generated in terms of the normal modes of the transition state (x,y), the mode x with imaginary frequency ω* being the reaction coordinate and y representing 3N-7 (transverse) modes.The kinetic energy operator is taken to be diagonal.The motion along the reaction coordinate is described by a double-minimum potential and the transverse modes are treated as independent harmonic oscillators, coupled linearly to the reaction coordinate.The transverse modes that are active in the process belong to two distinct groups: modes that undergo effective displacement between the initial and final configurations, or "antisymmetric" modes, normally suppress tunneling; modes with no such displacement, or "symmetric" modes, enhance tunneling.Generally, every active mode gives both contributions unless tunneling occurs between two identical configurations [17].Since roomtemperature experiments are typical for the dissociation process at hand, two contributions to the rate must be included: a classical contribution for hydrogen atoms with enough energy for over-barrier transfer and a tunneling contribution for lower-energy atoms.In the DOIT 1.2 program [18], used for the present dynamics calculations, the classical rate constant is calculated by standard Transition State Theory (TST) and the tunneling rate constant by AIM: where Here 0 U is the barrier height corrected for zero-point energy (ZPE) changes between the reactant state and the transition state, 0 R Ω is the effective frequency of tunneling in the reactant state and S I (T) is the multidimensional instanton action (hereafter in units of ) for low-dimensional model systems [17,29] 0 ( ) ( ) ( ) 1 ( ) In this expression S I 0 (T) is the instanton action of an effective one-dimensional motion along the reaction coordinate which accounts for coupling to the "fast" modes; these modes are treated in the adiabatic approximation.Coupling between the reaction coordinate and the "slow" modes yields the correction terms δ α,s : these modes are treated in the sudden approximation.The "antisymmetric" modes (or antisymmetric components thereof) contribute a Franck-Condon factor s a a e α δ − ∑ to the transition and thus lower the rate of transfer; this effect is similar to that of a condensed phase environment.The δ s parameters, which refer to "symmetric" modes (or components thereof), shorten the tunneling distance and lower the barrier, which facilitates transfer.The factor α s ≤ 1 describes the modulation of the antisymmetric correction by the symmetric coupling.The division of the transverse modes into "fast" vibrations treated in the adiabatic approximation and "slow" vibrations treated in the sudden approximation is based on a specific parameter, called zeta-factor, which depends on their frequencies and on the strength of their couplings to the tunneling mode [29].
All parameters needed for the dynamics can be easily obtained from the energetics and Hessians of the transition state and minima, which the DOIT code reads directly from standard quantum-chemical output.This allows evaluation of both, the tunneling and classical components of the rate constant in Eq. 1 and thus assessment of the relative importance of the two mechanisms.Tunneling is usually associated with large deuterium kinetic isotope effects (KIEs), measured as the ratio of the rate constants for H and D transfer.In the absence of tunneling, the KIE is governed by the effect of deuterium substitution on the difference in zero-point energy correction If transfer occurs by tunneling that is coupled to the skeletal vibrations, the KIEs exhibit unusual features.Deuterium substitution increases the tunneling path, but the main increase of 0 I S is due to the doubling of the mass of the tunneling particle, as 0 I S is proportional to m .This effect is mitigated by coupling of the tunneling motion to skeletal modes that narrow and lower the barrier ("symmetric" modes).As a result, the rate constant no longer depends exponentially on m , as in the onedimensional case, but exhibits a weaker mass-dependence, as the transition probability is affected by slow motions which are isotope-independent.This effect of reduced KIEs in multimensional tunneling is general, depends only on the strength of coupling and will always be present in complex systems.
Thus hydrogen-atom tunneling in such systems may manifest itself by two seemingly contradictory signatures: tunneling rates that are higher and KIEs that are similar to, or even lower than their classical counterparts.This effect was found in numerous application of the AIM to proton or H-atom tunneling [17,[19][20][21][22][23][24][25][26][27][28], most notably in a recent study of a model of the enzyme carbonic anhydrase II [28], where three protons move concertedly along a bridge formed by water molecules and the observed KIE of 3-4 at room temperature appears to be anomalously weak in face of the accumulating ZPE corrections.

Results: Structures, Energetics and Vibrational Frequencies
Three structures corresponding to hydrogen-atom transfer within a given dimer (N → Si 1 ) were optimized with the B3LYP /6-31G* method.These are NH 3 on Si(100) surface (Fig. 1a); NH 2 and H attached to Si atoms of the same surface dimer (Fig. 1b); and the transition state between these two (Fig. 1c).The hydrogen atoms which saturate the dangling bonds of the Si atoms and replace the inner layers of the crystal were frozen during optimizations in the positions of those atoms.Table 1 lists the energies and some of the structural parameters connected with the proton tranfer.In addition, the energy of separated Si(100) cluster and NH 3 molecule is given.One sees in Fig. 1b that the NH 2 configuration is of a gauche type and not a symmetric one, as assumed in many studies [30].Our calculations predict rather large exothermicity of the dissociative adsorption of ammonia on Si (100) surface.This effect (23.6 kcal/mol) is smaller than that obtained in some other theoretical treatments [10,33] (31.5 and 42.0 kcal/mol, respectively) but somewhat larger than those obtained in Refs.[34,35] (18 and 22 kcal/mol).The calculated barrier to the transfer is rather large, 20.9 kcal/mol, as compared with the values of 15 kcal/mol [10,33], 14 kcal/mol [34] and 24 kcal/mol [35].This indicates that our cluster model reproduces rather well results of much larger computational efforts.
The energy of the NH 3 -Si(100) complex is 35.85 kcal/mol lower than that of separated ammonia molecule and Si(100) cluster; no correction for basis set superposition error is included in this value.This calculated adsorption energy is close to that obtained by Fattal et al. [33] (33 kcal/mol) and somewhat larger than 29 kcal/mol, 28 kcal/mol and 25 kcal/mol obtained in Refs.[35], [34], and [10], respectively.
There is an alternative route for dissociation that is encompassed by our model.In this case the dissociating hydrogen moves to a silicon atom of an adjacent dimer (Si 4 ).The structures of the resulting NH 2 + H Si cluster and of the corresponding transition state are depicted in Fig. 2. The relevant energies and structural parameters for this alternative pathway are also included in Table 1.
The possibility of dissociation of NH 3 in Si(100) surface with the hydrogen atom transferred to the nearest neigbour Si dimer (N → Si 4 ) has not been considered theoretically before.Our results indicate that this route can contribute to the process as well.The structure of the product of this path is depicted in Fig. 2a, and the transition state structure is given in Fig. 2b.The energies and some geometrical parameters of these structures are collected in Table 1.On the other hand, the much smaller exothermicity of this reaction (13 kcal/mol vs. 23 kcal/mol) will tend to slow this process, whereas we also note that the Si 2 SiSi 4 angle at the transition state configuration closes to 99 degrees from the initial value of 106 degrees.We will discuss the rates of dissociation for the two paths, called intraand inter-dimer paths, as functions of temperature in the following section.
Table 2 collects some of the calculated, unscaled frequencies and compares them with other theoretical results and experimental data for NH 3 and ND 3 species adsorbed on the Si(100) surface; all NH stretch, NH scissor and NH wag frequencies obtained in our calculations are listed.There is an overall good agreement with other theoretical and experimental results.A large discrepancy was found between the calculated frequencies of NH scissoring obtained in Ref. [10].The experimental value is in a better agreement with our frequency, especially if we remember that unscaled calculated values are given in Table 2.The other large discrepancy occurs for the calculated values of the Si-N stretch of NH 3 on the (100) surface: our value of 479 cm -1 is considerably smaller than the 630 cm -1 value obtained by Srivastava et al. [10].Unfortunately, there is no experimental data for this particular vibration.The transition state structure of Fig. 1c is characterized by an imaginary frequency of 1472i cm -1 .Ref. [21] 4 exp.Ref. [7] 5 exp.Ref. [12] We find that the dissociation of the second hydrogen atom, which leaves a NH group, is a highly endoenergetic process.The calculated energy of the NH + H + H is 15.7 kcal/mol above that of the undissociated NH 3 on the same cluster, which makes the second step of dissociation highly improbable.

Rates of Dissociation and the Isotope Effect
In order to assess the importance of tunneling, we used the Approximation Instaton Method (AIM) and the DOIT 1.2 code [18] to calculate the rate constants of NH 3 and ND 3 dissociation on the model silicon cluster.Below we present the rate constants and the KIE's for the two dissociation pathways outlined in Section 3.

Intra-Dimer Pathway
If the dissociation follows this path, the reactant is the structure presented in Fig. 1a, the product is the structure in Fig. 1b, and the transition state that in Fig. 1c.In the transition state configuration, the distance between the transferred atom and the donor (N) and acceptor (Si) is 1.50 Å and 1.80 Å, respectively.As this process has large exothermicity, the tunneling barrier is strongly asymmetric, with a "late" transition state, its width along the reaction coordinate being ∆x = 1.141/0.725Å amu ½ on the reactant/product side, respectively.
The reaction coordinate, which is the mode with imaginary frequency, is purely hydrogenic, as illustrated by the value of ω* H/D = 1472i/1061i cm -1 and Fig. 3.The coupling to the surface motions is dominated by "antisymmetric" modes, the strongest being that to a mode with frequency of about 211 cm -1 ; in this case the motion of the donor and acceptor of the transferred atom leads to asymmetry of the tunneling potential and thus to reduced transfer probability.The coupling to "symmetric" modes of this type is to one with frequency of about 432 cm -1 ; in this case the motion of the donor and acceptor shortens the tunneling distance, thus enhancing tunneling.This effect results in reductuion of the donor-acceptor distance which is shortened from 3.75 Å in the reactant to 3.17 Å in the transition-state configuration.These two modes are also illustrated in Fig. 3.
Table 3 lists the calculated rate constants at three different temperatures (100 K, 200 K and 300 K) for the dissociation of NH 3 and ND 3 molecules on the Si(100) model cluster.Overall multidimensional rates, as defined by Eqs.(1-3); one-dimensional rates, in which the participation of transverse modes is neglected; and classical over-the-barrier rates (2) are listed.We see that tunneling dominates the dissociation process at these temperatures, since the classical components in Eq. ( 1) are smaller than Table 3. Rate constants (in s -1 ) and kinetic isotope effects (KIE) k(NH 3 )/k(ND 3 ) for the two pathways of dissociation of the NH 3 /Si(100  their tunneling counterparts even at room temperature.Also the deuterium effect is quite large, as the barrier is high and wide.The KIE is a combination of primary and secondary effects, the primary part taking over with increasing temperature, as illustrated by the fact that the rate of dissociation of NDH 2 , where only the transferred hydrogen has been deuterated, is very close in value to that of ND 3 .

Inter-Dimer Pathway
This path yields an (intermediate) product illustrated in Fig. 2a and is characterized by the transition state depicted in Fig. 2b.The motion along the reaction coordinate is similar to the one found above, with ω* = 1361i/990i cm -1 ; the tunneling barrier again is strongly asymmetric, but with an "early" transition state, as ∆x = 0.799/1.179Å amu ½ on the reactant/product side, respectively.The coupling picture is also very different, as the silicon atom accepting the hydrogen now belongs to a different dimer.In contrast to the intra-dimer pathway, the coupling is dominated by "symmetric" modes, which in this case is much stronger than the suppressing effect, the latter being altogether weaker and less localized.The strongest-coupled modes of both types are of similar frequencies to those found earlier, but involve different motions, as illustrated in Fig. 4. Similar to the intra-dimer pathway is the effect of shortening of the donor-acceptor distance, which is reduced from 3.45 Å in the reactant to 2.91 Å in the transition state.The rate constants of dissociation of NH 3 and ND 3 are also listed in Table 3.We could not determine with sufficient accuracy the tunneling rates of ND 3 dissociation, as the symmetric coupling in this case is strong and the AIM results are unreliable; the values of these rates and the corresponding KIEs, given in Table 3 for convenience, are thus only qualitative.
The KIEs of the classical mechanism given by Eq. ( 5) are defined by the barrier heights for H/D, corrected for zero-point energy.For the two pathways these equal 17.9/16.5kcal/mol and 17.6/16.3kcal/mol, respectively, and at room temperature yield similar KIEs of the order of 10.If there were no coupling to the surface motions, the rates would be given by the one-dimensional rate constants in Table 3 and the corresponding KIEs would have been of the same order.This value is reduced by the symmetric coupling, especially for the inter-dimer pathway, where this coupling is significant.And although in this case the calculated reduction of the KIEs is overestimated by the AIM, qualitatively the effect is there, as illustrated by the data in Table 3.

Discussion
We have reported the kinetics of dissociation of a single ammonia molecule on a model cluster representing the Si(100) surface.Two silicon dimers with dangling chemical bonds are included in the model, which yields structures that are in agreement with other calculations treating the Si crystal with the slab model.Also the energetics -adsorption energy, exotermicity and barriers for the dissociation process are in a close agreement with the previous theoretical treatments.
The completely new result of this study is the first principles calculation of the rate constants of dissociation of NH 3 and ND 3 molecules adsorbed on the Si(100) surface.This calculation takes fully into account the quantum mechanical nature of the process of tunneling of a light particle such as H or D atom through a potential barrier.Two different paths of dissociation were studied.The conventional path, tacitly assumed so far as the only one possible, is the one where the hydrogen atom is transferred to a Si atom in the same dimer in which the original NH 3 molecule was adsorbed; as a result the NH 2 group and the H atom reside on the same dimer.Such a structure corresponds to the absolute minimum on the potential energy surface of the model cluster.The calculated exothermicity of this process is 23.6 kcal/mol, which means that the total energy released by adsorption of an ammonia molecule and its subsequent dissociation is ca 59 kcal/mol.We have shown that the process is dominated by the tunneling mechanism even at room temperature, the calculated rate of ca 5 s -1 being in agreement with the experimental observation that all NH 3 molecules are dissociated at that temperature when coverage of the Si(100) surface with ammonia is low.We predict that replacement of NH 3 with ND 3 will slow down the process of dissociation by an order of magnitude at 300 K.The second path is the one in which the hydrogen atom is transferred to the nearest Si atom located in the adjacent dimer.This leads to a metastable state, some 10 kcal/mol above the absolute minimum, which is separated from the latter by a large energy barrier of ca 50 kcal/mol and thus should be quite stable at low temperatures.
We calculate that, at room temperature, this alternative path of dissociation is about an order of magnitude faster than the conventional one.Thus, for low ammonia coverage of the Si(100) surface, one should observe Si-H ans Si-N bonds on adjacent dimers in the same row.If, at higher coverage, ammonia molecules are adsorbed on silicones in the same row, the second path is blocked and dissociation within dimers is the only possible channel.However, if ammonia molecules alternate between rows of subsequent Si dimers, then both paths for dissociation remain open at higher coverage.
We found also that both mechanisms of dissociation are strongly hindered at very low temperatures by surface vibrations in which Si atoms swing in the same direction either along the dimeric bonds [21] (211 cm -1 , first path) or along rows of atoms (202 cm -1 , second path).However, as temperature rises, a second vibration takes over which brings closer the donor (N) and acceptor (Si) atoms, and thus promotes the transfer, leading to a sharp increase of the rate of dissociation.This vibration (432 cm -1 , first path; 417 cm -1 , second path) is far more active for the inter-dimer mechanism, leading to the domination of this process at higher temperatures.The role of the coupling on the thermal excitations of the process, as well as the typical for tunneling no-Arrhenius behaviour of the rate constant is illustrated in Fig. 5.
Let us note in passing that our results predict very low rates of tunneling at 100 K for both paths.
Thus, formation of a Si-NH 2 bond after ammonia deposition at 90 K should be a very rare event.The fact that it is still observed indicates that the energy released during ammonia molecule adsorption with the formation of a Si-NH 3 bond (35.85 kcal/mol) may leave some bound NH 3 molecules with vibrationally excited N-H stretching modes.An estimate shows that excitation of one vibrational quantum of the N-H stretch of a bond undergoing subsequent dissociation would increase the rate at 100 K by ca 6 orders of magnitude, making it plausible that some of the adsorbed molecules may dissociate from vibrationally excited states reached by vibrational relaxation following the adsoption process.Finally, we should mention that the tunneling process is very localized and that it is affected mainly by the surface atoms.Thus, replacing the fixed hydrogen atoms of the cluster by particles with the mass of silicon atom changes the calculated rate constants less than twice, indicating again that the present cluster model forms a reliable approximation for localized processes taking place on the Si(100) surface.
model.This buckling diminishes again upon NH 3 dissociation, resulting in the binding of the hydrogen atom to another Si atom in the suface dimer.The Si -Si and N -Si distances are close to those calculated by other authors (see ref.[10] and references therein) and reproduce very well the marked shortening of the N -Si distance upon dissociation.

Figure 1 .
Figure 1.Calculated structures of ammonia adsorbed on a cluster-simulated (100) surface of silicon.Two different projections are given: a) adsorbed NH 3 , b) dissociated NH 3 with the H atom on the other Si atom of the silicon dimer, c) transition state for the dissociation process.

Figure 2 .
Figure 2. The calculated structure of: a) dissociated NH 3 with the H atom on the next silicon dimer, b) the corresponding transition state.

Figure 3 .
Figure 3. Normal modes most important for the hydrogen transfer within a silicon dimer.211 cm -1the most active hindering mode, 432 cm -1 -the most active helping mode, 1472i cm -1 -the reaction mode.

Figure 4 .
Figure 4. Normal modes most important for the hydrogen transfer between adjacen dimers.202 cm -1the most active hindering mode, 417 cm -1 -the most active helping mode, 1361i cm -1 -the reaction mode.

Figure 5 .
Figure 5. Calculated overall rate constants (thick lines) for the dissociation of a single ammonia molecule on the model Si cluster, for the intra-dimer (solid lines) and inter-dimer (broken lines) pathways.The thin lines represent the contribution of the classical over-the-barrier mechanism.
The calculated Si-NH 2 distance of 1.75 Å is in a very good agreement with the value of ca.1.73 Å obtained by Franco et al. by the scanned-energy mode photo-electron diffraction technique [31,32].The structure of the Si-NH 2 group is pyramidal with the dihedral angle SiNHH being -137.2degrees.Chemical binding of an ammonia molecule causes attenuation of the magnitude of bunckling of adjacent Si dimers, as defined by the dihedral angle Si 1 Si 2 Si 3 Si 4 , calculated to be -31.1 degrees for our

Table 1 .
Geometrical parameters (Å and degrees) and energies of the clusters presented in Figs.1 and 2.

Table 2 .
Selected calculated and experimental frequencies (cm -1 ) of the adsorbed NH 3 and it dissociative product on the Si(100) surface.