# Electric field dependent structural and vibrational properties of the Si(100)-H(2×1) surface and its implications for STM induced hydrogen desorption

## I INTRODUCTION

STM induced desorption of hydrogen(H) from the monohydride Si(100) surface offers the possibility of lithography with atomic resolution. Investigations of the desorption mechanism have established the dependence of the desorption rate on the bias voltage, tunnel current and H isotope. At high positive biases, $`V_b>4`$ V, the experimental results are consistent with electron induced desorption due to direct excitation of the Si-H bond by a single electron. At negative and low positive biases the desorption rates show power-law dependencies on the electron current consistent with a multi-electron process, and the measured desorption rates are in quantitative agreement with first-principle calculations.

The desorption by multi-electron scattering is only possible because the H stretch frequency has a long lifetime. The long lifetime is a result of a vibrational quantum too low for coupling with electron-hole excitations, while well above the Si phonon spectrum, and can thus only transfer energy to the substrate via a multi-phonon process. At room temperature experimental estimates of the lifetime due to this process are $`\tau 10`$ ns. However, a local excitation is not an eigenmode and will decay into a H surface phonon by a coherent process. This decay is several orders of magnitude faster than multi-phonon energy relaxation, and must therefore be included in the theoretical models. It has been proposed by Persson and Avouris that the vibrational Stark shifts due to the electric field from the tip can localize vibrational modes below the tip. The localized modes may still transfer energy laterally by incoherent diffusion (the Förster mechanism) but it was found that this decay channel is also reduced by the Stark shifts. However, the work assumed that the electric field from the STM tip is localized on a single H atom below the tip, and this is not the case for realistic tip geometries.

In this work we present results for the vibrational properties of the Si(100)-H(2$`\times `$1) surface in the presence of the electric field from a more realistic model for the STM tip. The tip is described by sphere of radius, $`R_t=500`$ Å, with a protrusion of atomic dimensions, and we determine the electric field by solving the Poisson’s equation numerically. To obtain the effect on the H vibrations we set up a phonon Hamiltonian with parameters obtained from a first principles calculation of the vibrational properties of the Si(100)-H(2$`\times `$1) surface in an external electric field. We find that the electric field does give rise to a localized vibrational state below the tip; however, its lifetime is very short (10 ps) due to incoherent exciton motion. However, we find that the anharmonicity of the Si-H bond potential reduces the lateral energy transfer of higher excited excitations($`n>1`$) of the Si-H bond. We present first principle calculations of the desorption rate taking this effect into account, and find that two phonon excitations play an important role in the desorption process.

The organization of the paper is the following. In Section II we describe the first principles method which in section II A is used to calculate the zero field atomic structure and Si-H stretch frequencies of the Si(100)-H(2$`\times `$1) surface. The electric field dependence of the frequencies is calculated in Section II B. In Section III we introduce a simple dipole-dipole interaction model for the Si-H stretch phonon band, and uses it to find localized vibrational states in the presence of electric fields from different STM tip geometries. In Section IV we calculate the lifetimes of the localized states due to incoherent lateral diffusion. In section V the lifetimes are used to model STM induced desorption. Section VI summaries the results.

## II Structure and vibrational properties of the Si(100)-H(2$`\times `$1) surface.

In this section we calculate the vibrational frequencies and the dipole-dipole coupling matrix elements of the H vibrations on the Si(100)-H(2$`\times `$1) surface. In subsection A we present calculations for the unperturbed Si(100)-H(2$`\times `$1) surface, and the vibrational and structural shifts due to an external planar field are calculated in subsection B.

The first principles calculations are based on density functional theory within the Generalized Gradient Approximation (GGA) for the exchange-correlation energy. Since we only consider filled shell systems, the calculations are all non-spin-polarized. Ultra-soft pseudo potentials constructed from a scalar-relativistic all-electron calculation are used to describe H and silicon(Si). The wave functions are expanded in a plane-wave basis set with a kinetic-energy cutoff of 20 Ry, and with this choice absolute energies are converged better than $`0.5`$ mRy/atom.

With this approach we find a Si lattice constant of 5.47(5.43)Å and bulk modulus of 0.89(0.97) Mbar (parentheses show experimental values). For the H<sub>2</sub> molecule we obtain a bond length of 0.754(0.741)Å, a binding energy (including zero point motion) of 4.22(4.52) eV and a vibrational frequency of 4404(4399)cm<sup>-1</sup>. Generally the comparison with experiment is excellent, and similar theoretical values have been found in other studies using the GGA.

