5432
J. Phys. Chem. C 2007, 111, 5432-5440
Breaking Relaxing Bonds at a H:Si(100)-(2 × 1) Surface with Infrared Laser Pulses G. K. Paramonov,† Ivan Andrianov,‡ and Peter Saalfrank*‡ Institut fu¨r Chemie, UniVersita¨t Potsdam, Karl-Liebknecht-Str. 24-25, D-14476 Potsdam-Golm, Germany, and Institute of Physics, National Academy of Sciences of Belarus, Independence aVe. 70, 220602 Minsk, Belarus ReceiVed: NoVember 23, 2006; In Final Form: February 20, 2007
The time-dependent self-consistent field (TDSCF) approach is developed and applied to study the free relaxation and the infrared-laser-driven desorption of H atoms from a H:Si(100)-(2 × 1) surface. The method has made it possible to treat, by means of computer simulation within the Schro¨dinger wave function formalism, as many as 536 degrees of freedom of a H-Si cluster. The model used includes the H-Si stretching and the Si-Si-H bending motion of the adsorbed H atom, coupled nonlinearly to 534 degrees of freedom of the substrate, which are represented by the corresponding normal modes. It is shown that the bending vibrations relax on a picosecond time scale due to vibration-phonon coupling, in agreement with a previous perturbative study. It is also found that almost complete desorption of H atoms can be accomplished by optimally chosen, strong infrared (IR) laser pulses of 0.5 ps duration. The desorption rate is systematically studied as a function of the laser field strength for the case of plateau-type IR laser fields of 10 ps duration.
1. Introduction The problem of vibrational excitation and relaxation of adsorbates near surfaces is of significant importance in surface science.1,2 Especially the system considered here, namely hydrogen adsorbed on a Si(100)-(2 × 1) surface, has attracted great interest in spectroscopy,1 scanning tunneling microscope (STM)-induced desorption3,4 and switching,5 and UV/vis laser desorption.6 For a related system, H/D:Si(111)-(7 × 7), for which similar relaxation behavior is found,1 new experimental studies have demonstrated complete isotope and mode selectivity in the infrared-laser-induced desorption of H2 and D2 from the surface.7,8 In this case, the exciting infrared (IR) radiation was resonant with the Si-H vibrational mode, the latter being many times smaller than the Si-H bond energy. Yet, a scaling of the desorption rate proportional to only the square of the field intensity was found, rather than a high power law dependence typical for multiphoton processes. The reason for this behavior is not yet known. In our own previous work, we have investigated the problem of vibrational relaxation of H on H:Si(100)-(2 × 1) within a perturbation theory framework, developing a general model for evaluation of vibration-phonon coupling in this system.9,10 In this work, a two-mode subsystem was defined, with Si-H stretching modes and a Si-Si-H bending mode φ, which was coupled nonlinearly to a bath of up to 534 harmonic oscillators, the latter obtained from normal-mode analysis of a Si144H36 cluster; see Figure 1 for illustration. Using the Golden Rule of time-dependent perturbation theory and one- and two-“phonon” transitions within the bath, vibrational lifetimes were obtained for the Si-H stretch vibration on the nanosecond (ns) time scale and for the Si-Si-H bending vibration on the picosecond (ps) time scale. In a subsequent paper, we went beyond perturbation theory by investigating the dynamics of the complete system with the multiconfigurational time-dependent Hartree (MCTDH) method,11,12 which allowed † ‡
National Academy of Sciences of Belarus. Universita¨t Potsdam.
us to perform quantum mechanical calculations with up to 52 degrees of freedom in total.13 Again, for the Si-Si-H bending modes, a “vibrational lifetime” in the ps regime was found, with a nonexponential population decay of the excited modes. Finally, in ref 14, the possibility of IR pulse excitation of selected modes and even selected vibrational levels was investigated using Markovian and non-Markovian open-system density matrix theories and pulses either of sin2 shape or freely optimized within optimal control theory (OCT) for dissipative systems.15 It was found that even with relaxation and intermode coupling accounted for, mode-selective vibrational excitation by shaped IR pulses is possible. The more difficult task to achieve desorption, in this case of H atoms, was not attempted. In the present work, a nonperturbative system-bath approach based upon the time-dependent self-consistent field (TDSCF) form of quantum dynamics will be employed to study the problem of IR-induced desorption of H atoms from the fully H-covered Si(100)-(2 × 1) surface. The method does not rely on the assumptions of perturbation theory; yet, it allows one to include a large set of bath oscillators (in this case, all normal modes of a 180 atom cluster) in an (approximate) quantum mechanical fashion. The paper is organized as follows. In Section 2, we outline the physical model to be considered below, that is, an anharmonic, two-mode subsystem coupled to a large set of harmonic bath oscillators. In Section, 3 equations of motion within the TDSCF scheme will be derived, and numerical techniques to solve them will be outlined. In Section 4, we study first the free relaxation of H vibrations at the Si(100) surface due to vibration-phonon coupling and then present results on IRinduced desorption of H from H:Si(100)-(2 × 1). For that purpose, either single 0.5 ps pulses or 10 ps plateau-type pulses will be adopted. A final section, 5, concludes and summarizes our work. 2. The Model The overall Hamiltonian characterizing the combined total system, a H atom adsorbed on a Si dimer on a reconstructed,
10.1021/jp067796u CCC: $37.00 © 2007 American Chemical Society Published on Web 03/22/2007
Breaking Relaxing Bonds at a H:Si(100)-(2×1) Surface
J. Phys. Chem. C, Vol. 111, No. 14, 2007 5433
Figure 1. System coordinates (left) and cluster model (right) used in this work.
H-covered Si(100)-(2 × 1) surface and being excited by the laser field, can be written in the form
∂ ip ΞS(r,φ) ) H ˜ S(r,φ)ΞS(r,φ) ∂t
H ˆ )H ˆS + H ˆ SF + H ˆ SB + H ˆB
where the wave function ΞS(r,φ) is normalized with the new volume element. The modified system Hamiltonian is
(1)
where H ˆ S is the system Hamiltonian, H ˆ SF and H ˆ SB represent the system-field and the system-bath couplings, respectively, and H ˆ B is the bath Hamiltonian. 2.1. The System. The system Hamiltonian in polar coordinates r and φ reads
H ˜ S(r,φ) ) -
p 2 ∂2 p 2 ∂2 +V ˜ S(r,φ) 2mH ∂r2 2m r2 ∂φ2
(6)
(7)
H
where the modified 2D potential energy surface
p 2 ∂2 p2 ∂ p2 ∂ 2 + VS(r,φ) H ˆ S(r,φ) ) 2mH ∂r2 2mHr ∂r 2m r2 ∂φ2 H (2) where mH is the mass of the H atom, the polar coordinates r and φ represent the H-Si stretching and the Si-Si-H bending motion, respectively, and V(r,φ) is the two-dimensional (2D) potential energy surface of the adsorbed H atom, which is chosen as in ref 10
VS(r,φ) ) D[e-2R(r-r0) - 2e-R(r-r0)] + 2 2 k D + e-β(r-r0) (φ-φ0) (3) 2
Here, D ) 0.125 Eh (3.40 eV), R ) 0.83 a-1 0 , r0 ) 2.8395 a0 (1.264 Å), k ) 0.128 Eh/rad2, β ) 0.2 a-2 , and φ0 ) 1.9648 0 rad (112.6°). In passing, we note that, within a Golden Rule treatment of vibrational relaxation of the Si-H stretch mode, the inclusion of the second bending angle as an additional “system mode” would not significantly alter the relaxation dynamics, while the inclusion of the φ mode as a system mode is essential.10 The wave function ΘS(r,φ) of the time-dependent Schro¨dinger equation for the system
∂ ˆ S(r,φ)ΘS(r,φ) ip ΘS(r,φ) ) H ∂t
(4)
is normalized with the volume element dτ ) rdrdφ. It is more suitable to eliminate the second term in the right-hand side of eq 2 and to use the elementary volume dτ ) drdφ. The wellknown substitution
ΘS(r,φ) )
1 ΞS(r,φ) xr
yields the following equation of motion for the system
(5)
V ˜ S(r,φ) ) VS(r,φ) -
1 8mHr2
(8)
contains a potential-like term due to the transformation (eq 5). The wave function, the 2D potential energy surface, and the operators of the system Hamiltonian (see eqs 6-8) have been represented on a 128 × 64-point space grid for the r and φ coordinates, respectively, with r ∈ [1.3 a0, 8.0 a0] and φ ∈ [-1.0 rad, 4.92 rad]. The laser field is assumed to be linearly polarized perpendicular to the surface, that is, perpendicular to the Si dimers of the H:Si(100)-(2 × 1) surface. Its interaction with the system is described within the semiclassical electric dipole approximation by the interaction Hamiltonian
H ˆ SF(r,φ,t) ) -µS(r,φ)Ez(t)
(9)
where Ez(t) is the perpendicular component of the laser field. Further, µS(r,φ) is the z component of the dipole moment of the system, which has been chosen in the following form
µS(r,φ) ) A0 + A1(r - r0)e-A2(r-r0) + A3(φ - φ0) (10) with parameters A0 ) 0.4001 ea0, A1 ) -1.2587 e, A2 ) 0.3175 a-1 0 , A3 ) -0.4734 ea0/rad, and the equilibrium values r0 and φ0 are as those in eq 3. The parametrization (eq 10) resulted from B3LYP/6-31G** 16 cluster calculations, as outlined elsewhere.14 2.2. Zero-Order States of the System. It is suitable to analyze the vibrational relaxation and desorption in terms of the zero-order states representing (nonstationary) states of the system with explicitly defined excitation of both degrees of freedom r and φ. The 1D zero-order states of the system, 0 0 〉 and |φS,n 〉, are defined by the 1D zero-order Hamilton|ψS,m ians of the system
5434 J. Phys. Chem. C, Vol. 111, No. 14, 2007
Paramonov et al.
p2 ∂2 +V ˜ S(r,φ)|φ)φ0 2mH ∂r2
(11)
p 2 ∂2 +V ˜ S(r,φ)|r)r0 2mHr20 ∂φ2
(12)
0 H ˆ S,r (r) ) -
and 0 (φ) ) H ˆ S,φ
which are obtained from the total 2D system Hamiltonian H ˜ S(r,φ) (see eq 7) by fixing one of the degrees of freedom at the equilibrium value and thus neglecting the kinetic and 0 potential coupling between them. The zero-order states |ψS,m 〉 0 0 0 and |φS,n〉, as well as the respective energies ES,m and ES,n, are obtained from the time-independent Schro¨dinger equations 0 0 0 0 (r)|ψS,m 〉 ) ES,m |ψS,m 〉, H ˆ S,r
0 0 0 0 H ˆ S,φ (φ)|φS,n 〉 ) ES,n |φS,n 〉 (13)
where the quantum numbers m and n indicate m quanta in the stretching degree of freedom r and n quanta in the bending degree of freedom φ. The eigenproblems (eq 13) have been solved by the Fourier grid Hamiltonian method of BalintKurti.19 It was found that the 2D potential energy surface V ˜ S(r,φ) 0 〉 of the stretching supports 19 bound vibrational states |ψS,m degree of freedom r (which are of primary importance for the 0 desorption) and more than 30 bound bending states |φS,n 〉 (it was sufficient to use only the 19 lowest ones). In the following, we shall use the notation 0 0 〉|φS,n 〉 ) (m,n) |ψS,m
(14)
for the 2D direct-product zero-order states of the system, where two ordered quantum numbers, m and n, define the excitation of the r and φ degrees of freedom, respectively. The frequency of the fundamental one-photon transition (0,0) f (1,0) in the stretching mode (r) is ωr ) 2042.6 cm-1. Similarly, the frequency of the fundamental one-photon transition (0,0) f (0,1) in the bending mode (φ) is ωφ ) 645.1 cm-1. They are in good agreement with the results presented in ref 10, where the 2D Schro¨dinger equation with full 2D coupling was solved (and the small correction term in eq 8 neglected), indicating that intermode coupling within the system is weak. 2.3. The Bath. The bath Hamiltonian, based on the normalmode analysis of a cluster Si144H36 slab model (see Figure 1 and ref 10), reads NB
H ˆ B({qi}) )
∑i
(
-
p 2 ∂2
1 + Miω2i q2i 2Mi ∂q2 2 i
)
(15)
where Mi and ωi are the masses and frequencies of the normal modes of the bath, respectively. Since 180 atoms were included in the normal-mode analysis, the overall number of vibrational normal modes is NB ) 534. To obtain the normal modes, use was made of a semiempirical bond order potential by Dyson and Smith.17,18 The bath modes extend up to a calculated Debye frequency of 579 cm-1, in good accord with other sources.4 For the system-bath interaction Hamiltonian, the following form10 is assumed NB
H ˆ SB(r,φ,{qi}) )
∑i
λi(r,φ)qi +
1
NB NB
∑i ∑j
2
Λij(r,φ)qiqj (16)
where qi represents linear combinations of Cartesian displace-
ments of the individual atoms of the H-Si cluster, and λi and Λij are coupling functions that are obtained from the force field; see ref 10 for details. The system-bath Hamiltonian (eq 16) thus takes into account both one- and two-phonon transitions in the bath. In the following, the dynamics of the bath degrees of freedom will be treated by making use of the exact analytical solutions for the linearly driven harmonic oscillators (see, e.g., refs 2022 and references therein). Therefore, it is suitable to introduce the mass-weighted bath coordinate, Qi, as well as one- and twophonon coupling constants, λ˜ i and Λ ˜ ij, as follows
Qi ) xMiqi
λ˜ i(r,φ) )
λi(r,φ)
Λ ˜ ij(r,φ) )
xMi
Λij(r,φ)
xMixMj
(17)
Then, the modified Hamiltonians for the bath and the systembath interaction take the following forms
H ˜ˆ B({Qi}) )
NB
(
1
1
∑i 2Pˆ 2i + 2ω2i Q2i
)
(18)
where Pˆ i is the momentum operator adjusted to the Qi coordinate of the bath, and
H ˜ˆ SB(r,φ,{Qi}) )
NB
∑i
λ˜ i(r,φ)Qi +
1
NB NB
∑i ∑j Λ˜ ij(r,φ)QiQj
2
(19)
In passing, we note that possible IR excitation of the bath is neglected below, that is, no field-bath coupling term appears in the total Hamiltonian. While phonon excitation is possible in principle, this assumption is reasonable due to the fact that the adopted field frequencies are far off-resonant with the lowfrequency bath modes. Also, while lateral energy transfer of vibrational energy out of an isolated Si-H vibration into Si-H phonon bands is formally included, in our cluster approach, the corresponding matrix elements (in a Golden Rule-type treatment) were found to be small,10 that is, Fo¨rster transfer is not efficient in our model. 3. Equations of Motion and Techniques The time-dependent self-consistent field (TDSCF) method is often referred to as the time-dependent Hartree (TDH) method (see, e.g., ref 23 and references therein). This approach has been used in the present study because it has made it possible to treat, by means of computer simulation within the Schro¨dinger wave function formalism, as many as 536 degrees of freedom of the H-Si cluster: two degrees of freedom of the adsorbed H atom (r and φ) and all normal modes of the substrate. It follows from the previous section that the overall Hamiltonian for the H-Si cluster H ˜ )H ˜S + H ˜ SF + H ˜ SB + H ˜ B (see eq 1) can be written in the following form
H ˜ (r,φ,{Qi},t) )
[
-
p2 ∂2
p2
-
2mH ∂r2 NB
∑ i)1
∂2
2mHr2 ∂φ2
{
]
+ VS(r,φ) - µS(r,φ)Ez(t) +
1 2 1 2 Pˆ i + [ωi + Λ ˜ ii(r,φ)]Q2i + 2 2
[
λ˜ i(r,φ) +
1 2
]}
Λ ˜ ij(r,φ)Qj Qi ∑ j*i
(20)
Breaking Relaxing Bonds at a H:Si(100)-(2×1) Surface In the TDSCF approach, the overall wave function Ψ(r,φ,{Qi},t) is written as a single Hartree product
{
χi(Qi,t) ∏ i)1
∑ k*i
∂2
+
2mHr2 ∂φ2
|
〉
{
1 ˜ kk(r,φ)|ΦS〉Q2k |χk〉 + 〈χk|Pˆ 2k + 〈ΦS|Λ 2
[
〈ΦS|λ˜ k(r,φ)|ΦS〉 +
1
]
∑ 〈ΦS|Λ˜ kl(r,φ)|ΦS〉〈χl|Ql|χl〉 〈χk|Qk|χk〉
2 l*k,i
N
B B ∂ ˜ˆ | χi〉ΦS ip ΦS ) 〈 χi|H ∂t i)1 i)1
∏
p2
-
2mH ∂r2
(22)
They have, in general, the following form N
p2 ∂2
VS(r,φ) - µS(r,φ)Ez(t) ΦS +
(21)
and the TDSCF equations of motion are derived from the Dirac-Frenkel variational principle24,25
〈δΨ|ip∂t∂ - H˜ |Ψ〉 ) 0
|
κi(t) ) ΦS -
NB
Ψ(r,φ,{Qi},t) ) ΦS(r,φ,t)
J. Phys. Chem. C, Vol. 111, No. 14, 2007 5435
∏
(23)
}
(27)
respectively. They are of no importance in the following and easily eliminated by the phase transformations
[
ΦS(r,φ,t) ) exp for the wave function of the system ΦS(r,φ,t) and
i p
∫0t κS(t′)dt′]ΨS(r,φ,t)
(28)
and
∂
ip χi ) 〈ΦS ∂t
χj|H ˜ˆ |ΦS∏χj〉χi ∏ j*i j*i
i ) 1, 2, ..., NB
(24)
for the normal modes of the surface. With the overall Hamiltonian given by eq 20, the TDSCF equations of motion are specified as follows
∂ ip ΦS(r,φ,t) ) ∂t
{[
2
-
p
∂
2
-
2mH ∂r
2
µS(r,φ)Ez(t) +
i)1
(
λ˜ i(r,φ) +
1 2
∂
2
2
2mHr ∂φ
] ∑[ NB
p
2
2
+ VS(r,φ) -
i p
∫0t κi(t′)dt′]ψi(Qi,t)
i ) 1, 2, ..., NB (29)
The final TDSCF equations of motion can then be written the following form
[
p 2 ∂2 p2 ∂2 ∂ + VSCF(r,φ,t) ip ΨS(r,φ,t) ) ∂t 2mH ∂r2 2m r2 ∂φ2 H
]
µS(r,φ)Ez(t) ΨS(r,φ,t) (30)
1 Λ ˜ ii(r,φ)〈χi|Q2i |χi〉 + 2
)
[
χi(Qi,t) ) exp -
Λ ˜ ij(r,φ)〈χj|Qj|χj〉 〈χi|Qi|χi〉 ∑ j*i
]}
ΦS(r,φ,t) +
for the adsorbed H atom, where the time-dependent selfconsistent field potential VSCF(r,φ,t) is defined as NB
VSCF(r,φ,t) ) VS(r,φ) +
[
κS(t)ΦS(r,φ,t) (25)
λ˜ i(r,φ) +
for the wave function of the system (the adsorbed H atom), and
1 2
∑ i)1
{
1 Λ ˜ ii(r,φ)〈ψi|Q2i |ψi〉 + 2
]
Λ ˜ ij(r,φ)〈ψj|Qj|ψj〉 〈ψi|Qi|ψi〉 ∑ j*i
}
(31)
and
∂ ip χi(Qi,t) ) ∂t 1 2 1 2 ˜ ii(r,φ)|ΦS〉]Q2i + 〈ΦS|λ˜ i(r,φ)|ΦS〉 + Pˆ i + [ωi + 〈ΦS|Λ 2 2 1 〈ΦS|Λ ˜ ij(r,φ)|ΦS〉〈χj|Qj|χj〉 Qi χi(Qi,t) + κi(t)χi(Qi,t) 2 j*i i ) 1, 2, ..., NB (26)
1 1 ∂ ip ψi(Qi,t) ) Pˆ 2i + Ω2i (t)Q2i + Fi(t)Qi ψi(Qi,t) ∂t 2 2 i ) 1, 2, ..., NB (32)
for the normal modes of the substrate. The time-only-dependent functions κS(t) and κi(t) in eqs 25 and 26 are defined by
and
{
]}
∑
NB
κS(t) ) and
1
[
〈χi|Pˆ 2i |χi〉 ∑ i)1 2
[
]
for the substrate normal modes. Here, the time-dependent “harmonic frequencies” Ωi(t) and the time-dependent “intracluster linearly driving forces” Fi(t) are defined as
Ωi(t) ) [ω2i + 〈ΨS|Λ ˜ ii(r,φ)|ΨS〉]1/2
Fi(t) ) 〈ΨS|λ˜ i(r,φ)|ΨS〉 +
1 2
(33)
〈ΨS|Λ ˜ ij(r,φ)|ΨS〉〈ψj|Qj|ψj〉 ∑ j*i
(34)
respectively. Note that the time dependence of VSCF(r,φ,t) (eq 31), Ωi(t) (eq 33), and Fi(t) (eq 34) come from the time dependence of wave functions ψi(t) and ΨS(t).
5436 J. Phys. Chem. C, Vol. 111, No. 14, 2007
Paramonov et al.
The wave functions and the operators of the equation of motion for the adsorbed H atom (eq 30) have been represented on the 128 × 64-point space grid specified in Section 2.1. The time propagation has been accomplished by a split operator method,26,27 with the fast Fourier transform procedure used to obtain the spatial derivatives at each time grid point. The time step of propagation, ∆t, ranged from 0.2 to 0.6 p/Eh. The time-dependent populations of the bound states (m,n) of the adsorbed H atom are defined by projection of the timedependent wave function ΨS(t) on the 2D direct-product zeroorder states (eq 14) 0 0 P(m,n)(t) ) |〈φS,n (φ)ψS,m (r)|ΨS(t)(r,φ)〉r,φ|2
(35)
At large values of the r coordinate, the absorbing boundary conditions have been provided by an imaginary smooth optical potential Vopt(r)
{[
Vopt(r) ) -iV0 exp
(rmax - ropt) 3 12 (r - r )2
]}
2
(36)
Figure 2. Vibrational relaxation of φ-excited states (0,n), for n ) 1, 2, and 3, due to vibration-phonon coupling.
if r g ropt, and Vopt(r) ) 0 otherwise (for further details, see refs 28-30 and references therein). The results presented below have been calculated at V0 ) 0.083 Eh and ropt ) 7.2 a0. The time-dependent desorption probability PD(t) has been defined by the time- and φ-integrated outgoing flux as follows
field being switched off shows that P(0,0)(t) > 0.9996 up to t ) 10 ps if ΨS(r,φ,t ) 0) ) ΨS,G(r,φ) (see eqs 35 and 14). In contrast, in the present subsection, the wave function with n quanta in the φ zero-order state
PD(t) )
p mH
∫0t dt′ ∫φφ
max
min
opt
[
]|
∂ΨS(r,φ,t′) ∂r
dφ Im Ψ*S(r,φ,t′)
Ψn,0(r,φ,{Qi}) ) r)ropt
(37) The equations of motion of the bath modes (eq 32) represent the linearly driven harmonic oscillators. The analytical solutions of these equations are well-known (see, for example, refs 20 and 21 and references therein). Therefore, all time-dependent expectation values, 〈Ψi(t)|Qi|Ψi(t)〉 and 〈Ψi(t)|Q2i |Ψi(t)〉, are reduced to the tabulated integrals.31 This remarkably facilitates the efficiency of the numerical routine. A similar approach has been already used in previous work.22 4. Results and Discussion 4.1. Free Decay of Vibrationally Excited States. As a first application of the TDSCF method, the free decay of an excited H vibration due to vibration-phonon coupling was studied. For this purpose, we first assumed that all degrees of freedom were in their ground states. Then, the “true” eigenfunction of the ground state of the entire cluster, ΨG(r,φ,{Qi}), was calculated by the relaxation method,32 that is, propagating the TDSCF equations of motion in imaginary time, with the laser field being switched off, Ez(t) ) 0. Thus obtained, the stationary, ground states |ΨS,G〉 and |ψi,G〉, i ) 1, 2, ..., NB correspond to the overall ground state |ΨG〉 of the entire H-Si cluster represented by the wave function NB
ΨG(r,φ,{Qi}) ) ΨS,G(r,φ)
ψi,G(Qi) ∏ i)1
NB
(38)
The ground state (eq 38) has been used in Sections 4.2 and 4.3 as the initial one in all simulations of the laser-driven desorption of the H atom from the Si(100) surface. Note that the ground state of the system |ΨS,G〉 is very well represented by the 2D direct-product zero-order state (m ) 0,n ) 0); propagation of the TDSE equations of motions in the real time with the laser
0 ψ0S,0(r)φS,n (φ)
ψi,G(Qi) ∏ i)1
(39)
0 (φ) is was used as an initial state. Again, the product ψ0S,0(r)φS,n not an exact eigenfunction of the system Hamiltonian (eq 7) but comes close to it if n is not too large. Using exact eigenfunctions of eq 7 as the zero-order states has only a very small influence on the results, as demonstrated below. Initial excitation of the r mode was not studied because, in this case, the relaxation is very slow, nanoseconds according to refs 1 and 10. Using φ-excited states, (0,n) n ) 1,2,3, we note from Figure 2 that these “decay” on a ps time scale due to coupling to substrate phonons. Shown in the figure are the populations of the zero-order states as a function of time; see eq 35. States which are not shown have only very little population. It is also observed that (i) the relaxation process is nonexponential, (ii) the faster the relaxation process, the larger n is, and, finally, (iii) recurrences occur at longer times. Due to the recurrences and the nonexponential decay, it is not straightforward to define a lifetime. We rather define a half-lifetime, T1/2(n) for state (0,n), the time after which the population in this state has decayed to 1/2. This half-lifetime can be compared with the one obtained from the perturbative model of ref 10, where the same model has been used, however, within Fermi’s Golden Rule of perturbation theory. In this case, one has a strictly exponential decay, leading to lifetimes τ(n), which are related to half-lifetimes by T1/2(n) ) τ(n) ln 2. From Table 1, where the half-lifetimes are shown, we note that the agreement between the perturbative and nonperturbative treatments is almost perfect as far as half-lifetimes are concerned. This demonstrates that, even for the faster decaying φ-excited modes, the system-bath coupling is obviously weak enough to make a perturbative approach appropriate. Nevertheless, the good agreement between perturbative and nonperturbative theories remains astonishing, in particular also since, in ref 10, slightly different initial states (namely exact eigenstates of the twodimensional-system Hamiltonian) and even a slightly different
Breaking Relaxing Bonds at a H:Si(100)-(2×1) Surface
J. Phys. Chem. C, Vol. 111, No. 14, 2007 5437
TABLE 1: Half-lifetimes of Low-Lying Vibrational Energy Levels (0,n)a
state
T1/2(n) (TDSCF) NB ) 534 (ps)
T1/2(n) (PT) NB ) 534 (ps)
T1/2(n) (MCTDH) NB ) 50 (ps)
(0,1) (0,2) (0,3)
0.92 (0.95/0.91) 0.475 (0.49/0.47) 0.335 (0.34/0.32)
0.94 0.48 0.33
3.54 1.13 0.87
a The first entry is the TDSCF result with 534 oscillators, the second one is a Golden-Rule treatment with 534 oscillators,10 and the last entry is from a MCTDH calculation with 50 bath oscillators.13 The first values in brackets of the second column are for TDSCF with exact eigenfunctions of the system hamiltonian (eq 7) used as zero-order states rather than the product of one-dimensional states (eq 14). The second numbers in brackets are for direct products of one-dimensional states computed without the additional potential term in eq 8.
system Hamiltonian (namely without the additional potential term in eq 8) have been used. Obviously, the adopted approximations work well. This is also further demonstrated by the second column of Table 1, where the first numbers in brackets denote half-lifetimes when calculated with the exact eigenfunctions of the two-dimensional-system Hamiltonian (eq 7). We have also used direct products of one-dimensional states computed without the additional potential term in eq 8, whose half-lifetimes are shown as the second numbers in brackets. For all cases, the differences are small; the agreement between TDSCF and PT remains excellent. The faster decay of higher excited states according to Table 1 follows, in good approximation, the law
T1/2(n) )
T1/2(n ) 1) n
(40)
Equation 40 holds exactly in time-dependent perturbation theory if a harmonically idealized system vibration-couples linearly to the bath of oscillators.10 In ref 10, the approximate relationship (eq 40) has been demonstrated for the bending mode up to n ≈ 10 (see Table III in that reference). In Table 1, we also show the results of the recent MCTDH study on the same system.13 This method is also nonperturbative and gives, when properly converged, an exact solution to the time-dependent Schro¨dinger equation but is more expensive than the TDSCF method. In ref 13, only 50 bath oscillators could be computationally afforded. The resulting half-lifetimes, are by a factor of about 3, longer in this case, showing that a much larger number of bath oscillators is needed to converge the lifetimes. (It must also be said, however, that the decay in MCTDH is very nonexponential, making the half-lifetime a less suitable measure for the relaxation.13) Despite the good agreement between perturbation theory and TDSCF for the half-lifetimes, differences remain at early times (nonexponential behavior for TDSCF) and at later times (nonmonotonic behavior of the decaying state in TDSCF due to recurrences). The recurrences are due to the finite bath size and appear earlier for higher excited initial states. appear the higher excited the initial state was. More interesting is the fact that, in contrast to perturbation theory, where the system energy eventually decays to zero and all populations end up in the (0,0) state, the nonperturbative approach leads to a “thermal ensemble” with finite temperature. In fact, after 10 ps or so, at least the φ-excited states are Boltzmann populated. This is demonstrated in Figure 3, which shows the populations of the lowest seven φ-excited states, (0,0)-(0,6), after 10 ps of propagation time, when the same three initial states as those in
Figure 3. Populations P(0,n)(10 ps) for various φ-excited states, for different initial states.
Figure 2 were used. In Figure 3, we plot P(0,n) (10 ps) as a function of n. Since the φ-excited states are very harmonic (E(0,1) - E(0,0) ) 645 cm-1, E(0,6) - E(0,5) ) 627 cm-1), the perfect exponential behavior in Figure 3 reflects the Boltzmann law P(0,n) ∝ exp(-E(0,n)/kBT). A temperature, T, can be extracted for each initial state, T ≈ 800 K for (0,1), T ≈ 1500 K for (0,2), and T ≈ 2200 K for (0,3), showing that the φ mode is very hot after 10 ps. These temperatures do not apply for the entire cluster (and also not for the r mode), though, that is, thermal equilibrium is not reached because the φ mode is only strongly coupled to a subset of cluster modes. 4.2. Desorption Controlled by Strong IR Laser Pulses. We now study the possibility of desorbing single H atoms from a Si(100)-(2 × 1) surface using strong infrared laser pulses. We adopt single IR pulses which are linearly polarized perpendicular to the Si(100) surface, having a sin2 shape which starts at t)0
Ez(t) ) E0 sin2(πt/tp)sin(ωlt)
(41)
where E0 is the amplitude of the pulse, tp is the pulse duration, and ωl is the laser carrier frequency. In the following, the laser pulse duration is fixed to tp ) 0.5 ps, while the pulse amplitude E0 and its carrier frequency ωl are optimized such as to provide the maximum desorption probability PD at the end of the pulse. In the first stage of the optimization procedure, the amplitude of the IR laser pulse is fixed to E0 ) 0.0462 Eh/(ea0), or 237.4 MV/cm. This corresponds to a peak intensity Imax ) (1/2)0cE20 of 7.48 × 1013 W/cm2. With this high intensity, multiphoton transitions are possible. We calculate the desorption probability at the end of the pulse, PD(t ) tp), for various laser frequencies ωl corresponding to multiphoton transitions in the stretching mode (r) as follows
(0,0) f m photons f (m,0)
m ) 1, 2, ..., 18
(42)
0 that is, pωl(m) ) (ES,m - E0S,0)/m. Due to the anharmonicity of the potential along r, ωl(m) decreases gradually from ωl(1) ) ωr ) 2043 cm-1 with increasing m. For example, for m ) 10, ωl(10) ) 1672 cm-1. The results presented in Figure 4 show, first of all, that under the chosen conditions desorption is possible, with a large probability, and that the optimal laser frequencies correspond to the 5 photon, 7 photon, and 10 photon resonance transitions in the stretching mode. The desorption probability is drastically reduced at m < 4 and m > 12. The r mode has a vibrational lifetime in the ns regime, and even states closer to the desorption
5438 J. Phys. Chem. C, Vol. 111, No. 14, 2007
Paramonov et al.
Figure 4. Desorption probability at the end of the 0.5 ps laser pulse as a function of the discrete laser frequency corresponding to the exact m photon resonance transition (0,0) f (m,0) in the stretching mode of H-Si. The laser pulse amplitude is E0 ) 237.4 MV/cm.
Figure 5. Desorption probability at the end of the 0.5 ps laser pulse as a function of the amplitude of the pulse, E0. The laser carrier frequency, ωl ) 1672.34 cm-1, corresponds to the exact 10 photon resonance transition (0,0) f (10,0) in the stretching mode of H-Si.
energy of 0.125 Eh (3.4 eV) survive long enough for efficient ladder climbing with high-intensity laser pulses. The optimal frequency ωl(10) ) 1672 cm-1 is a good compromise to account for the reduced energy spacing of the levels close to the desorption continuum and the larger energies that are needed to climb the lowest steps of the vibrational ladder. As a fine detail, we note that the desorption probabilities are comparatively low for m ) 6 and 8. This is mainly due to nonlinear multiphoton transitions33 of the type
(k,0) f m′ photons f (l,0)
(43)
where the number of absorbed photons m′ is not equal to the number of excited vibrational quanta (l - k), in contrast with the linear multiphoton transitions defined by eq 42. Nonlinear multiphoton transitions always exist in anharmonic vibrational ladders being pumped on a linear multiphoton resonance and can play both a positive and a negative role.33 For example, in the case of m ) 8, or the 8 photon pumping of the linear transition (0,0) f 8 photons f (8,0), there exists a 5 photon quasiresonant nonlinear transition (8,0) f 5 photons f (16,0), being detuned from the exact resonance by 2 cm-1 only. Therefore, the excitation pathway is initially as follows: (0,0) f 8 photons f (8,0) f 5 photons f (16,0) f ... f continuum. Approximately in the middle of the laser pulse, the state (16,0) acquires a substantial population (up to 20%), and the backward population flow down the vibrational manifold begins to compete against desorption, which explains the lower desorption probability. For example, the final population of the ground state (0,0) at the end of the laser pulse is 0.105. A similar situation takes place for m ) 6. In contrast, for m ) 7, for example, nonlinear transitions are found to open efficient channels into the continuum, while population flow down the vibrational ladder plays only a minor role. In the second stage of the optimization, the laser carrier frequency ωl is chosen to be in the exact 10 photon resonance with the (0,0) f (10,0) transition (i.e., ωl ) 1672 cm-1), and the desorption probability at the end of the laser pulse is calculated for various amplitudes E0 at the same pulse duration tp ) 0.5 ps. The results presented in Figure 5 demonstrate that, with the optimally chosen strong IR laser pulses, the desorption of the H atom from the Si(100) surface can be accomplished with almost 100% probability on the subpicosecond time scale, with intense laser pulses. 4.3. Desorption with Plateau-Type IR Laser Fields. In this section, the desorption rate of the H atom from the Si(100) surface is systematically studied as a function of the laser field strength for the case of long, plateau-type IR laser fields. The plateau-type IR laser field is, again, assumed to be linearly
Figure 6. Time-resolved desorption probability as a function of the amplitude of the pulse, E0 (in units of Eh/ea0). The laser carrier frequency, ωl ) 2042.6 cm-1, corresponds to the exact 1 photon resonance transition (0,0) f (1,0) in the stretching mode of H-Si.
polarized perpendicular to the Si(100) surface and has the following form
Ez(t) )
{
E0 sin2(πt/2t1)sin (ωlt) (0 e t e t1) E0 sin(ωlt) (t g t1)
(44)
The field amplitude increases in the interval 0 e t e t1 and remains constant (E0) at t g t1. The results presented below are calculated at the raising time of t1 ) 0.3 ps. The overall time of excitation is t ) 10 ps, and the laser carrier frequency is fixed to ωl ) 2042.6 cm-1, which corresponds to the exact 1 photon resonance for the (0,0) f (1,0) transition in the stretching mode r. With the longer pulses, appreciable desorption probabilities can be achieved with smaller field intensities. The situation is closer to the experimental one of ref 7, where a long “train” of picosecond IR pulses has been used to desorb molecular hydrogen from H:Si(111). In Figure 6, we show the probability for IR desorption of H atoms from H:Si(100) for different field amplitudes E0. From the figure, several observations can be made. (i) It takes a longer time for the less intense fields to cause desorption. (ii) The desorption yield depends strongly on the field amplitude E0; more quantitatively, one finds PD(10 ps) ) 1.57 × 10-7 for E0 ) 0.009232 Eh/a0, PD(10 ps) ) 4.26 × 10-2 for E0 ) 0.0018464 Eh/a0, and PD(10 ps) ) 0.749 for E0 ) 0.027696 Eh/a0. (iii) The PD(t) curve has a nontrivial shape, having an oscillatory, step-like behavior (see, for example, E0 ) 0.013848 Eh/a0) and saturation phenomena (for the larger E0). This indicates a complicated nonequilibrium situation, an interpreta-
Breaking Relaxing Bonds at a H:Si(100)-(2×1) Surface
J. Phys. Chem. C, Vol. 111, No. 14, 2007 5439
Figure 7. Zero-order state populations P(m,0) (circles, solid line) and P(m,1) (squares, dashed line) after t ) 10 ps for E0 ) 0.018464 Eh/a0.
averaged over certain time intervals, or desorption probabilities at times other than those shown. When plotting these quantities versus intensity, we find, in all cases, similar behavior as that demonstrated in Figure 8, that is, high average exponents, k, around 7-10. Since direct dipole coupling is at work, k can be interpreted as the effective number of photons required for photodesorption. In a truncated harmonic oscillator model, this would correspond to the number of bound states. For the anharmonic, dissipative system considered here, k is smaller than the number of bound states along r (eq 19), but it is still large. In particular, k is significantly larger than the value of k ≈ 2 found in ref 7 for the desorption of H2 from H:Si(111). Hence, the surprising scaling with IR intensity in ref 7 remains unexplained. It is also clear, however, that H desorption from Si(100) is different from H2 desorption from Si(111) and that different laser pulses were applied in both cases. It will be interesting to model H2 desorption from H-covered silicon surfaces under conditions similar to those of ref 7. 5. Summary and Conclusions
Figure 8. Scaling of PD(2 ps) (solid line, circles) and PD(10 ps) (solid line, squares) with the field amplitude squared, E0. The two dashed lines are power law fits ∝ I6.60 (upper curve) and ∝ I7.68 (lower curve), respectively.
tion which is supported by the population of bound states along the vibrational ladder. In Figure 7, the zero-order state populations P(m,0) and P(m,1) are shown after t ) 10 ps for the case with E0 ) 0.018464 Eh/a0. The figure demonstrates that the higher levels are generally less populated; however, it is clearly impossible to assign, in a meaningful way, a temperature to the adsorbatesurface mode r, as is frequently done for the interpretation of STM-induced surface reactions where similar ladder climbing occurs.34 The lack of a temperature is due to the weak coupling of r to any other modes, in contrast to φ, as argued in Section 4.1. From the figure, we also note that states (m,1) play a much smaller role than (m,0) (and in general, φ-excited levels carry a very small population), which emphasizes the highly unidirectional reaction pathway (along r) with little r-φ coupling. For the plateau pulses, we have also monitored the desorption yield as a function of laser intensity, I. In Figure 8, we plot PD after two different times, t ) 2 and t ) 10 ps, as a function of E20 (which is proportional to the laser intensity). The shorter time is not plagued by saturation effects, which occur for t ) 10 ps at larger field strengths (see Figure 6). Since, in Figure 8, a double-logarithmic scale is used, a linear curve would correspond to a power law
PD ∝ I k
(45)
According to Figure 8, neither PD(2 ps) nor PD(10 ps) follow a strict power law; rather, the curves are non-monotonic. Still, a power-law-like trend is observed as indicated by the dashed lines, which are fits to the computed curves assuming the functional form (eq 45). The best fits are PD(2 ps) ∝ I7.68 and PD(10 ps) ∝ I6.60. We have also investigated other quantities such as desorption rates RD(t) ) dPD(t)/dt at various times, rates
In summary, the TDSCF method has been applied to study the desorption of H atoms from a Si(100) surface induced by IR laser pulses. The general methodology is of the systembath type but is nonperturbative, that is, an anharmonic subsystem is coupled nonlinearly to a set of harmonic oscillators with one- and two-phonon transitions taken into account. There are many situations in physical chemistry where similar situations apply, not only in surface science but also in other condensed phases and for large molecules. For the present example, perturbation theory works well, but in other instances, this may not be the case, and a full solution of the Schro¨dinger equation, even if only approximate, is required. At any rate, the solution of the Schro¨dinger equation, if possible, is superior to approximate treatments such as open-system density matrix theory, and details of the dynamics at early and late times can be very different, even for H/Si(100). As far as the IR-induced desorption process is concerned, we have shown that it should be possible to selectively excite the Si-H stretch mode in the presence of a multitude of other modes, which eventually leads to desorption. The desorption yield can be very high with intense laser pulses. Even under moderate, lower intensity conditions, desorption of H atoms occurs. Then, the desorption yield, the rate, or both increase roughly according to a high-power law with intensity I, supporting a multiphoton mechanism. Acknowledgment. We acknowledge generous support from the Deutsche Forschungsgemeinschaft (through Project Sa 547/7 and, in part, through Sfb 658) and from the Fonds der Chemischen Industrie. G.K.P. acknowledges the hospitality of J. Manz (Berlin) and access to computers. References and Notes (1) Guyot-Sionnest, P.; Lin, P.; Miller, E. J. Chem. Phys. 1995, 102, 4269. (2) Lass, K.; Han, X.; Hasselbrink, E. J. Chem. Phys. 2005, 123, 051102. (3) Foley, E.; Kam, A.; Lyding, J.; Avouris, P. Phys. ReV. Lett. 1998, 80, 1336. (4) Persson, B.; Avouris, P. Surf. Sci. 1997, 390, 45. (5) Stokbro, K.; Quaade, U.; Lin, R.; Thirstrup, C.; Grey, F. Faraday Discuss. 2000, 117, 231. (6) Vondrak, T.; Zhu, X.-Y. Phys. ReV. Lett. 1999, 82, 1967. (7) Liu, Z.; Feldman, L.; Tolk, N.; Zhang, Z.; Cohen, P. Science 2006, 312, 1024. (8) Tully, J. Science 2006, 312, 1004.
5440 J. Phys. Chem. C, Vol. 111, No. 14, 2007 (9) Andrianov, I.; Saalfrank, P. Chem. Phys. Lett. 2001, 350, 191. (10) Andrianov, I.; Saalfrank, P. J. Chem. Phys. 2006, 124, 034710. (11) Meyer, H.-D.; Manthe, U.; Cederbaum, L. Chem. Phys. Lett. 1990, 165, 73. (12) Beck, M. H.; Ja¨ckle, A.; Worth, G. A.; Meyer, H.-D. Phys. Rep. 2000, 324, 1. (13) Andrianov, I.; Saalfrank, P. Chem. Phys. Lett. 2006, 433, 91. (14) Paramonov, G. K.; Beyvers, S.; Andrianov, I.; Saalfrank, P. Phys. ReV. B 2007, 75, 045405. (15) Ohtsuki, Y.; Zhu, W.; Rabitz, H. J. Chem. Phys. 1999, 110, 9825. (16) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (17) Dyson, A.; Smith, P. Mol. Phys. 1999, 96, 1491. (18) Sbraccia, C.; Silvestrelli, P.; Ancilotto, F. Surf. Sci. 2002, 516, 147. (19) Marston, C. C.; Balint-Kurti, G. G. J. Chem. Phys. 1989, 91, 3571. (20) Pechukas, P.; Light, J. C. J. Chem. Phys. 1966, 44, 3897. (21) Stone, J.; Thiele, E.; Goodman, M. F. J. Chem. Phys. 1975, 63, 2936.
Paramonov et al. (22) Paramonov, G. K.; Naundorf, H.; Ku¨hn, O. Eur. Phys. J. D 2001, 14, 205. (23) Gerber, R. B.; Buch, V.; Ratner, M. A. J. Chem. Phys. 1982, 77, 3022. (24) Dirac, P. A. M. Proc. Cambridge Philos. Soc. 1930, 26, 376. (25) Frenkel, J. WaVe Mechanics; Clarendon Press: Oxford, U.K., 1934. (26) Feit, M. D.; Fleck, J. A. J. Comput. Phys. 1982, 47, 412. (27) Feit, M. D.; Fleck, J. A. J. Chem. Phys. 1982, 78, 301. (28) Vibo´k, A Ä .; Balint-Kurti, G. G. J. Chem. Phys. 1992, 96, 7615. (29) Vibo´k, A Ä .; Balint-Kurti, G. G. J. Phys. Chem. 1992, 96, 8712. (30) Kaluzˇa, M.; Muckerman, J. T.; Gross, P.; Rabitz, H. J. Chem. Phys. 1994, 100, 4211. (31) Gradshteyn, I. S.; Ryzhik, I. M. Tables of Integrals, Series, and Products; Academic Press: New York, 1980. (32) Kosloff, R.; Tal-Ezer, H. Chem. Phys. Lett. 1986, 127, 223. (33) Paramonov, G. K. Chem. Phys. 1993, 177, 169. (34) Gao, S. Phys. ReV. B 1997, 55, 1876 and references therein.