### A  Zero field properties

To model the monohydride Si(100) surface at zero field we use a (2$`\times `$1) slab with 12 Si atoms and 6 H atoms. The atoms at the bottom surface are bulk like, and their dangling bonds are saturated with H atoms. The two surfaces are separated by a 7.5 Å vacuum region, and we use the dipole correction in order to describe the different workfunction of the two surfaces. The surface is insulating and we use 2 $`k`$-points in the irreducible part of the Brillouin-zone (BZ) for the BZ integrations. Test calculations with more dense meshes show that BZ integration errors are negligible small.

Figure 1 shows the atomic structure after relaxation of the H atoms and the 4 upper Si layers. Positions of the H atoms and the two upper Si layers compare well with other studies. For the third and fourth layer Si atoms we find a small asymmetric relaxation, and to our knowledge such relaxations have not been included in previous studies.

To obtain the dynamical matrix of the Si-H stretch frequency we make H displacements of $`\pm 0.07,0.14,\mathrm{},0.35`$ Å in the Si-H bond direction and fit a sixth order polynomial to the data points. Since the Si-H stretch frequencies are four times higher than Si bulk frequencies, the Si substrate acts like a solid wall, and we therefore use the H mass $`\mathrm{M}_\mathrm{H}`$ when calculating frequencies from the dynamical matrix. With this approach we have calculated the frequency of the symmetric stretch $`\omega _\mathrm{s}`$ and the asymmetric stretch $`\omega _\mathrm{a}`$ at three high symmetry points in the surface Brillouin zone. In Table 1 the results are listed together with $`\mathrm{\Gamma }`$ point frequencies(parenthesis) obtained by infrared spectroscopy and the comparison between theory and experiment is excellent, especially we note that theory correctly predicts the splitting, $`\omega _\mathrm{s}\omega _\mathrm{a}=11`$ cm<sup>-1</sup>.

We next investigate the bonding potential of the H atoms. In Fig. 2 the solid circles show the hydrogen energy, $`E_\mathrm{H}`$, for the Si-H bond lengths, $`z`$, used in the calculation of the dynamical matrix. The data are accurately described by a Morse potential

$$E_\mathrm{H}(z)=E_\mathrm{d}(e^{2\alpha (zz_0)}2e^{\alpha (zz_0)}),$$
(1)   
and from a least-squares fit we obtain the frequency $`\omega _0=0.26`$ eV($`\alpha =1.57`$ Å<sup>-1</sup>), equilibrium bond length $`z_0=1.50`$ Å, and desorption barrier $`E_\mathrm{d}=3.4`$ eV. The extrapolated desorption energy coincides with the surface energy without a H atom plus the energy of a spin-polarized H atom. The triangles in Fig. 2 show the total energy for large values of $`zz_0`$. When the interaction between the H atom and the surface becomes weak the electrons start to spin-polarize and the data points in Fig. 2 show that spin polarization effects become important for $`z>2.5`$ Å.

The inset shows the change in the surface dipole, $`\mathrm{\Delta }p`$ as a function of $`z`$ (the positive direction is from Si to H). The surface dipole increases almost linearly with $`z`$ and the dynamic dipole moment is $`\gamma =p/z=0.6`$ Debye/Å (=0.13 e). Modelling the surface dipole by an effective charge $`e^{}`$ on the H atom and its image charge -$`e^{}`$, we find $`e^{}0.07e`$. The sign of this charge transfer from Si to H is in accordance with a higher electronegativity of H relative to Si.

### B Electric field dependent properties

To model the surface in a planar external field we use a (2$`\times `$1) slab with 24 Si atoms, 2 H atoms, and a vacuum region of 10 Å and the external field is modelled using the method of Ref. . The Si atoms at the back side of the slab are not passivated by H atoms, and dangling bonds on these atoms can donate free electrons and holes. In this way we take into account the effect of mobile carriers. Other computational details are identical to those for the zero field calculation.

Curves in Fig. 3 show the field dependence of the equilibrium Si-H bond length $`z`$, the $`\mathrm{\Gamma }`$ point symmetric Si-H stretch frequency $`\omega _\mathrm{s}`$ and the $`\mathrm{\Gamma }`$ point symmetric–asymmetric splitting $`\mathrm{\Delta }\omega =\omega _\mathrm{s}\omega _\mathrm{a}`$. We first notice that all three quantities have an extremum at $`1.5`$ V/Å. This behaviour can be described by a simple Si-H tight-binding model with a field dependent H on-site element. The extrema occurs at the field where the H and Si on-site levels are in resonance, since at resonance the Si-H bond is strongest, and therefore the bond length minimal and the vibrational frequency maximal. Furthermore, at resonance the H dynamic dipole vanishes, and therefore also the part of $`\mathrm{\Delta }\omega `$ caused by H-H dipole interactions.

The three solid lines in Fig. 3 show second order polynomials obtained by least squares fit to the data. The interpolated zero field values of $`z_0`$ and $`\omega _\mathrm{s}`$ agrees exactly with those obtained in section II A, while the interpolated $`\mathrm{\Delta }\omega `$ zero field value is slightly off. Taking into account the quite different slabs used for the two calculations we find the agreement fully satisfactory, and note that the difference can be taken as a measure of the accuracy of the approach.

Recently, the electric field dependent properties of the H/Si(111)(1$`\times `$1) surface were calculated by Akpati et al., and they found Stark shifts $`30`$ percent larger than in the present calculation, and the extremum in bond length and frequency appears for a field of 1 V/Å. The agreement with the present calculation seems reasonable, bearing in mind that the Stark shifts are for different crystallographic surface directions. However, part of the difference might be due to the use of a cluster geometry and the local spin density approximation in Ref. . The present study is based on a slab geometry and the GGA. We expect that the thick slab geometry better describes the electric field induced polarization of the surface.

## III STM induced Stark localization

In this section we will model the collective modes of the Si-H stretch vibrations by a set of local oscillators interacting through dipole forces, and use this model to calculate the Stark localization in the external electric field from a STM tip. In Fig. 4 is shown the lattice sites of the oscillators, corresponding to the positions of the H atoms in the (2$`\times `$1) cell. Each oscillator is described by a local frequency $`\omega _i`$ and a dynamic dipole moment $`\gamma _i`$. The Hamiltonian of the system is given by

$$H_{ij}=\mathrm{}\omega _i\delta _{ij}+\frac{\chi _i\chi _j}{|𝐫_i𝐫_j|^3}(1\delta _{ij}),$$
(2)   
where $`r_i`$ is the position of oscillator $`i`$, and $`\chi _i=\sqrt{\mathrm{}/2\mathrm{M}_\mathrm{H}\omega _i}\gamma _i`$.

We first consider the zero field case of identical oscillators with parameters $`\gamma _0`$, $`\omega _0`$. We find the dispersion from a numerical Fourier transform, and at the $`\mathrm{\Gamma }`$ point we obtain the Hamiltonian

$$H(\mathrm{\Gamma })=\left(\begin{array}{cc}\omega _0+4.05V_0& 5.12V_0\\ 5.12V_0& \omega _0+4.05V_0\end{array}\right),$$
(3)   
where $`V_0=(\gamma _0^2)/(\mathrm{M}_\mathrm{H}\omega _0a_{100}^3)`$ and $`a_{100}=3.87`$ Å is the surface lattice constant. The two eigenmodes are $`\omega _\mathrm{s}(\mathrm{\Gamma })=\omega _0+9.17V_0`$, and $`\omega _\mathrm{a}(\mathrm{\Gamma })=\omega _01.07V_0`$. Using the calculated values of $`\gamma _0`$ and $`\omega _0`$ from section 2A, we obtain $`V_0=0.07`$ meV and thereby $`\mathrm{\Delta }\omega =5`$ cm<sup>-1</sup>. Dipole-dipole interactions can therefore only account for half of the dispersion obtained in the the frozen phonon calculation. We suggest that the remainder of the splitting is due to a short-range electronic interaction. This electronic interaction gives rise to an additional splitting, and explains why $`\mathrm{\Delta }\omega >0`$ at $`E=1.5`$ V/Å(see Fig. 3c) even though the dipole-dipole interaction vanishes at this field.

To simplify the calculations we will in the following use the dipole-dipole interaction model(Eq. (2)) to describe all the interactions, and determine field dependent parameters $`\omega _\mathrm{E}`$ and $`\gamma _\mathrm{E}`$ by relating the $`\mathrm{\Gamma }`$-point eigenmodes of the model to the calculated frozen phonon values. In this way we approximate the effect of the short-range electronic interactions by long-range dipole forces. To test the accuracy of this approximation we have used the model to calculate $`\omega _\mathrm{s}`$ and $`\omega _\mathrm{a}`$ at the J and J’ point in the surface Brillouin zone. The result is shown in Table 1 and the comparison with the first principles calculation is reasonable.

We next model the electric field below the STM tip. Usually it is found that the tip has a curvature in the range $`100`$ Å – $`1000`$ Å and it is generally accepted that the atomic resolution arises from a small protrusion or a single atom sticking out of the tip. We use the geometry in Fig. 5 to model such a tip. Two parameters, the tip curvature, $`R_t`$, and the protrusion size $`s`$, determine the tip geometry and we present results for parameters in the range, $`R_t=100`$$`500`$ Å and $`s=0`$–9 Å. For the tip-sample distance we use, $`h=3`$$`7`$ Å, which is the typical distance range in STM lithography experiments. To find the electric field below the tip the Poisson’s equation is solved numerically using ANSYS finite element analysis. Curves in Fig. 6 show the radial electric field at the surface for a potential difference of 5 V between the tip and the surface. Curves in Fig. 6a show the result when there is no protrusion on the tip ($`s=0`$ Å), and for this geometry the electric field attains its half value at $`r\sqrt{Rh}`$. In Fig. 6b results are for a tip with a protrusion of size $`s`$, and the protrusion gives rise to a reduced electric field below the tip and it decays rapidly around the tip apex. The small protrusion changes the electric field of the tip very little, and localization of the electric field is most pronounced for the large protrusion. Curves in Fig. 6c show the electric field from the geometry with $`R_t=500`$ Å and $`s=6`$ Å at three tip-surface separations $`h=3,6`$, and $`9`$ Å. The curves show that the field becomes more localized when the tip approaches the surface.

To determine the vibrational states below the tip in the presence of the electric field, we set up the Hamiltonian in Eq. (2) for a finite cluster including sites up to a cutoff radius $`r_{cut}`$ and diagonalizes it numerically to find the eigenmodes $`\psi _\alpha `$ and frequencies $`\omega _\alpha `$. There may be several localized modes, but we are only interested in the localized state with the largest projection $`p`$ at the site directly below the STM tip ($`r=0`$). This state is determined using

$$p=\underset{\alpha }{\mathrm{max}}[|\psi _\alpha |0|^2],$$
(4)   
where the maximum is over the eigenmodes with frequency outside the phonon band, $`\omega _\alpha \omega _0[2.8V_0,9.2V_0]`$. For the spatial electric fields considered in this paper the value of $`p`$ is converged for cluster sizes $`r_{cut}=50`$$`100`$ Å.

Curves in Fig. 7 show $`p`$ and the corresponding vibrational frequency $`\omega _p`$ when the surface is subject to the fields of Fig. 6a. The “local E” curve corresponds to the geometry of Persson and Avouris where the electric field is localized at $`r=0`$. In this case a localized state is split of the phonon band at all negative fields, while at positive bias a threshold field of $`0.12`$ V/Å is needed to obtain localization. For typical fields in H desorption experiments, 0.5–1 V/Å, the state is completely localized at the site below the tip ($`p=1`$) and $`\omega _p`$ is similar to the frequency, $`\omega _\mathrm{E}`$, of the local oscillator at $`r=0`$. In the case of a tip with radius, $`R_t`$, a localized state exists for nearly all fields, i.e. at positive bias the threshold field is 0.03 V/Å. The lower positive threshold field compared to the “local E” case is obtained because the mode is a superposition of several sites with $`\omega _i\omega _\mathrm{E}`$. For typical fields in desorption experiments the mode has a substantial weight, $`p0.3`$, at the site below the tip.

In Fig. 8 we show the effect of a small protrusion on the STM tip. In this case the spatial localization is improved, and for fields 0.5–1 V/Å we have $`p0.8`$. Thus we confirm the results of Persson and Avouris, that there exists a localized mode in the region below the tip, however, it is not completely localized at a single site.

## IV Decay of the localized vibration

Consider an STM experiment where a tunneling electron scatters inelastically with the H atom below the tip and the H atom is excited into the $`n=1`$ vibrational state of the stretch mode. We now consider the decay of such an excitation. There are three important time scales, the coherent transfer time, $`\tau _c`$, the phase relaxation time $`\tau _{ph}`$ and the energy relaxation time $`\tau _{en}`$. The coherent transfer time is the time it takes for the local excitation to be transfered into the localized eigenmode below the tip, $`\tau _c\mathrm{}/\mathrm{\Delta }\omega 0.5`$ ps. Next the eigenmode looses its phase due to coupling with a $`200`$ cm<sup>-1</sup> Si phonon and the phase relaxation time has been measured to be $`\tau _{ph}8`$ ps at room temperature, and $`\tau _{ph}75`$ ps at 100 K. Finally the energy of the mode will decay into the Si substrate via a coupling with three Si-H bending modes(600 cm<sup>-1</sup>) and one 300 cm<sup>-1</sup> Si phonon. The time scale for this process is $`\tau _{en}10`$ ns at room temperature.

In the previous section we found a localized eigenmode with $`p0.8`$. The excitation at $`r=0`$ will be a superposition of this mode and more extended states. After $`\tau _c`$ the extended states have diffused away, thus $`80`$% of the initial excitation is in the localized eigenmode, and the total probability of finding the initial excitation at $`r=0`$ is $`p^20.6`$. For $`t>\tau _{ph}`$ the excitation can diffuse away to the neighbouring H atoms due to dipole-dipole couplings. This is the so-called Försters mechanism for incoherent diffusion, and in the following we will calculate the incoherent diffusion rate, $`w`$, using Försters formula

$`w`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{i0}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}|H_{0i}|^2A_0^0(\omega )A_i^0(\omega )𝑑\omega .`$ (5)   
In this equation $`A_i^0`$ is the spectral function at site $`i`$ for noninteracting H modes ($`H_{ij}^0\delta _{ij}`$), but including the coupling with substrate phonons which gives rise to the phase relaxation. The spectral functions are obtained from the noninteracting retarded Greens functions

$`G_i^0(t)`$ $`=`$ $`i\mathrm{\Theta }(t)[\widehat{c}_i(t),\widehat{c}_i^{}(0)],`$ (6)   
$`A_i^0(\omega )`$ $`=`$ $`2\mathrm{I}\mathrm{m}G_i^0(\omega ),`$ (7)   
where $`\widehat{c}_i^{}`$ and $`\widehat{c}_i`$ are local creation and annihilation operators of the stretch mode. The phase relaxation can be described approximately by the Hamiltonian

$$H_{ii}^0=\mathrm{}(\omega _i+\delta \omega \widehat{n}_i)\widehat{c}_i^{}\widehat{c}_i,$$
(8)   
where $`n_i=\widehat{b}_i^{}\widehat{b}_i`$ is the projected occupation operator of the $`\mathrm{\Omega }=200`$ cm<sup>-1</sup> Si phonon, and $`\delta \omega `$ is the change in the local frequency when the Si phonon is excited from level $`n`$ to $`n+1`$. The correlation functions of $`n_i`$ have been calculated by Persson et al.

$`n(t)`$ $`=`$ $`n_\mathrm{B}(\mathrm{\Omega }),`$ (9)   
$`n(t)n(0)`$ $`=`$ $`n(1+n)e^{\eta t}+n^2.`$ (10)   
The friction parameter $`\eta `$ describes the damping of the Si phonon, and $`n_\mathrm{B}(\omega )=(e^{\beta \omega }1)^1`$ is the Bose occupation number and $`\beta `$ the inverse temperature.

We now use the Matsubara formalism to obtain $`G_i^0`$ from an perturbation expansion in $`\delta \omega `$. We only consider the two lowest order diagrams shown in Fig. 9, and the corresponding self energies are

$`\mathrm{\Sigma }_i^{(1)}`$ $`=`$ $`\mathrm{}\delta \omega n,`$ (11)   
$`\mathrm{\Sigma }_i^{(2)}`$ $`=`$ $`\mathrm{}\delta \omega ^2n(1+n){\displaystyle \frac{1}{\omega \omega _i+i\eta }}.`$ (12)   
The $`\mathrm{\Sigma }^{(1)}`$ term gives rise to a small frequency shift, while the $`\mathrm{\Sigma }^{(2)}`$ leads to a damping of the mode. Considering only the latter term, we find

$`A_i^0(\omega )`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }_i(\omega )/\mathrm{}}{(\omega \omega _i)^2+\mathrm{\Gamma }_i(\omega )^2/4}},`$ (13)   
$`\mathrm{\Gamma }_i(\omega )`$ $`=`$ $`{\displaystyle \frac{2\delta \omega ^2}{\eta }}{\displaystyle \frac{n_\mathrm{B}(\mathrm{\Omega })[n_\mathrm{B}(\mathrm{\Omega })+1]}{(\omega \omega _i)^2/\eta ^2+1}},`$ (14)   
where $`\mathrm{\Gamma }_i(\omega _i)`$ is the phase relaxation rate and we have used $`\mathrm{\Gamma }_i(\omega _i)\eta `$. Thus the spectral function resembles a Lorentzian with Full Width at Half Maximum(FWHM) $`\mathrm{\Gamma }_i(\omega _i)`$ for $`\omega \omega _0`$ and it decays as $`(\omega \omega _0)^4`$ in the tails. From the experimental dephasing lifetimes we obtain $`\mathrm{\Gamma }_i(\omega _i)=1/\tau _{ph}0.7`$ cm<sup>-1</sup> at room temperature. We estimate the coupling strength using $`\delta \omega _i\mathrm{\Omega }\omega _i/4E_\mathrm{d}=4`$ cm<sup>-1</sup>, and the friction parameter can then be determined from $`\eta =2\delta \omega ^2n_\mathrm{B}(\mathrm{\Omega })(1+n_\mathrm{B}(\mathrm{\Omega }))/2\mathrm{\Gamma }_i(\omega _i)50`$ cm<sup>-1</sup>. The values of $`\delta \omega `$ and $`\eta `$ obtained in this way are similar to the measured room temperature values for Si(111).

To obtain the diffusion rate we perform the integration in Eq. (5) thus obtaining

$$w\frac{4}{\mathrm{}^2}\underset{i0}{}\frac{\chi _0^2\chi _i^2}{r_i^6}\frac{\mathrm{\Gamma }_0(\omega _0)+\mathrm{\Gamma }_0(\omega _i)}{[\omega _i\omega _0]^2+[\mathrm{\Gamma }_0(\omega _0)+\mathrm{\Gamma }_0(\omega _i)]^2/4}.$$
(15)   
In the case where $`\mathrm{\Gamma }_0(\omega _i)=\mathrm{\Gamma }_0(\omega _0)`$ this result is similar to that of Ref. .

Curves in Fig. 10 show the values of $`w`$ as obtained from Eq. (15) when the surface is subject to the same electric fields as in Fig. 7. The solid line corresponds to the electric field model of Persson and Avouris and similar to Ref.  we find $`w5\times 10^9`$ s<sup>-1</sup> for typical STM fields. The other curves in Fig. 10 and the curves in Fig. 11 show that for more realistic models of the tip electric field the value of $`w`$ is more than one order of magnitude larger, and a typical value in an STM experiment is $`w10^{11}`$ s<sup>-1</sup>. Thus, the $`n=1`$ vibrational excitation at $`r=0`$ will diffuse away very fast to the nearest neighbour sites in contrast to the result of Persson and Avouris. The reason for this is that for a realistic STM geometry the electric field at $`r=0`$ is not very different from the nearest neighbour sites and there is a large diffusion rate into these sites.

For the decay of the $`n>1`$ excitation we have to take into account the anharmonicity of the Si-H bond potential. In section II A it was shown that the bond potential of the H atom is well described by a Morse potential. The eigenstates of a Morse potential is given by

$$\mathrm{}\omega (n)=E_\mathrm{d}\left[1\frac{\alpha \mathrm{}}{\sqrt{2\mathrm{M}_\mathrm{H}E_\mathrm{d}}}(n+\frac{1}{2})\right]^2,$$
(16)   
where $`n`$ takes positive integral values from zero to the greatest value for which $`n+\frac{1}{2}<\sqrt{2\mathrm{M}_\mathrm{H}E_\mathrm{d}}/\alpha \mathrm{}`$. For the H potential $`n=0,1,\mathrm{}24`$ and $`\omega (n)=0.129,0.378,0.618,0.847,\mathrm{}`$ eV. The anharmonicity is substantial and $`\omega (2)\omega (1)=\omega (1)\omega (0)+U`$, where $`U=0.010`$ eV. The frequency of the $`n=2`$ state is outside the phonon band, and this gives rise to a localization of the state. The diffusion rate of this state can be estimated from Eq. (15) by using $`\omega _0+U`$ for the frequency at site $`0`$, and the result of such a calculation is shown by the three lower curves in figure 11. The value of $`w`$ is of the same order of magnitude as the room temperature energy relaxation rate ($`10^8`$ s<sup>-1</sup>). For $`n>2`$ the relaxation rate is $`10^8`$ s<sup>-1</sup>. Thus, it is mainly the lifetime of the $`n=1`$ excitation which is affected by incoherent diffusion. In the next section we will investigate the effect of the reduced lifetime of the $`n=1`$ excitation on STM induced desorption.

## V Calculation of the desorption rate

In this section we will calculate the desorption rate, $`R`$, of the H atom below the STM tip, due to electron inelastic scattering through dipole coupling or by resonance coupling with the Si-H $`5\sigma `$ and $`6\sigma ^{}`$ resonances. The fraction of electrons which scatters inelastically through dipole coupling is given by $`f_{in}^{dip}(\chi _0/ea_0)^20.001`$. The theoretical model we use for calculating the inelastic current, $`I_{nn+N}`$, due to resonance coupling has been described in Ref. . In those works we only considered resonance coupling and decay through energy relaxation with $`w_{\mathrm{en}}=1/\tau _{\mathrm{en}}=10^8`$ s<sup>-1</sup>, and dotted lines in Fig. 12 correspond to those results. The dashed lines show the result of including dipole coupling and the little difference between the dotted and dashed lines justify the neglect of dipole coupling in our previous studies. The solid lines in Fig. 12 show the result of including both dipole coupling and lateral diffusion of the $`n=1`$ excitation with $`w=10^{11}`$ s<sup>-1</sup>. Defining $`\mu (w)=R(w)/R(w=0)`$ as the suppression of the desorption due to lateral diffusion of the excitation, we find $`\mu 0.1`$$`0.3`$ at negative bias and $`\mu 0.02`$$`0.08`$ at positive bias. Using $`w=10^{10}`$ s<sup>-1</sup> or $`w=10^{12}`$ s<sup>-1</sup> changes $`\mu `$ less than 10 percent. This is quite different from the model of Persson and Avouris where $`\mu w/w_{\mathrm{en}}`$. The reason is that in our model we include multiple phonon excitations, i.e. we use $`N=1,2,3`$ in the calculation of the inelastic current. When the lateral diffusion rate of the $`n=1`$ level is large, the desorption proceeds via a direct excitation from $`n=0`$ to $`n=2`$. At negative biases $`<5`$ V the rate of double excitations relative to single excitations is $`I_{nn+2}/I_{nn+1}=0.07`$$`0.15\times (n+2)`$, while at positive biases $`>2`$ V it is $`0.015`$$`0.04\times (n+2)`$. Thus, the larger $`\mu `$ at negative bias relative to positive bias is due to a higher probability of a multiple excitation.

## VI Summary

We have studied the effect of electric field on incoherent lateral diffusion of vibrational excitations and its implications for STM induced desorption of H from Si(100)-H(2$`\times `$1). We calculated the electric field at the surface for realistic STM tip geometries and determined the field dependent vibrational properties of the H overlayer based on first principles calculations of vibrational Stark shifts and dipole-dipole interaction matrix elements. We found that the electric field will localize the vibrational states below the STM tip, however, the lifetime of the $`n=1`$ excitations is short ($`10`$ ps) due to incoherent diffusion. The diffusion of higher level excitations $`n>1`$ is suppressed due to anharmonic frequency shifts. The damping of the STM induced desorption of H due to the lateral escape of the $`n=1`$ excitation depends on the fraction of multiple phonon excitation events relative to one phonon events in the inelastic scattering process. At low positive biases we find a damping of the desorption rate by $`\mu 0.02`$$`0.08`$, while at negative bias $`\mu 0.1`$$`0.3`$, reflecting the higher probability of inelastic scattering events with a multiple phonon excitation at the negative biases. There are no adjustable parameters in our model and the calculated desorption rates are in quantitative agreement with measured desorption rates.

## ACKNOWLEDGMENTS

I acknowledge Jan Tue Rasmussen for making the ANSYS finite-element calculations, and thank Ben Yu-Kuang Hu, U. Quaade and F. Grey for valuable discussions and careful reading of the manuscript. This work was supported by the Danish Ministries of Industry and Research through project No. 9800466 and the use of national computer resources was supported by the Danish Research Councils.