- Email: [email protected]
Si(111) 1×1 surfaces
Chemical Physics ELSEVIER
Chemical Physics 200 (1995) 357-368
Vibrational energy relaxation dynamics of Si-H stretching modes on stepped H/Si ( 111 ) 1 × 1 surfaces Y i n g - C h i e h S u n l, H u a d o n g G a i 2, G r e g o r y A . V o t h Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323, USA Received 12 May 1995
Abstract The vibrational energy relaxation rates of excited Si-H stretching modes on the monohydride step of miscut H/Si( I 11 ) 1 × 1 surfaces are calculated using Bloch-Redfield theory combined with classical molecular dynamics (MD) simulation. The structure and vibrational frequencies of the surface are first investigated using the Car-Parrinello ab initio MD method. The calculated Si-Si-H bending frequencies and relaxed structures are then used to refine the empirical potential for the classical MD simulations. The lifetime of the excited Si-H stretching mode at the step is found to be shorter than the modes on the terrace. Both the magnitude and the trend of the calculated results agree well with the experimental measurement on the 9 ° monohydride stepped surface. The vibrational relaxation rate of the Si-H stretching modes on the 15° monohydride stepped surface are also calculated and predicted to have a slightly shorter lifetime than for the 9 ° surface.
1. Introduction The area around steps on solid surfaces is one of the most interesting zones for many physical and chemical processes at such interfaces. For example, the steps on etched surfaces usually act as the nucleation sites in the crystal growth process. The appearance of the steps on the surface may break chemical bonds of adsorbates and produce new chemical species on the surface [ 1 ]. These processes can also be substantially influenced by the adsorbate dynamics. A specific case of the latter is the hydrogen dynamics at the step of
t Present address: School of Pharmacy, 926-S, Department of Pharmaceutical Chemistry, Universityof California, San Francisco, San Francisco, CA 94143, USA. 2 Present address: EnvironmentalMolecular Sciences Laboratory, Battelle/Pacific Northwest Laboratories, KI-90, Richland, WA 99352, USA.
etched H / S i ( 111 ) surfaces which may play a crucial role in the production of semiconductor devices. Experimentally, significant progress has been made to reveal the hydrogen dynamics on Si( 111 ) surfaces [ 2 - 1 6 ] . Many aspects, for instance, the vibrational frequencies, growth geometries etc., of the stepped H / S i ( 111 ) surfaces were uncovered by various experimental methods. [3,9-12,14] The studies of the H-terminated silicon surfaces with infrared spectroscopy have been summarized in a review article by Chabal [ 13]. Because of the development of sum frequency generation(SFG) spectroscopy (see, e.g., Refs. [ 17,18 ] ), the vibrational modes can be resolved in both the frequency and time domains. For example, the Si-H stretching frequency on the flat H-terminated surface was measured using SFG spectroscopy and its lifetime was found to be about 1.0 ns at 300 K [4]. Due to the existence of the steps on the etched Si( 111 ) surface, the involvement of additional energy
0301-0104/95/$09.50 Q 1995 Elsevier Science B.V. All fights reserved SSD! 0301-0104(95)00270-7
358
Y.-C. Sun et al./Chemical Physics 200 (1995) 357-368
transfer pathways complicates the vibrational energy relaxation processes of the excited Si-H stretching modes on the stepped H / S i ( I 11 ) surfaces. In recent years, experiments have been carried out [15,16] to probe the vibrations and energy transfer mechanisms on stepped H/Si( 111 ) 1 x I surfaces using pump/probe SFG spectroscopy. The frequencies of the Si-H stretching modes at the step and on the terrace, as well as their lifetimes, were determined. On the 9 ° stepped monohydride surface, the Si-H stretching mode at the step was found to be 2089 c m - J . In addition, the measurements with two-color pump/probe spectroscopy have provided extensive information on the energy flow pathways on these surfaces. The spectral effects due to the renormalization of vibrational modes after excitation have also been calculated in a recent paper by Kuhnke et al. [ 19]. The latter effects account for the unusual spectral transients observed in the pump-probe measurements. Although the picture of energy flow on the dihydride H-terminated H/Si( 111 ) 1 × 1 surface was mapped out using on a kinetic model [ 16], there remains an ambiguity for the relaxation rates for some energy transfer channels of the monohydride stepped H / S i ( 111 ) 1 x I surface due to the various possible effects from coupling to surface phonons. The energy relaxation mechanism of Si-H stretching modes on the flat H/Si( 111 ) 1 x 1 surface has been extensively examined [2-16]. The long lifetime is presumably related to the large frequency mismatch between the Si-H stretching modes and the rest of the phonon modes on the surface. The temperature dependence of the lifetime suggests that the dominant relaxation mechanism is a four-phonon process [5]. A previous theoretical study of this system [20] was based on the quantum Redfield theory of relaxation [21,22] and obtained good agreement between the measured and calculated lifetimes for the Si-H stretching mode on the flat H / S i ( 111 ) 1 x 1 surface. The relaxation rate of the single excited Si-H surface bond into the substrate was determined by calculating the population evolution of vibrational states which, in turn, are influenced by the coupling to the substrate phonons. Specifically, the fluctuating force from the substrate calculated from classical molecular dynamics (MD) and projected along the Si-H bond coordinate contains the information on the latter couplings. On the monohydride stepped H/Si( 111 ) 1 × 1 sur-
face, experimental evidence suggests that the hydrogen atoms at the step terminate the surface in a similar way to the hydrogen on the terrace [ 12]. The silicon atoms at the step are therefore tetrahedrallycoordinated. The direct application of our simulation approach [20] for the relaxation of Si-H stretching modes of flat surfaces is therefore reasonable, enabling us to calculate the Si-H vibrational energy lifetimes both at the step and on the terrace. On the fiat surface, a large portion of the Si-H stretch vibrational energy slowly relaxes into a state consisting of three Si-Si-H bending phonons plus a surface phonon. This behavior occurs since the Si-H stretching motion is kinetically coupled with the adjacent Si-Si-H bending mode [20]. Because of this coupling, the bending frequency has a critical effect in the stretching mode relaxation. While the Si-Si-H bending frequency on the flat H/Si( 111 ) 1 × 1 surface has been measured to be 640 cm - l [23,24], to our knowledge the bending frequency at the monohydride step has not been measured. It was therefore necessary to first estimate this bending frequency and refine our empirical potential [20] before the lifetime of the Si-H stretching mode at the step could be calculated using the BlochRedfield/MD methodology. Besides the bending frequency, the geometrical structure of the stepped surface may also have a significant effect on the relaxation rate of the Si-H stretching modes. While the structure of the dihydride stepped H / S i ( 111 ) 1 × 1 surface has been studied both experimentally and theoretically [ 12,25,26], the structure of the stepped monohydride H/Si( 111 ) 1 × 1 surface has not been examined by direct calculation. Fortunately, both the Si-SiH bending frequency at the monohydride step and the step structure can be examined using ab initio methods. Accordingly, the first-principles Car-Parrinello molecular dynamics (CP-MD) method [27], which has been widely used to obtain dynamical and structural information on Si(l 11) surfaces, [28-32] was first implemented to study the monohydride stepped H/Si( 111 ) 1 x I surface. (See also Ref. [26] for an ab initio study of the dihydride stepped H/Si( 111 ) 1 x 1 surface.) In our CP-MD calculations (cf. Section 3), the bending frequency and structure were obtained using the ab initio simulation, and then this information was used to refine the empirical potential used in the Bloch-Redfield/classical MD simulation method
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
[20]. The sections of this paper are organized as follows: Section 2 briefly describes the Bloch-Redfield/MD methodology as it was applied in the present calculations. Next, Section 3 contains the description of the Car-Parrinello molecular dynamics simulations to estimate the geometry of the stepped surface as well as the vibrational frequencies shifts at the step. The refinement of the empirical potential from the results of the Car-Parrinello method simulation and the resulting Bloch-Redfield/MD simulations are then described in Section 4. A discussion of the results is given in Section 5 and the factors which affect the relaxation rate of the Si-H stretching modes are analyzed. Concluding remarks are given in Section 6.
2. Theory The description of Bloch-Redfield theory [21,3335 ] as it applies to the relaxation of an excited stretching state on a surface [20] is given briefly as follows: The total Hamiltonian is written as H = Hs(q) +
Ha(Q) + Hint(q, Q ) ,
(2.1)
where q and H s ( q ) are the coordinate and Hamiltonian of the quantum Si-H oscillator, respectively. The remaining "bath" coordinates and Hamiltonian are represented by Q and HB(Q), respectively, while the interaction between the system and the bath is described by the term Hint(q, Q). Several assumptions are made in the derivation of the Redfield equations [ 34]. First, it is assumed that the interaction Hamiltonian can be written a s H i n t ( q , Q ) = qrF(Q). This assumption is reasonable since at 300 K the dominant energy transfer should proceed one quantum at a time. 3 The second assumption is that the bath remains close 3 This assumption amounts to the linear response assumption. Since the Si-H oscillators are somewhat anharmonic, coupling terms like q2G(Q) could in principle contribute to the relaxation rate since even for one-phonon processes the matrix element (0lq211) (nlG(Q)In') could be non-zero. Though it is difficult to estimate the size of such couplings, it is unlikely that they are large in the present case because the Si-H oscillator is nearly harmonic in the lower regions of the well (i.e., (01q211) is very small). Furthermore, if such terms were significant the present one-phonon linear response-based approach would be incapable of describing the relaxation of the v = 2 Si-H mode which is not the case [20].
359
to thermal equilibrium during the relaxation process. This is valid because the amount of the energy transferred from the excited stretching states to the bath is small relative to the total energy of the bath and equilibrium for the bath can be re-established quickly after a small perturbation. The third assumption is that the relaxation time of the oscillator is much longer than the force-force autocorrelation time of the bath. This condition can be justified in that the time scale of the correlation function in the bath is within the picosecond range which is much shorter than the nanosecond lifetime of the excited Si-H stretching mode on the H / S i ( 111 ) 1 x 1 surface. The dynamical information of the Si-H system coupled to the bath is contained in the time evolution of the reduced density matrix elements for the system. The evolution of the reduced density matrix elements is governed by a series of linearly coupled differential equations which can be derived under the above assumptions according to the Bloch-Redfield theory. These equations are given by
dpit (t) dt
Z Riljk ei(°~'J+°~kDtpjk(t) , j,k pit(t) is the reduced density
(2.2)
where matrix element of the Si-H oscillator stretching states Ii) and Il) in the interaction representation, and o)i and o)l are the eigenfrequencies of the li) and Il) states respectively. The Riljk t e r m s are given by
r
r
(2.3) where the coupling coefficients between the density matrix elements are given by
Fk+tij x
1 h 2 ( 1 + e#O~,j) +oo
/.
dte -'~°'jt
(klqll) (ilqlj)
(6F(O)6F(t))cl
(2.4)
--00
and
F~O = x
1 h2(l + e-#~*,)
/.
(klqll)(ilqlJ)
dte -'~°'a (~F(O)~F(t))c~.
--(X)
(2.5)
360
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
In the above equations, (klqll) is the transition matrix element between Si-H stretching states [k) and [l), t~kt is the Kronecker delta, and the dynamical function 8F(t) is the classical fluctuating force on the Si-H bond arising from the coupling to the other modes of the system. The use of a classical force autocorrelation function in Eq. (2.4) and Eq. (2.5) is an approximation to a fully quantum mechanical Redfield theory. A standard symmetrization procedure [ 36] is employed in the derivation of quantum Redfield equations and then the resulting correlation function is replaced with the classical correlation function. It has been shown in a recent study by Bader and Berne [37] that one can obtain the Fourier transform of the quantum force autocorrelation function by multiplying the Fourier transform of the classical force autocorrelation function by a frequency-dependent factor in the case that the system mode is linearly coupled with a harmonic bath. This factor accounts for the quantum thermal population of the bath. However, due to the nonlinear coupling [20] between the Si-H stretching mode, the bending modes, and the surface phonons, the classical force autocorrelation function is employed in the Redfield equations in the present calculations since it is not clear how to generalize the Bader-Berne procedure. A newly proposed quantum dynamics method [ 38-41 ] should provide a better approximation to the exact quantum correlation function than the classical limit, particularly at low temperatures. Such an approach will be the focus of future research.
3. Car-Parrinelio simulation: geometry and frequencies of vibrational modes on the stepped surfaces As discussed in the previous study [20], the dominant vibrational energy relaxation pathway of the Sill stretching modes on the flat H/Si( 111 ) 1 × 1 surface relies on three Si-Si-H bending phonons and one surface phonon to establish resonance [5,20]. Because of this, the value of the bending frequency has a significant effect on the lifetime of the stretching mode on the terrace. The Si-Si-H bending frequency on the terrace of the flat H / S i ( 111 ) 1 x 1 surface was found to be 640 c m - I with electron energy loss spectroscopy [ 23,24]. Since the hydrogen atoms at the steps on the monohydride stepped surface terminate the surface at
the tetrahedral positions in the same way as the hydrogen atoms on the terrace, the Si-Si-H bending motion is therefore expected to play a similar key role in the relaxation process of the step Si-H stretching mode. Besides the bending motion, the equilibrium geometry of the stepped monohydride H / S i ( 111 ) 1 x 1 surface could also have a significant effect on the power spectrum of the fluctuating force along the Si-H bond and must be examined as well. In order to probe the bending frequency and geometry of the monohydride step, Car-Parrinello molecular dynamics (CP-MD) simulations of two miscut stepped surfaces were carried out. It has been shown that many dynamical and structural characteristics of the silicon surfaces can be calculated using this ab initio simulation method [28-31,42]. For instance, the features of the de-reconstruction of H / S i ( 111 )2× 1 surface [28] have been studied using the CP-MD method. Moreover, the nonlinear vibrational mode coupling on the flat H/Si( 111 ) 1 × 1 surface have been studied [ 32]. In the present calculations, the wave functions of the valence electrons were expanded in a plane wave basis set. The simulation supercell was periodically replicated in the x, y, and z directions. The potential energy between the nuclei and electrons was calculated using the nonlocal, norm-conserving pseudopotentials of Bachelet et al. [43] The Kleinman-Bylander [44] approach was employed for the nonlocal terms for both Si and H atoms. The Perdew-Zunger parameters for the exchange-correlation energy expression [ 45] which was formulated by Ceperley and Alder [46] was calculated within the local density approximation. The interaction between the nuclei was calculated with the standard Ewald sum energy approach [47]. An Ewald convergence parameter [48] of r/ = 1.0 was determined to be sufficient for the study of the stretching and bending mode frequencies. Due to the two-dimensional periodicity of the the substrate, the monohydride stepped surfaces needed to be constructed in a proper way to fit the twodimensional periodic boundary conditions in the directions defined by x t and y' directions in Fig. 1 (the z ~ direction is normal to the slab surface). The procedure to construct the 9 ° cut monohydride stepped H / S i ( I I 1 ) surface is described as follows: First, a slab of eight layers of silicon, which models the substrate, was placed in the simulation cell at bulk crys-
Y.-C. Sun et al. / Chemical Physics 200 (1995) 357-368
Fig. 1. The coordinate system (top view) for the flat H/Si( I 11 ) 1x 1 surface. The x I and yl directions are shown, while the z t axis is normal to the slab surface. The filled circles are hydrogen atoms, while the open circles are silicon atoms.
Z
~
. . . . . . . . . . . . . . . . . . . . .
l
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
,~y,
.
Fig. 2. A side view of the etched silicon surface where about half of the silicon atoms in the first two layers are removed for both sides of the slab. The dashed line shows the surface before etching. tal positions. Second, the surface o f the substrate was "etched" as shown in Fig. 2. A b o u t half of the silicon atoms o f the first two layers on both sides of the model slab were removed. The etching process results in a step height of 3.14/~. Third, the Si substrate was rotated 10.0 ° about the x ~ axis to satisfy the periodic boundary conditions in the x and y directions as shown by the dashed line in Fig. 3. Finally, a twodimensional rectangular periodic boundary condition which had the length of 18.0 /~ along the y direction was applied in the x - y - z coordinate system. The multiplicity o f 3.84 A and 1 8 . 0 / k could be adopted for the extended periodic boundary condition in the x and y directions, respectively. The supercell then periodically duplicates itself in the x and y directions to model the 9 ° monohydride stepped H / S i ( 111 ) surface. The hydrogen atoms were added to the both sides
361
Fig. 3. A side view of the s~pped silicon surface. The etched surface shown in Fig. 2 is now rotated 10.0° about the x p axis. The dashed lines represent the boundary of the supercell. of the slab at the tetrahedral positions to terminate its top and bottom surfaces. A calculation o f a smaller superceii system was first carried out to examine the relaxed structure o f 9 ° surface. The supercell was set up to have a side length o f 7.68/~ and 18.0/~ in the x and y directions, respectively. Due to the periodicity of the supercell in the z direction, the inter-slab length was set to be 23.0/~ to reduce the interaction between two neighboring slabs to a negligible amount. All silicon atoms were placed in their crystal positions and the hydrogen atoms were placed so that the S i - H bonds had a bond length of 1.48/~. The supercell o f this model 9 ° monohydride surface contained 64 Si atoms and 24 H atoms with a terrace length of 17.7/~. The wave functions o f valence electrons were expressed in terms of the plane wave basis coefficients in k-space. The kinetic energy cutoff was set to be 10 Ry. All of the plane wave coefficients were initialized randomly. The wave functions of the ground state were found by first minimizing the electronic energy with all the nuclei fixed at their initial positions. The nuclei were then allowed to move to find their minimum configurations. The nuclei were also randomly displaced within a radius of 0 . 3 / ~ from these positions and allowed to relax again in order to preclude being trapped in a local minimum. Several features of the relaxed structure on the 9 ° H / S i ( l l l ) I x l surface are illustrated in Fig. 4. The bond length of the H-Si bond on the terrace was 1.55 A, while it was found to be 1.53/~ at the step. The equilibrium structure is slightly distorted to the direction into vacuum as shown in the Fig. 4. The displacement
362
Y.-C. Sun et al. / Chemical Physics 200 (1995) 357-368
1-53A/H
H
°
2.0
Fourier Transform of -,
1.5 .~1.0
0.5
Fig. 4. The characteristics of the relaxed structure of the stepped surface according to the Car-Parrinello calculations. from bulk crystal positions is within 0.02 (,~). Similar results were found in a study of H / S i ( 111 ) 2 × 1 surface [28]. To check convergence, the simulation superceil was then expanded along the y direction to include a total number of 176 atoms. The characteristics of the relaxed structure remained unchanged. When the kinetic energy cutoff for the electronic wave function coefficients was increased to 12 Ry, the S i - H bond length decreased to 1.52 and 1.51 ~, on the terrace and step respectively. In addition, the supercell length in the z direction was increased to reduce the interaction between two neighbor model substrates in the z-direction. The features of the relaxed structure remained unchanged. From these calculations, it was concluded that there is very little difference in the S i l l bond length between the terrace and the step. The vibrational frequencies of the modes on the surface were also examined by Fourier transforming a trajectory from the CP-MD simulation [32]. By applying a Fourier transform to the S i - H bond distance, one can determine the a n h a r m o n i c S i - H stretching frequency at a given temperature in accordance with the Heisenberg Correspondence Principle [32]. The bending frequency can likewise be obtained by applying a Fourier transform to the hydrogen atom momentum or velocity in the x and y directions. The CP-MD simulation was initiated from the relaxed structure in the supercell. The "coordinates" and "velocities" of the electronic coefficients which represent the wave functions in k-space were given their initial values from the aforementioned minimization procedure. The fictitious mass of the electronic wave functions was set to 400 a.u., and 0.19 fs was used for the integration
00
50 0
10'00
1500 .....
2000
(era-') Fig. 5. The Fourier transform of the Si-H bond distance coordinates at the step and on the terrace. The solid line is the spectrum for the step Si-H stretch while the dashed line is the spectrum of the terrace Si-H stretch. step size during the simulation. A short molecular dynamics run was first carried out to rescale the temperature of the nuclei to 300 K. The final condition of the equilibrated system in the superceil was thus saved as the initial condition of the subsequent CP-MD simulation. The latter consisted of a 4 ps trajectory to collect the M D data. The S i - H stretching mode at the step and on the terrace were first examined. The calculated spectrum of the S i - H bond distance is shown in Fig. 5. A peak at 1966 cm - j for the step S i - H stretching mode is shifted ~ 20 c m - ] to the blue of the peak for the terrace S i - H mode. In comparing with the measured stretching frequencies, which are 2089 and 2084 c m - J [ 15,16] for the step and terrace S i - H modes, respectively, the trend in the calculated results are consistent with the experimental results in that the step S i l l mode has higher frequency than the terrace S i - H mode. For the bending frequency, the calculated spectrum of the hydrogen coordinate in the y direction is shown in Fig. 6. The peak for the step bending mode at 590 cm -~ is ~ 20 c m - I higher than on the terrace. On the terrace, both the calculated stretching and bending frequencies exhibit an overall shift to the red of the measured results. This is consistent with a number of studies which have found that density functional theory within the local density approximation underestimates the vibrational frequencies for some systems [30,49,50]. As the measured terrace bending
363
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
4. Calculations and results: lifetimes of the Si-H stretching modes
'ourier Transform of <6Y(t)6Y(O)) 1.0 F
0.8 ~--.0.6 c~ 0.4 0.2 i
O0
250
500
750
1000 1250
500
(cm -1) Fig. 6. The Fourier transform of the hydrogen y coordinates at the step and on the terrace. The solid line is the spectrum of the step Si-Si-H bend, while the dashed line is the spectrum of the terrace Si-Si-H bend.
frequency is 640 cm -1, the actual Si-Si-H bending frequency at the step should be higher than 640 c m - l in view of the calculated trend. The empirical Si-SiH bending potential at the step which was used in the molecular dynamics simulation for the lifetime calculations was therefore obtained by appropriately scaling the relevant bending force constant so the step bending frequency was given by 640 cm-1 times a factor of 1.04. The latter factor is the ratio of the calculated bending frequencies at the step to on the terrace. It should be noted that the CP-MD procedure was necessary to help refine the empirical potential [20] since the experimental step bending frequency is not known. Besides the 9 ° surface, a 15 ° surface was also examined. The model 15° surface was constructed following the procedures described previously for the 9 ° surface, except the x axis rotation was changed to 15.4 °. After the periodic boundary condition were imposed, the terrace length was found to be 11.1 A. With side lengths of 7.68 and 11.5/~ along the x and y directions, respectively, the CP-MD simulation superceil contained 56 Si atoms and 16 H atoms. The same procedure was carried out to examine the relaxed structure of the 15 ° stepped surface as described earlier for the 9 ° surface. The features of the relaxed structure for the 15° surface were found to be the same as for the 9 ° surface.
With the refined empirical potential from the CPMD simulations in hand, a molecular dynamics simulation was carried out using the empirical potential to obtain the fluctuating force autocorrelation functions necessary for the Bloch-Redfield equations which, in turn, yield the relaxation rate of excited Si-H stretching modes on the surfaces. The lifetimes of the vibrationally excited states were calculated by monitoring the population evolution represented by the timedependent diagonal reduced density matrix elements. The time evolution of these elements is determined by the linearly coupled Redfield equations, Eq. (2.2). The coupling terms in Eq. (2.2) were obtained from the power spectrum of the fluctuating force along a particular Si-H bond. During the simulations, the Sill bond was constrained to its equilibrium bond length and the force along the bond was calculated according to the equation FSi-H =
/-/,Si-H [m~,'Fsi -
mH'FH] • nSi-H ,
(4.1)
where msi(mH) and Fsi(FH) are the mass and force on the Si (H) atom. The vector nSM-I is the unit vector along the Si-H bond coordinate, and the quantity/ZSi-H is the reduced mass of the Si-H oscillator. The empirical potential used in the simulations is described in Ref. [20] aside from the modifications described in the previous section and below. This potential is separable in the Si-H stretching and Si-Si-H bending bond-angle coordinates, but these motions are kinetically coupled [20]. The accuracy of this potential for describing the bulk and surface phonon spectrum is discussed in Ref. [20] and references cited therein. While the potential parameters on the terrace of the stepped surface were taken to be the same as in a previous study on the flat H/Si( 111 ) 1 x I surface [20], the step Si-H Morse parameters were chosen to reproduce the measured frequency, 2089 cm -I . The Si-Si-H bending potential at the step was scaled to give a frequency of 660 c m - J as discussed in Section 3. Partial charges of +0.014e and -0.014e were assigned to the hydrogen and first-layer silicon atoms to account for the dipole-dipole interaction between the Si-H oscillators. Due to the high frequencies of the Si-H stretching modes which range from 2070 to 2090 cm-J on
364
Y.-c. Sun et al./ Chemical Physics 200 (1995)357-368
the stepped surface, they decouple from the low frequency surface phonon. The B 1 mode, which is not observed in the one-color picosecond pump/probe SFG spectroscopy [ 15,16], may arise due to the coupling among vibrational modes. In order to examine the phonon coupling effect, a molecular dynamics simulation of 36 ps was carried out without constraining the Si-H bond length. The distances between the Si and H of several Si-H bonds were collected and their spectrum was obtained by applying a Fourier transform. The spectrum of the stretching mode at the step was found to have a single peak at 2089 c m - ' . The spectrum of two Si-H modes on the terrace were examined as well. The Si-H mode which is only one Si-Si bond away from the step Si-H mode and the Si-H mode at the middle of the terrace have the same peak at 2084 cm -j . From this result, it was concluded that there is no Si-H stretching phonon coupling in the classical simulation. The experimentally observed rise of the B1 mode on the monohydride surface may therefore be due to a quantum tunnel splitting. For the calculation of the Si-H stretch lifetime, the model superceli for the stepped H / S i ( I I 1 ) surface was set up in the same procedure described in Section 3, but the supercell size was taken to be larger than in the CP-MD simulations, containing 440 Si atoms and 120 H atoms. The Si atoms were initially placed at their crystal positions, and the H atoms were placed so as to reproduce the known Si-H bond length of 1.48 ~. The atomic velocities were selected from a Boltzmann distribution corresponding to 600 K, which is twice the desired equilibrium temperature of 300 K, and each trajectory was equilibrated for 2 ps to reach the target temperature. The trajectory information was then collected for 36 ps. The Si-H stretching mode at the step of the 9 ° monohydride H/Si( 111 ) surface was studied first. During these force calculations, one step Si-H bond was constrained to its equilibrium bond length. The force along the bond was then collected from the trajectories to calculate the power spectrum of the fluctuating force. A total of 40 independent trajectories were run and their spectra were averaged. A fourth-order Runge-Kutta-Gill integrator was used to integrate Newton's equations with the Si-H bond length fixed through a holonomic constraint. According to the Wiener-Khinchin theorem, one can directly calculate the power spectrum of the force instead of calculating the force-force autocorrelation
Fourier Transform of <$F(L)6F(O)> '1
103
,
r
101 10-1 i0
10-5 10-7 10-9
2089 . 500
0
,
i
,
1000 1500 2000 2500 3000 ( c m -1) F~7-
Fig. 7. The Fourier transform of the fluctuating force autocorrelation function projected along a step Si-H bond. The solid line is the position of the v - 1 --* 0 fundamental transition frequency, O)10
-- 2089
cm
-j
.
function and then applying the Fourier transform, i.e. [20], +oo
f
dte -i°~jt (6F(O)6F(t))cl
--00 2
= iim 1 +lfrr e_iO,ot~$F(t ) 7--.oo ~ dt
.,
(4.2)
To reduce the error due to using finite-time data, a fourterm Blackman-Harris [51 ] window was applied to the force before the Fourier transform. The calculated power spectra of the fluctuating force on the Si-H stretching modes at the step and on the terrace are shown in Figs. 7 and 8. Comparing the two, one notices that the splitting of the overtones of the bend (i.e., at ,,~ 1300 and ,~ 1950 cm -I in Fig. 8) is not seen for the step Si-H mode in Fig. 7. The splitting on the terrace is presumably due to the strong coupling between the Si-Si-H bend mode and a surface phonon. On the stepped surface, the nearest Si-H mode to the step Si-H mode is only one Si-Si bond away while all Si-H modes are separated by two Si-Si bonds on the terrace. The Si-Si-H bending motion at the step may not couple with the surface phonon as strongly as on the terrace because of this geometrical feature. As outlined earlier, the lifetime of the Si-H stretching modes can be obtained by solving the Redfield
Y.-C. Sun et al. I Chemical Physics 200 (1995) 357-368
Fourier Transform of <6F L)6F(O)> -- , ,
103
10 - I
~
lO 1o-5 1°-7' 10-90-
2094, 500
j
10'00 1500 20'00 25'00 3000
(cm F,~ J Fig. 8. The Fourier transform of the fluctuating force autocorrelation function projected along a terrace S i - H bond. The solid line is the position o f the v = I ~ 0 fundamental transition frequency, oJm - 2084 cm - I .
equations, Eq. (2.2). Since these equations are linear, the eigensolution of the coupling coefficient matrix and the initial condition determine the evolution of all the density matrix elements. Although all the density matrix elements are coupled, the eigensolution of a sixteen-by-sixteen matrix was found to be sufficient for converging the time evolution of the first excited state probability. From the slope of the v = 1 state probability versus time one finds that the curve is very close to an exponential decay with a lifetime of 0.7 us. The exponential decay of the first excited state population is mainly due to one term in the Redfield matrix, given by [20] 21(01qll)l 2 F+ool + r~-om = h2(1 + e-~,o,o) +oo
x
I --0(1
dte -"°'°' ( 6 F ( O ) 6 F ( t ) } d ,
•
(4.3)
the inverse of which is the lifetime. The lifetime of the Si-H stretching mode on the terrace was calculated as well and found to be 1.8 ns which is essentially the same as the calculated Si-H stretching mode lifetime for the fiat H / S i ( 111 ) 1 × 1 surface [20]. The appearance of the steps on the H / S i ( 111 ) 1 x 1 surfaces therefore does not change the lifetime of the terrace Si-H modes significantly.
365
As described previously, the Si atoms in the simulation were given the crystal equilibrium geometry. In order to examine the effect of surface structural relaxation on the vibrational relaxation rate of the Sill stretching modes, the slightly relaxed structure of the 9 ° monohydride H / S i ( 111 ) 1 x 1 surface obtained from the Car-Pamnello molecular dynamics calculations was used for the equilibrium positions of the H and Si atoms. The lifetime of the step Si-H mode was found to differ from the lifetime for the unrelaxed structure by only 0.1 us. The lifetime of the terrace mode did not change by any significant amount. In the one-color spectroscopy of the monohydride stepped H/Si( 111 ) 1 x 1 surface, the step and terrace Si-H stretching modes were found to have frequencies of 2089 and 2084 cm - t , respectively [ 15,16]. In addition to these two modes, another peak at 2071 c m - t which was labeled the B 1 mode arises in the two-color excitation spectrum. Due to the weak classical coupling found among the stretching modes as discussed above, this mode might be plausibly assigned to Si-H stretching modes which are only one Si-Si bond away from the step Si-H stretching modes. As a test, these Si-H modes were assigned to have the measured frequency of the B1 mode, 2071 cm -1. The molecular dynamics simulations were then repeated to calculate the lifetime of the Si-H stretching modes at the step and on the terrace. Lifetimes of 0.7 and 1.7 ns were obtained for the step and terrace Si-H modes, respectively. The assignment of the B 1 mode to the nearby Si-H modes apparently does not significantly change the lifetimes of the step or terrace Si-H modes. In all cases of the 9 ° monohydride stepped H / S i ( 111 ) 1 x 1 surface studied in the calculations, the lifetime of the step Si-H mode was found to be 2.5 times shorter than for the terrace mode. The experimentally measured values [ 15,16] are 0.4 ns and 1.0 ns for the step and terrace Si-H mode lifetimes, respectively, Note that the experimental ratio is 2.5, in exact agreement with the calculated trend. In addition to 9 ° monohydride stepped surface, the vibrational relaxation rates of Si-H stretching modes on the 15 ° monohydride stepped H/Si( 111 ) 1 x I surface were also calculated at 300 K. The 15 ° surface can be characterized by its terrace length of 11.1 /~. A model 15 ° monohydride H / S i ( 111 ) 1 x 1 surface was constructed in the same way described in Section 3. In the molecular dynamics simulation with the em-
366
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
pirical potential for the 15 ° stepped surface, the superceli contained 340 Si atoms and 80 H atoms. The power spectrum for the step and terrace Si-H stretching modes were calculated using the same procedure as described previously for the 9 ° surface. The lifetimes for the step and terrace Si-H stretching modes were found to be 1.7 ns and 0.6 ns, respectively. For this surface, the step Si-H mode is therefore predicted to have a shorter lifetime relative to the terrace mode than for the 9 ° surface. However, the difference is very small.
5. Discussion The energy transfer mechanism of the Si-H stretching modes into the substrate on the flat H/Si( 111 ) 1 x 1 surface has been discussed in a previous study [20]. The vibrational energy of the excited Si-H stretching mode relaxes into the bulk with the assistance of three Si-Si-H bending phonons and one surface phonon. However, in the empirical potential there are no potential energy coupling terms between the stretching and bending modes [20]. The stretching modes are instead coupled to the adjacent bending modes through the effective mass matrix G [52], the kinetic energy being given by ½PtGP, where P represents the column matrix of the momenta conjugate to the corresponding internal coordinate. The elements of the G matrix, which are functions of internal bond distances angles, provide the nonlinear coupling between the stretching and its adjacent bending mode "kinetically" [ 20,53]. The analysis of the G matrix elements for the internal coordinate system uncovers the origin of the nonlinear coupling between the bending and stretching motions. For example, one of the G matrix terms ~(3) v~4, for a single Si-Si-H bend effective mass is given by [52] G~3) 4'~ = m H I r~2H + ms~lrs-~2Si
gas phase molecules [53] and the hydrogens on the H/Si( 111 ) 1 x I surface [5,8,20]. On the stepped monohydride H/Si( 111 ) 1 x 1 surface, the termination of the monohydride at the step is located at the tetrahedral position in a similar fashion to the terrace. The kinetic coupling between the Si-H stretching modes and their adjacent Si-Si-H bending modes is therefore expected to be similar. A comparison of the calculated power spectrum at the step (Fig. 7) to that on the terrace (Fig. 8) confirms the energy relaxation of the Si-H stretching mode at the step proceeds by the same mechanism as the stretching modes on the terrace. The Si-H stretching mode on the step of the monohydride stepped H/Si( 111 ) 1 x 1 surface is nearly in resonance with the second overtone of the bending mode, the difference being made up by a surface phonon. The shorter lifetime of step Si-H mode is because of the somewhat higher Si-Si-H bending frequency at the step. As illustrated in the previous study on the flat H / S i ( 11 i ) 1 x 1 surface [20], a higher SiSi-H bending frequency results in a shorter lifetime of the Si-H stretching mode because of a more favorable resonance condition. The calculated lifetimes of the step Si-H mode is shorter than the terrace Si-H mode. In the simulations described earlier, a Si-H bond was constrained to its equilibrium bond length in the computation of the force fluctuations along that coordinate. The constrained molecular dynamics simulation employed the Lagrange undetermined multiplier method. In this method, the force on the particles is calculated first in the potential field and then the constraint force is added to keep the bond length constrained. The undetermined multiplier, ,~, can be expressed exactly in terms of dynamical quantities at any instant along the trajectory [54]. After ,~ is evaluated, the constraint force on the H atom and Si atom can then be calculated. Following Eq. (4.1), it is found that the force along the Si-H bond can be expressed as
+ ms-~l (r~_2H + rs-~_2Si-- 2 c o s ~ b rSi-H - ' r Si-Si} -l ' (5.1) where rsi_ H and rsi_si are the Si-H and Si-Si bond length coordinates, respectively. Several other elements in the G matrix also provide coupling between the bonds and the bends. Some time ago, it was shown that this coupling is an important factor for the intramolecular energy relaxation mechanism in
fbond = Fp + Fc +/ZSi-H I/'l---~, r0
(5.2)
where Fp is the force in the potential field before adding the constraint force. The second term, Fc, is the "conservation" force fixing the bond length,/" is the time derivative of the bond vector, and r0 is the fixed equilibrium bond length. It is straightforward to show that Fc is equal and opposite to the force Fp. The
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
third term in Eq. (5.2) is the centrifugal force on the Si-H bond. The effect of bond length changes on the power spectrum of the force fluctuations was examined by artificially changing the equilibrium bond length of the Si-H oscillators to 1.6 ]k on the terrace with other system parameters remaining the same. The Si-H lifetime was found to increase to 4.2 ns. This longer lifetime can be understood from the expression of the force in Eq. (5.2) since the force is inversely proportional to the fixed bond length of the constrained bond. However, the slight difference of the actual Sill bond length at the step to that on the terrace is not a contributing factor to the different lifetimes.
6. Concluding remarks The vibrational energy lifetimes of excited Si-H modes on stepped monohydride H / S i ( 111 ) 1 × 1 surfaces have been calculated using the Bloch-Redfield relaxation theory combined with molecular dynamics simulation [20]. The Si-H lifetimes on both the steps and the terrace were calculated. The relaxation rate was calculated from the power spectrum of the fluctuating force along one Si-H bond in the molecular dynamics simulation. The calculated lifetimes of the Si-H stretching modes at the step were found to be shorter (0.7 ns) than on the terrace ( 1.8 ns) for both the 9 ° and 15 ° miscut monohydride surfaces. Both the magnitude and the trend of the calculated values agree well with the experimentally measured lifetimes of 0.4 and 1.0 ns for the step and terrace Si-H stretching modes, respectively. The shorter lifetime of the step Si-H stretching mode arises because the Si-SiH bending frequency at the step is higher than on the terrace according to the results of an ab initio CarParrinello simulation. The higher bending frequency results in a shorter lifetime for the Si-H modes due to a more favorable resonance condition. The coupling mechanism in the calculations arises from a kinetic coupling between the bending and stretching motions, and this mechanism appears to account for most of the relaxation behavior.
367
Acknowledgement This research was supported in part by the National Science Foundation Materials Research Laboratory (MRL) Program, under Grant No. DMR9120668. Additional support was provided through a David and Lucile Packard Fellowship in Science and Engineering. Many of the computations reported in this paper were carried out on the University of Pennsylvania MRL Central Facility computers. The authors thank Dr. Charles Ursenbach for his suggestions on the preparation of the manuscript.
References [1] G.A. Somorjai, Chemistry in two dimensions: surfaces (Cornell Univ. Press, Ithaca, 1981 ). [2] Y.J. Chabal, Surface Sci. Repts. 8 (1988) 211. [3] YJ. Chabal, G.S. Higashi and K. Raghavachari, J. Vac. Sci. Technol. A. 7 (1989) 2104. [4] E Guyot-Sionnest, E Dumas, Y.J. Chabal and G. Higashi, Phys. Rev. Letters 64 (1990) 2156. [5] P. Guyot-Sionnest, P. Dumas and Y.J. Chabal, J. Electron Spectry. Relat. Phenom. 54/55 (1990) 27. [6] P. Guyot-Sionnest, Phys. Rev. Letters 66 ( 1991 ) 1489. [7] P. Guyot-Sionnest, Phys. Rev. Letters 67 (1991) 2323. [8] P. Dumas, Y.J. Chabal and G.S. Higashi, Phys. Rev. Letters 65 (1990) 1124. [9] P. Jakob and Y.J. Chabal, J. Chem. Phys. 95 ( 1991 ) 2897. [ 10] P. Jakob, Y.J. Chabal, K. Raghavachari, P. Dumas and S.B. Christman, Surface Sci. 285 (1993) 251. [ 11 ] P. Jakob, Y.J. Chabal, K. Raghavachari and S.B. Christman, Phys. Rev. B 47 (1993) 6839. [ 12] P. Jakob, YJ. Chabal, K. Raghavachari, R.S. Becker and A.J. Becke, Surface Sci. 275 (1992) 407. [ 13] Y.J. Chabal, J. Mol. Struct. 292 (1993) 65. [14] M.A. Hines, Y.J. Chabal, T.D. Harris and A. Harris, Phys. Rev. Letters 71 (1993) 2280. [15] M. Morin, P. Jakob, N.J. Levions, Y.J. Chabal and A.L. Hams, J. Chem. Phys. 96 (1992) 6203. ll6] K. Kuhnke, M. Morin, P. Jakob, N.J. Levions, YJ. Chabal and A.L. Hams, J. Chem. Phys. 99 (1993) 6114. [17] Y.R. Shen, Nature 337 (1989) 519. [18] X.D. Zhu, H. Suhr and Y.R. Shen, Phys. Rev. B 35 (1987) 3047. [ 19] K. Kuhnke, A.L. Hams, Y.J. Chabal, P. Jakob and M, Morin, J. Chem. Phys. 100 (1994) 6896. [20] H. Gai and G.A. Voth, J. Chem. Phys. 99 (1993) 740. [21] EE. Figueirido and R.M. Levy, J. Chem. Phys. 97 (1992) 703. [22] J.M. Jean, R.A. Friesner and G.R. Fleming, J. Chem. Phys. 96 (1992) 5827. [23] H. Kobayashi, K. Edamoto, M. Onchi and M. Nishijima, J. Chem. Phys. 78 (1983) 7429.
368
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
[24] H. Froitzheim, U. KOhler and H. Lammering, Surface Sci. 149 (1985) 537. [25] K. Raghavachari, P. Jakob and Y.J. Chabal, Chem. Phys. Letters 206 (1993) 156. [26] X.-E Li, D. Vanderbilt and R.D. King-Smith, Phys. Rev. B 50 (1994) 4637. [27] R. Car and M. Parrinello, Phys. Rev. Letters 55 (1985) 2471. [28] E Ancilotto and A. Selloni, Phys. Rev. Letters 68 (1992) 2640. [29] E Ancilotto, W. Andreoni, A. Selloni, R. Car and M. Parrinello, Phys. Rev. Letters 65 (1990) 3148. [30] I. Stich, M.C. Payne, R.D. King-Smith, J.-S. Lin and L.J. Clarke, Phys. Rev. Letters 68 (1992) 1351. [31] K.D. Brommer, M. Needles, B.E Larson and J.D. Joannopoulos, Phys. Rev. Letters 68 (1992) 1355. [32] H. Gai and G.A. Voth, J. Chem. Phys. 101 (1994) 1734. [33] R.K. Wangsness and E Bloch, Phys. Rev. 89 (1953) 728. [34] A.G. Redlield, IBM J. Res. Devel. 1 (1957) 19. [35] C.E Slichter, Principles of magnetic resonance (Harper and Row, New York, 1963). [361 D. Oxtoby, Advan. Chem. Phys. 47 (1981) 487. [37] J.S. Bader and B.J. Berne, J. Chem. Phys. 100 (1994) 8359. [38] J. Cao and G.A. Voth, J. Chem. Phys. 99 (1993) 10070. [39] J. Cao and G.A. Voth, J. Chem. Phys. 100 (1994) 5106. [40] J. Cao and G.A. Voth, J. Chem. Phys. 101 (1994) 6157.
[41] J. Cao and G.A. Voth, J. Chem. Phys. 101 (1994) 6168. [42] K. Laasonen, private communication. [43] G.B. Bachelet, D.R. Hamann and M. Schluter, Phys. Rev. B 26 (1982) 4199. [44] L. Kleinman and D.M. Bylander, Phys. Rev. Letters 48 (1992) 1425. [45] J. Perdew and A. Zunger, Phys. Rev. B 23 ( 1981 ) 5048. [46] D.M. Cepedey and B.J. Alder, Phys. Rev. Letters 45 (1980) 566. [47] G. Gaili and A. Pasquarello, in: Computer simulation in chemical physics, NATO ASI Series C, Vol. 397, eds. M.P. Allen and D.J. Tildesley (Kluwer Academic Publishers, Dordrecht, 1993) pp. 285-287. [48] D.K. Remler and P.A. Madden, Mol. Phys. 70 (1990) 921. [49] X.-P. Li and D. Vanderbilt, Phys. Rev. Letters 69 (1992) 2543. [50] E. Kaxiras and J.D. Joannopoulos, Phys. Rev. B 37 (1988) 8842. [511 F.J. Harris, Proc. IEEE 66 (1978) 51. [52] E.B. Wilson Jr., J.C. Decius and P.C. Cross, Molecular vibrations (Dover, New York, 1980). [53] E.L. Sibert III, W.P. Reinhardt and J.T. Hynes, J. Chem. Phys. 81 (1984) 1115. [54] D.L. Thompson and L.M. Raft, in: Theory of chemical reaction dynamics, Vol. 3 (CRC Press, Boca Raton, 1984).
Chemical Physics 200 (1995) 357-368
Vibrational energy relaxation dynamics of Si-H stretching modes on stepped H/Si ( 111 ) 1 × 1 surfaces Y i n g - C h i e h S u n l, H u a d o n g G a i 2, G r e g o r y A . V o t h Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323, USA Received 12 May 1995
Abstract The vibrational energy relaxation rates of excited Si-H stretching modes on the monohydride step of miscut H/Si( I 11 ) 1 × 1 surfaces are calculated using Bloch-Redfield theory combined with classical molecular dynamics (MD) simulation. The structure and vibrational frequencies of the surface are first investigated using the Car-Parrinello ab initio MD method. The calculated Si-Si-H bending frequencies and relaxed structures are then used to refine the empirical potential for the classical MD simulations. The lifetime of the excited Si-H stretching mode at the step is found to be shorter than the modes on the terrace. Both the magnitude and the trend of the calculated results agree well with the experimental measurement on the 9 ° monohydride stepped surface. The vibrational relaxation rate of the Si-H stretching modes on the 15° monohydride stepped surface are also calculated and predicted to have a slightly shorter lifetime than for the 9 ° surface.
1. Introduction The area around steps on solid surfaces is one of the most interesting zones for many physical and chemical processes at such interfaces. For example, the steps on etched surfaces usually act as the nucleation sites in the crystal growth process. The appearance of the steps on the surface may break chemical bonds of adsorbates and produce new chemical species on the surface [ 1 ]. These processes can also be substantially influenced by the adsorbate dynamics. A specific case of the latter is the hydrogen dynamics at the step of
t Present address: School of Pharmacy, 926-S, Department of Pharmaceutical Chemistry, Universityof California, San Francisco, San Francisco, CA 94143, USA. 2 Present address: EnvironmentalMolecular Sciences Laboratory, Battelle/Pacific Northwest Laboratories, KI-90, Richland, WA 99352, USA.
etched H / S i ( 111 ) surfaces which may play a crucial role in the production of semiconductor devices. Experimentally, significant progress has been made to reveal the hydrogen dynamics on Si( 111 ) surfaces [ 2 - 1 6 ] . Many aspects, for instance, the vibrational frequencies, growth geometries etc., of the stepped H / S i ( 111 ) surfaces were uncovered by various experimental methods. [3,9-12,14] The studies of the H-terminated silicon surfaces with infrared spectroscopy have been summarized in a review article by Chabal [ 13]. Because of the development of sum frequency generation(SFG) spectroscopy (see, e.g., Refs. [ 17,18 ] ), the vibrational modes can be resolved in both the frequency and time domains. For example, the Si-H stretching frequency on the flat H-terminated surface was measured using SFG spectroscopy and its lifetime was found to be about 1.0 ns at 300 K [4]. Due to the existence of the steps on the etched Si( 111 ) surface, the involvement of additional energy
0301-0104/95/$09.50 Q 1995 Elsevier Science B.V. All fights reserved SSD! 0301-0104(95)00270-7
358
Y.-C. Sun et al./Chemical Physics 200 (1995) 357-368
transfer pathways complicates the vibrational energy relaxation processes of the excited Si-H stretching modes on the stepped H / S i ( I 11 ) surfaces. In recent years, experiments have been carried out [15,16] to probe the vibrations and energy transfer mechanisms on stepped H/Si( 111 ) 1 x I surfaces using pump/probe SFG spectroscopy. The frequencies of the Si-H stretching modes at the step and on the terrace, as well as their lifetimes, were determined. On the 9 ° stepped monohydride surface, the Si-H stretching mode at the step was found to be 2089 c m - J . In addition, the measurements with two-color pump/probe spectroscopy have provided extensive information on the energy flow pathways on these surfaces. The spectral effects due to the renormalization of vibrational modes after excitation have also been calculated in a recent paper by Kuhnke et al. [ 19]. The latter effects account for the unusual spectral transients observed in the pump-probe measurements. Although the picture of energy flow on the dihydride H-terminated H/Si( 111 ) 1 × 1 surface was mapped out using on a kinetic model [ 16], there remains an ambiguity for the relaxation rates for some energy transfer channels of the monohydride stepped H / S i ( 111 ) 1 x I surface due to the various possible effects from coupling to surface phonons. The energy relaxation mechanism of Si-H stretching modes on the flat H/Si( 111 ) 1 x 1 surface has been extensively examined [2-16]. The long lifetime is presumably related to the large frequency mismatch between the Si-H stretching modes and the rest of the phonon modes on the surface. The temperature dependence of the lifetime suggests that the dominant relaxation mechanism is a four-phonon process [5]. A previous theoretical study of this system [20] was based on the quantum Redfield theory of relaxation [21,22] and obtained good agreement between the measured and calculated lifetimes for the Si-H stretching mode on the flat H / S i ( 111 ) 1 x 1 surface. The relaxation rate of the single excited Si-H surface bond into the substrate was determined by calculating the population evolution of vibrational states which, in turn, are influenced by the coupling to the substrate phonons. Specifically, the fluctuating force from the substrate calculated from classical molecular dynamics (MD) and projected along the Si-H bond coordinate contains the information on the latter couplings. On the monohydride stepped H/Si( 111 ) 1 × 1 sur-
face, experimental evidence suggests that the hydrogen atoms at the step terminate the surface in a similar way to the hydrogen on the terrace [ 12]. The silicon atoms at the step are therefore tetrahedrallycoordinated. The direct application of our simulation approach [20] for the relaxation of Si-H stretching modes of flat surfaces is therefore reasonable, enabling us to calculate the Si-H vibrational energy lifetimes both at the step and on the terrace. On the fiat surface, a large portion of the Si-H stretch vibrational energy slowly relaxes into a state consisting of three Si-Si-H bending phonons plus a surface phonon. This behavior occurs since the Si-H stretching motion is kinetically coupled with the adjacent Si-Si-H bending mode [20]. Because of this coupling, the bending frequency has a critical effect in the stretching mode relaxation. While the Si-Si-H bending frequency on the flat H/Si( 111 ) 1 × 1 surface has been measured to be 640 cm - l [23,24], to our knowledge the bending frequency at the monohydride step has not been measured. It was therefore necessary to first estimate this bending frequency and refine our empirical potential [20] before the lifetime of the Si-H stretching mode at the step could be calculated using the BlochRedfield/MD methodology. Besides the bending frequency, the geometrical structure of the stepped surface may also have a significant effect on the relaxation rate of the Si-H stretching modes. While the structure of the dihydride stepped H / S i ( 111 ) 1 × 1 surface has been studied both experimentally and theoretically [ 12,25,26], the structure of the stepped monohydride H/Si( 111 ) 1 × 1 surface has not been examined by direct calculation. Fortunately, both the Si-SiH bending frequency at the monohydride step and the step structure can be examined using ab initio methods. Accordingly, the first-principles Car-Parrinello molecular dynamics (CP-MD) method [27], which has been widely used to obtain dynamical and structural information on Si(l 11) surfaces, [28-32] was first implemented to study the monohydride stepped H/Si( 111 ) 1 x I surface. (See also Ref. [26] for an ab initio study of the dihydride stepped H/Si( 111 ) 1 x 1 surface.) In our CP-MD calculations (cf. Section 3), the bending frequency and structure were obtained using the ab initio simulation, and then this information was used to refine the empirical potential used in the Bloch-Redfield/classical MD simulation method
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
[20]. The sections of this paper are organized as follows: Section 2 briefly describes the Bloch-Redfield/MD methodology as it was applied in the present calculations. Next, Section 3 contains the description of the Car-Parrinello molecular dynamics simulations to estimate the geometry of the stepped surface as well as the vibrational frequencies shifts at the step. The refinement of the empirical potential from the results of the Car-Parrinello method simulation and the resulting Bloch-Redfield/MD simulations are then described in Section 4. A discussion of the results is given in Section 5 and the factors which affect the relaxation rate of the Si-H stretching modes are analyzed. Concluding remarks are given in Section 6.
2. Theory The description of Bloch-Redfield theory [21,3335 ] as it applies to the relaxation of an excited stretching state on a surface [20] is given briefly as follows: The total Hamiltonian is written as H = Hs(q) +
Ha(Q) + Hint(q, Q ) ,
(2.1)
where q and H s ( q ) are the coordinate and Hamiltonian of the quantum Si-H oscillator, respectively. The remaining "bath" coordinates and Hamiltonian are represented by Q and HB(Q), respectively, while the interaction between the system and the bath is described by the term Hint(q, Q). Several assumptions are made in the derivation of the Redfield equations [ 34]. First, it is assumed that the interaction Hamiltonian can be written a s H i n t ( q , Q ) = qrF(Q). This assumption is reasonable since at 300 K the dominant energy transfer should proceed one quantum at a time. 3 The second assumption is that the bath remains close 3 This assumption amounts to the linear response assumption. Since the Si-H oscillators are somewhat anharmonic, coupling terms like q2G(Q) could in principle contribute to the relaxation rate since even for one-phonon processes the matrix element (0lq211) (nlG(Q)In') could be non-zero. Though it is difficult to estimate the size of such couplings, it is unlikely that they are large in the present case because the Si-H oscillator is nearly harmonic in the lower regions of the well (i.e., (01q211) is very small). Furthermore, if such terms were significant the present one-phonon linear response-based approach would be incapable of describing the relaxation of the v = 2 Si-H mode which is not the case [20].
359
to thermal equilibrium during the relaxation process. This is valid because the amount of the energy transferred from the excited stretching states to the bath is small relative to the total energy of the bath and equilibrium for the bath can be re-established quickly after a small perturbation. The third assumption is that the relaxation time of the oscillator is much longer than the force-force autocorrelation time of the bath. This condition can be justified in that the time scale of the correlation function in the bath is within the picosecond range which is much shorter than the nanosecond lifetime of the excited Si-H stretching mode on the H / S i ( 111 ) 1 x 1 surface. The dynamical information of the Si-H system coupled to the bath is contained in the time evolution of the reduced density matrix elements for the system. The evolution of the reduced density matrix elements is governed by a series of linearly coupled differential equations which can be derived under the above assumptions according to the Bloch-Redfield theory. These equations are given by
dpit (t) dt
Z Riljk ei(°~'J+°~kDtpjk(t) , j,k pit(t) is the reduced density
(2.2)
where matrix element of the Si-H oscillator stretching states Ii) and Il) in the interaction representation, and o)i and o)l are the eigenfrequencies of the li) and Il) states respectively. The Riljk t e r m s are given by
r
r
(2.3) where the coupling coefficients between the density matrix elements are given by
Fk+tij x
1 h 2 ( 1 + e#O~,j) +oo
/.
dte -'~°'jt
(klqll) (ilqlj)
(6F(O)6F(t))cl
(2.4)
--00
and
F~O = x
1 h2(l + e-#~*,)
/.
(klqll)(ilqlJ)
dte -'~°'a (~F(O)~F(t))c~.
--(X)
(2.5)
360
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
In the above equations, (klqll) is the transition matrix element between Si-H stretching states [k) and [l), t~kt is the Kronecker delta, and the dynamical function 8F(t) is the classical fluctuating force on the Si-H bond arising from the coupling to the other modes of the system. The use of a classical force autocorrelation function in Eq. (2.4) and Eq. (2.5) is an approximation to a fully quantum mechanical Redfield theory. A standard symmetrization procedure [ 36] is employed in the derivation of quantum Redfield equations and then the resulting correlation function is replaced with the classical correlation function. It has been shown in a recent study by Bader and Berne [37] that one can obtain the Fourier transform of the quantum force autocorrelation function by multiplying the Fourier transform of the classical force autocorrelation function by a frequency-dependent factor in the case that the system mode is linearly coupled with a harmonic bath. This factor accounts for the quantum thermal population of the bath. However, due to the nonlinear coupling [20] between the Si-H stretching mode, the bending modes, and the surface phonons, the classical force autocorrelation function is employed in the Redfield equations in the present calculations since it is not clear how to generalize the Bader-Berne procedure. A newly proposed quantum dynamics method [ 38-41 ] should provide a better approximation to the exact quantum correlation function than the classical limit, particularly at low temperatures. Such an approach will be the focus of future research.
3. Car-Parrinelio simulation: geometry and frequencies of vibrational modes on the stepped surfaces As discussed in the previous study [20], the dominant vibrational energy relaxation pathway of the Sill stretching modes on the flat H/Si( 111 ) 1 × 1 surface relies on three Si-Si-H bending phonons and one surface phonon to establish resonance [5,20]. Because of this, the value of the bending frequency has a significant effect on the lifetime of the stretching mode on the terrace. The Si-Si-H bending frequency on the terrace of the flat H / S i ( 111 ) 1 x 1 surface was found to be 640 c m - I with electron energy loss spectroscopy [ 23,24]. Since the hydrogen atoms at the steps on the monohydride stepped surface terminate the surface at
the tetrahedral positions in the same way as the hydrogen atoms on the terrace, the Si-Si-H bending motion is therefore expected to play a similar key role in the relaxation process of the step Si-H stretching mode. Besides the bending motion, the equilibrium geometry of the stepped monohydride H / S i ( 111 ) 1 x 1 surface could also have a significant effect on the power spectrum of the fluctuating force along the Si-H bond and must be examined as well. In order to probe the bending frequency and geometry of the monohydride step, Car-Parrinello molecular dynamics (CP-MD) simulations of two miscut stepped surfaces were carried out. It has been shown that many dynamical and structural characteristics of the silicon surfaces can be calculated using this ab initio simulation method [28-31,42]. For instance, the features of the de-reconstruction of H / S i ( 111 )2× 1 surface [28] have been studied using the CP-MD method. Moreover, the nonlinear vibrational mode coupling on the flat H/Si( 111 ) 1 × 1 surface have been studied [ 32]. In the present calculations, the wave functions of the valence electrons were expanded in a plane wave basis set. The simulation supercell was periodically replicated in the x, y, and z directions. The potential energy between the nuclei and electrons was calculated using the nonlocal, norm-conserving pseudopotentials of Bachelet et al. [43] The Kleinman-Bylander [44] approach was employed for the nonlocal terms for both Si and H atoms. The Perdew-Zunger parameters for the exchange-correlation energy expression [ 45] which was formulated by Ceperley and Alder [46] was calculated within the local density approximation. The interaction between the nuclei was calculated with the standard Ewald sum energy approach [47]. An Ewald convergence parameter [48] of r/ = 1.0 was determined to be sufficient for the study of the stretching and bending mode frequencies. Due to the two-dimensional periodicity of the the substrate, the monohydride stepped surfaces needed to be constructed in a proper way to fit the twodimensional periodic boundary conditions in the directions defined by x t and y' directions in Fig. 1 (the z ~ direction is normal to the slab surface). The procedure to construct the 9 ° cut monohydride stepped H / S i ( I I 1 ) surface is described as follows: First, a slab of eight layers of silicon, which models the substrate, was placed in the simulation cell at bulk crys-
Y.-C. Sun et al. / Chemical Physics 200 (1995) 357-368
Fig. 1. The coordinate system (top view) for the flat H/Si( I 11 ) 1x 1 surface. The x I and yl directions are shown, while the z t axis is normal to the slab surface. The filled circles are hydrogen atoms, while the open circles are silicon atoms.
Z
~
. . . . . . . . . . . . . . . . . . . . .
l
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
,~y,
.
Fig. 2. A side view of the etched silicon surface where about half of the silicon atoms in the first two layers are removed for both sides of the slab. The dashed line shows the surface before etching. tal positions. Second, the surface o f the substrate was "etched" as shown in Fig. 2. A b o u t half of the silicon atoms o f the first two layers on both sides of the model slab were removed. The etching process results in a step height of 3.14/~. Third, the Si substrate was rotated 10.0 ° about the x ~ axis to satisfy the periodic boundary conditions in the x and y directions as shown by the dashed line in Fig. 3. Finally, a twodimensional rectangular periodic boundary condition which had the length of 18.0 /~ along the y direction was applied in the x - y - z coordinate system. The multiplicity o f 3.84 A and 1 8 . 0 / k could be adopted for the extended periodic boundary condition in the x and y directions, respectively. The supercell then periodically duplicates itself in the x and y directions to model the 9 ° monohydride stepped H / S i ( 111 ) surface. The hydrogen atoms were added to the both sides
361
Fig. 3. A side view of the s~pped silicon surface. The etched surface shown in Fig. 2 is now rotated 10.0° about the x p axis. The dashed lines represent the boundary of the supercell. of the slab at the tetrahedral positions to terminate its top and bottom surfaces. A calculation o f a smaller superceii system was first carried out to examine the relaxed structure o f 9 ° surface. The supercell was set up to have a side length o f 7.68/~ and 18.0/~ in the x and y directions, respectively. Due to the periodicity of the supercell in the z direction, the inter-slab length was set to be 23.0/~ to reduce the interaction between two neighboring slabs to a negligible amount. All silicon atoms were placed in their crystal positions and the hydrogen atoms were placed so that the S i - H bonds had a bond length of 1.48/~. The supercell o f this model 9 ° monohydride surface contained 64 Si atoms and 24 H atoms with a terrace length of 17.7/~. The wave functions o f valence electrons were expressed in terms of the plane wave basis coefficients in k-space. The kinetic energy cutoff was set to be 10 Ry. All of the plane wave coefficients were initialized randomly. The wave functions of the ground state were found by first minimizing the electronic energy with all the nuclei fixed at their initial positions. The nuclei were then allowed to move to find their minimum configurations. The nuclei were also randomly displaced within a radius of 0 . 3 / ~ from these positions and allowed to relax again in order to preclude being trapped in a local minimum. Several features of the relaxed structure on the 9 ° H / S i ( l l l ) I x l surface are illustrated in Fig. 4. The bond length of the H-Si bond on the terrace was 1.55 A, while it was found to be 1.53/~ at the step. The equilibrium structure is slightly distorted to the direction into vacuum as shown in the Fig. 4. The displacement
362
Y.-C. Sun et al. / Chemical Physics 200 (1995) 357-368
1-53A/H
H
°
2.0
Fourier Transform of -,
1.5 .~1.0
0.5
Fig. 4. The characteristics of the relaxed structure of the stepped surface according to the Car-Parrinello calculations. from bulk crystal positions is within 0.02 (,~). Similar results were found in a study of H / S i ( 111 ) 2 × 1 surface [28]. To check convergence, the simulation superceil was then expanded along the y direction to include a total number of 176 atoms. The characteristics of the relaxed structure remained unchanged. When the kinetic energy cutoff for the electronic wave function coefficients was increased to 12 Ry, the S i - H bond length decreased to 1.52 and 1.51 ~, on the terrace and step respectively. In addition, the supercell length in the z direction was increased to reduce the interaction between two neighbor model substrates in the z-direction. The features of the relaxed structure remained unchanged. From these calculations, it was concluded that there is very little difference in the S i l l bond length between the terrace and the step. The vibrational frequencies of the modes on the surface were also examined by Fourier transforming a trajectory from the CP-MD simulation [32]. By applying a Fourier transform to the S i - H bond distance, one can determine the a n h a r m o n i c S i - H stretching frequency at a given temperature in accordance with the Heisenberg Correspondence Principle [32]. The bending frequency can likewise be obtained by applying a Fourier transform to the hydrogen atom momentum or velocity in the x and y directions. The CP-MD simulation was initiated from the relaxed structure in the supercell. The "coordinates" and "velocities" of the electronic coefficients which represent the wave functions in k-space were given their initial values from the aforementioned minimization procedure. The fictitious mass of the electronic wave functions was set to 400 a.u., and 0.19 fs was used for the integration
00
50 0
10'00
1500 .....
2000
(era-') Fig. 5. The Fourier transform of the Si-H bond distance coordinates at the step and on the terrace. The solid line is the spectrum for the step Si-H stretch while the dashed line is the spectrum of the terrace Si-H stretch. step size during the simulation. A short molecular dynamics run was first carried out to rescale the temperature of the nuclei to 300 K. The final condition of the equilibrated system in the superceil was thus saved as the initial condition of the subsequent CP-MD simulation. The latter consisted of a 4 ps trajectory to collect the M D data. The S i - H stretching mode at the step and on the terrace were first examined. The calculated spectrum of the S i - H bond distance is shown in Fig. 5. A peak at 1966 cm - j for the step S i - H stretching mode is shifted ~ 20 c m - ] to the blue of the peak for the terrace S i - H mode. In comparing with the measured stretching frequencies, which are 2089 and 2084 c m - J [ 15,16] for the step and terrace S i - H modes, respectively, the trend in the calculated results are consistent with the experimental results in that the step S i l l mode has higher frequency than the terrace S i - H mode. For the bending frequency, the calculated spectrum of the hydrogen coordinate in the y direction is shown in Fig. 6. The peak for the step bending mode at 590 cm -~ is ~ 20 c m - I higher than on the terrace. On the terrace, both the calculated stretching and bending frequencies exhibit an overall shift to the red of the measured results. This is consistent with a number of studies which have found that density functional theory within the local density approximation underestimates the vibrational frequencies for some systems [30,49,50]. As the measured terrace bending
363
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
4. Calculations and results: lifetimes of the Si-H stretching modes
'ourier Transform of <6Y(t)6Y(O)) 1.0 F
0.8 ~--.0.6 c~ 0.4 0.2 i
O0
250
500
750
1000 1250
500
(cm -1) Fig. 6. The Fourier transform of the hydrogen y coordinates at the step and on the terrace. The solid line is the spectrum of the step Si-Si-H bend, while the dashed line is the spectrum of the terrace Si-Si-H bend.
frequency is 640 cm -1, the actual Si-Si-H bending frequency at the step should be higher than 640 c m - l in view of the calculated trend. The empirical Si-SiH bending potential at the step which was used in the molecular dynamics simulation for the lifetime calculations was therefore obtained by appropriately scaling the relevant bending force constant so the step bending frequency was given by 640 cm-1 times a factor of 1.04. The latter factor is the ratio of the calculated bending frequencies at the step to on the terrace. It should be noted that the CP-MD procedure was necessary to help refine the empirical potential [20] since the experimental step bending frequency is not known. Besides the 9 ° surface, a 15 ° surface was also examined. The model 15° surface was constructed following the procedures described previously for the 9 ° surface, except the x axis rotation was changed to 15.4 °. After the periodic boundary condition were imposed, the terrace length was found to be 11.1 A. With side lengths of 7.68 and 11.5/~ along the x and y directions, respectively, the CP-MD simulation superceil contained 56 Si atoms and 16 H atoms. The same procedure was carried out to examine the relaxed structure of the 15 ° stepped surface as described earlier for the 9 ° surface. The features of the relaxed structure for the 15° surface were found to be the same as for the 9 ° surface.
With the refined empirical potential from the CPMD simulations in hand, a molecular dynamics simulation was carried out using the empirical potential to obtain the fluctuating force autocorrelation functions necessary for the Bloch-Redfield equations which, in turn, yield the relaxation rate of excited Si-H stretching modes on the surfaces. The lifetimes of the vibrationally excited states were calculated by monitoring the population evolution represented by the timedependent diagonal reduced density matrix elements. The time evolution of these elements is determined by the linearly coupled Redfield equations, Eq. (2.2). The coupling terms in Eq. (2.2) were obtained from the power spectrum of the fluctuating force along a particular Si-H bond. During the simulations, the Sill bond was constrained to its equilibrium bond length and the force along the bond was calculated according to the equation FSi-H =
/-/,Si-H [m~,'Fsi -
mH'FH] • nSi-H ,
(4.1)
where msi(mH) and Fsi(FH) are the mass and force on the Si (H) atom. The vector nSM-I is the unit vector along the Si-H bond coordinate, and the quantity/ZSi-H is the reduced mass of the Si-H oscillator. The empirical potential used in the simulations is described in Ref. [20] aside from the modifications described in the previous section and below. This potential is separable in the Si-H stretching and Si-Si-H bending bond-angle coordinates, but these motions are kinetically coupled [20]. The accuracy of this potential for describing the bulk and surface phonon spectrum is discussed in Ref. [20] and references cited therein. While the potential parameters on the terrace of the stepped surface were taken to be the same as in a previous study on the flat H/Si( 111 ) 1 x I surface [20], the step Si-H Morse parameters were chosen to reproduce the measured frequency, 2089 cm -I . The Si-Si-H bending potential at the step was scaled to give a frequency of 660 c m - J as discussed in Section 3. Partial charges of +0.014e and -0.014e were assigned to the hydrogen and first-layer silicon atoms to account for the dipole-dipole interaction between the Si-H oscillators. Due to the high frequencies of the Si-H stretching modes which range from 2070 to 2090 cm-J on
364
Y.-c. Sun et al./ Chemical Physics 200 (1995)357-368
the stepped surface, they decouple from the low frequency surface phonon. The B 1 mode, which is not observed in the one-color picosecond pump/probe SFG spectroscopy [ 15,16], may arise due to the coupling among vibrational modes. In order to examine the phonon coupling effect, a molecular dynamics simulation of 36 ps was carried out without constraining the Si-H bond length. The distances between the Si and H of several Si-H bonds were collected and their spectrum was obtained by applying a Fourier transform. The spectrum of the stretching mode at the step was found to have a single peak at 2089 c m - ' . The spectrum of two Si-H modes on the terrace were examined as well. The Si-H mode which is only one Si-Si bond away from the step Si-H mode and the Si-H mode at the middle of the terrace have the same peak at 2084 cm -j . From this result, it was concluded that there is no Si-H stretching phonon coupling in the classical simulation. The experimentally observed rise of the B1 mode on the monohydride surface may therefore be due to a quantum tunnel splitting. For the calculation of the Si-H stretch lifetime, the model superceli for the stepped H / S i ( I I 1 ) surface was set up in the same procedure described in Section 3, but the supercell size was taken to be larger than in the CP-MD simulations, containing 440 Si atoms and 120 H atoms. The Si atoms were initially placed at their crystal positions, and the H atoms were placed so as to reproduce the known Si-H bond length of 1.48 ~. The atomic velocities were selected from a Boltzmann distribution corresponding to 600 K, which is twice the desired equilibrium temperature of 300 K, and each trajectory was equilibrated for 2 ps to reach the target temperature. The trajectory information was then collected for 36 ps. The Si-H stretching mode at the step of the 9 ° monohydride H/Si( 111 ) surface was studied first. During these force calculations, one step Si-H bond was constrained to its equilibrium bond length. The force along the bond was then collected from the trajectories to calculate the power spectrum of the fluctuating force. A total of 40 independent trajectories were run and their spectra were averaged. A fourth-order Runge-Kutta-Gill integrator was used to integrate Newton's equations with the Si-H bond length fixed through a holonomic constraint. According to the Wiener-Khinchin theorem, one can directly calculate the power spectrum of the force instead of calculating the force-force autocorrelation
Fourier Transform of <$F(L)6F(O)> '1
103
,
r
101 10-1 i0
10-5 10-7 10-9
2089 . 500
0
,
i
,
1000 1500 2000 2500 3000 ( c m -1) F~7-
Fig. 7. The Fourier transform of the fluctuating force autocorrelation function projected along a step Si-H bond. The solid line is the position of the v - 1 --* 0 fundamental transition frequency, O)10
-- 2089
cm
-j
.
function and then applying the Fourier transform, i.e. [20], +oo
f
dte -i°~jt (6F(O)6F(t))cl
--00 2
= iim 1 +lfrr e_iO,ot~$F(t ) 7--.oo ~ dt
.,
(4.2)
To reduce the error due to using finite-time data, a fourterm Blackman-Harris [51 ] window was applied to the force before the Fourier transform. The calculated power spectra of the fluctuating force on the Si-H stretching modes at the step and on the terrace are shown in Figs. 7 and 8. Comparing the two, one notices that the splitting of the overtones of the bend (i.e., at ,,~ 1300 and ,~ 1950 cm -I in Fig. 8) is not seen for the step Si-H mode in Fig. 7. The splitting on the terrace is presumably due to the strong coupling between the Si-Si-H bend mode and a surface phonon. On the stepped surface, the nearest Si-H mode to the step Si-H mode is only one Si-Si bond away while all Si-H modes are separated by two Si-Si bonds on the terrace. The Si-Si-H bending motion at the step may not couple with the surface phonon as strongly as on the terrace because of this geometrical feature. As outlined earlier, the lifetime of the Si-H stretching modes can be obtained by solving the Redfield
Y.-C. Sun et al. I Chemical Physics 200 (1995) 357-368
Fourier Transform of <6F L)6F(O)> -- , ,
103
10 - I
~
lO 1o-5 1°-7' 10-90-
2094, 500
j
10'00 1500 20'00 25'00 3000
(cm F,~ J Fig. 8. The Fourier transform of the fluctuating force autocorrelation function projected along a terrace S i - H bond. The solid line is the position o f the v = I ~ 0 fundamental transition frequency, oJm - 2084 cm - I .
equations, Eq. (2.2). Since these equations are linear, the eigensolution of the coupling coefficient matrix and the initial condition determine the evolution of all the density matrix elements. Although all the density matrix elements are coupled, the eigensolution of a sixteen-by-sixteen matrix was found to be sufficient for converging the time evolution of the first excited state probability. From the slope of the v = 1 state probability versus time one finds that the curve is very close to an exponential decay with a lifetime of 0.7 us. The exponential decay of the first excited state population is mainly due to one term in the Redfield matrix, given by [20] 21(01qll)l 2 F+ool + r~-om = h2(1 + e-~,o,o) +oo
x
I --0(1
dte -"°'°' ( 6 F ( O ) 6 F ( t ) } d ,
•
(4.3)
the inverse of which is the lifetime. The lifetime of the Si-H stretching mode on the terrace was calculated as well and found to be 1.8 ns which is essentially the same as the calculated Si-H stretching mode lifetime for the fiat H / S i ( 111 ) 1 × 1 surface [20]. The appearance of the steps on the H / S i ( 111 ) 1 x 1 surfaces therefore does not change the lifetime of the terrace Si-H modes significantly.
365
As described previously, the Si atoms in the simulation were given the crystal equilibrium geometry. In order to examine the effect of surface structural relaxation on the vibrational relaxation rate of the Sill stretching modes, the slightly relaxed structure of the 9 ° monohydride H / S i ( 111 ) 1 x 1 surface obtained from the Car-Pamnello molecular dynamics calculations was used for the equilibrium positions of the H and Si atoms. The lifetime of the step Si-H mode was found to differ from the lifetime for the unrelaxed structure by only 0.1 us. The lifetime of the terrace mode did not change by any significant amount. In the one-color spectroscopy of the monohydride stepped H/Si( 111 ) 1 x 1 surface, the step and terrace Si-H stretching modes were found to have frequencies of 2089 and 2084 cm - t , respectively [ 15,16]. In addition to these two modes, another peak at 2071 c m - t which was labeled the B 1 mode arises in the two-color excitation spectrum. Due to the weak classical coupling found among the stretching modes as discussed above, this mode might be plausibly assigned to Si-H stretching modes which are only one Si-Si bond away from the step Si-H stretching modes. As a test, these Si-H modes were assigned to have the measured frequency of the B1 mode, 2071 cm -1. The molecular dynamics simulations were then repeated to calculate the lifetime of the Si-H stretching modes at the step and on the terrace. Lifetimes of 0.7 and 1.7 ns were obtained for the step and terrace Si-H modes, respectively. The assignment of the B 1 mode to the nearby Si-H modes apparently does not significantly change the lifetimes of the step or terrace Si-H modes. In all cases of the 9 ° monohydride stepped H / S i ( 111 ) 1 x 1 surface studied in the calculations, the lifetime of the step Si-H mode was found to be 2.5 times shorter than for the terrace mode. The experimentally measured values [ 15,16] are 0.4 ns and 1.0 ns for the step and terrace Si-H mode lifetimes, respectively, Note that the experimental ratio is 2.5, in exact agreement with the calculated trend. In addition to 9 ° monohydride stepped surface, the vibrational relaxation rates of Si-H stretching modes on the 15 ° monohydride stepped H/Si( 111 ) 1 x I surface were also calculated at 300 K. The 15 ° surface can be characterized by its terrace length of 11.1 /~. A model 15 ° monohydride H / S i ( 111 ) 1 x 1 surface was constructed in the same way described in Section 3. In the molecular dynamics simulation with the em-
366
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
pirical potential for the 15 ° stepped surface, the superceli contained 340 Si atoms and 80 H atoms. The power spectrum for the step and terrace Si-H stretching modes were calculated using the same procedure as described previously for the 9 ° surface. The lifetimes for the step and terrace Si-H stretching modes were found to be 1.7 ns and 0.6 ns, respectively. For this surface, the step Si-H mode is therefore predicted to have a shorter lifetime relative to the terrace mode than for the 9 ° surface. However, the difference is very small.
5. Discussion The energy transfer mechanism of the Si-H stretching modes into the substrate on the flat H/Si( 111 ) 1 x 1 surface has been discussed in a previous study [20]. The vibrational energy of the excited Si-H stretching mode relaxes into the bulk with the assistance of three Si-Si-H bending phonons and one surface phonon. However, in the empirical potential there are no potential energy coupling terms between the stretching and bending modes [20]. The stretching modes are instead coupled to the adjacent bending modes through the effective mass matrix G [52], the kinetic energy being given by ½PtGP, where P represents the column matrix of the momenta conjugate to the corresponding internal coordinate. The elements of the G matrix, which are functions of internal bond distances angles, provide the nonlinear coupling between the stretching and its adjacent bending mode "kinetically" [ 20,53]. The analysis of the G matrix elements for the internal coordinate system uncovers the origin of the nonlinear coupling between the bending and stretching motions. For example, one of the G matrix terms ~(3) v~4, for a single Si-Si-H bend effective mass is given by [52] G~3) 4'~ = m H I r~2H + ms~lrs-~2Si
gas phase molecules [53] and the hydrogens on the H/Si( 111 ) 1 x I surface [5,8,20]. On the stepped monohydride H/Si( 111 ) 1 x 1 surface, the termination of the monohydride at the step is located at the tetrahedral position in a similar fashion to the terrace. The kinetic coupling between the Si-H stretching modes and their adjacent Si-Si-H bending modes is therefore expected to be similar. A comparison of the calculated power spectrum at the step (Fig. 7) to that on the terrace (Fig. 8) confirms the energy relaxation of the Si-H stretching mode at the step proceeds by the same mechanism as the stretching modes on the terrace. The Si-H stretching mode on the step of the monohydride stepped H/Si( 111 ) 1 x 1 surface is nearly in resonance with the second overtone of the bending mode, the difference being made up by a surface phonon. The shorter lifetime of step Si-H mode is because of the somewhat higher Si-Si-H bending frequency at the step. As illustrated in the previous study on the flat H / S i ( 11 i ) 1 x 1 surface [20], a higher SiSi-H bending frequency results in a shorter lifetime of the Si-H stretching mode because of a more favorable resonance condition. The calculated lifetimes of the step Si-H mode is shorter than the terrace Si-H mode. In the simulations described earlier, a Si-H bond was constrained to its equilibrium bond length in the computation of the force fluctuations along that coordinate. The constrained molecular dynamics simulation employed the Lagrange undetermined multiplier method. In this method, the force on the particles is calculated first in the potential field and then the constraint force is added to keep the bond length constrained. The undetermined multiplier, ,~, can be expressed exactly in terms of dynamical quantities at any instant along the trajectory [54]. After ,~ is evaluated, the constraint force on the H atom and Si atom can then be calculated. Following Eq. (4.1), it is found that the force along the Si-H bond can be expressed as
+ ms-~l (r~_2H + rs-~_2Si-- 2 c o s ~ b rSi-H - ' r Si-Si} -l ' (5.1) where rsi_ H and rsi_si are the Si-H and Si-Si bond length coordinates, respectively. Several other elements in the G matrix also provide coupling between the bonds and the bends. Some time ago, it was shown that this coupling is an important factor for the intramolecular energy relaxation mechanism in
fbond = Fp + Fc +/ZSi-H I/'l---~, r0
(5.2)
where Fp is the force in the potential field before adding the constraint force. The second term, Fc, is the "conservation" force fixing the bond length,/" is the time derivative of the bond vector, and r0 is the fixed equilibrium bond length. It is straightforward to show that Fc is equal and opposite to the force Fp. The
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
third term in Eq. (5.2) is the centrifugal force on the Si-H bond. The effect of bond length changes on the power spectrum of the force fluctuations was examined by artificially changing the equilibrium bond length of the Si-H oscillators to 1.6 ]k on the terrace with other system parameters remaining the same. The Si-H lifetime was found to increase to 4.2 ns. This longer lifetime can be understood from the expression of the force in Eq. (5.2) since the force is inversely proportional to the fixed bond length of the constrained bond. However, the slight difference of the actual Sill bond length at the step to that on the terrace is not a contributing factor to the different lifetimes.
6. Concluding remarks The vibrational energy lifetimes of excited Si-H modes on stepped monohydride H / S i ( 111 ) 1 × 1 surfaces have been calculated using the Bloch-Redfield relaxation theory combined with molecular dynamics simulation [20]. The Si-H lifetimes on both the steps and the terrace were calculated. The relaxation rate was calculated from the power spectrum of the fluctuating force along one Si-H bond in the molecular dynamics simulation. The calculated lifetimes of the Si-H stretching modes at the step were found to be shorter (0.7 ns) than on the terrace ( 1.8 ns) for both the 9 ° and 15 ° miscut monohydride surfaces. Both the magnitude and the trend of the calculated values agree well with the experimentally measured lifetimes of 0.4 and 1.0 ns for the step and terrace Si-H stretching modes, respectively. The shorter lifetime of the step Si-H stretching mode arises because the Si-SiH bending frequency at the step is higher than on the terrace according to the results of an ab initio CarParrinello simulation. The higher bending frequency results in a shorter lifetime for the Si-H modes due to a more favorable resonance condition. The coupling mechanism in the calculations arises from a kinetic coupling between the bending and stretching motions, and this mechanism appears to account for most of the relaxation behavior.
367
Acknowledgement This research was supported in part by the National Science Foundation Materials Research Laboratory (MRL) Program, under Grant No. DMR9120668. Additional support was provided through a David and Lucile Packard Fellowship in Science and Engineering. Many of the computations reported in this paper were carried out on the University of Pennsylvania MRL Central Facility computers. The authors thank Dr. Charles Ursenbach for his suggestions on the preparation of the manuscript.
References [1] G.A. Somorjai, Chemistry in two dimensions: surfaces (Cornell Univ. Press, Ithaca, 1981 ). [2] Y.J. Chabal, Surface Sci. Repts. 8 (1988) 211. [3] YJ. Chabal, G.S. Higashi and K. Raghavachari, J. Vac. Sci. Technol. A. 7 (1989) 2104. [4] E Guyot-Sionnest, E Dumas, Y.J. Chabal and G. Higashi, Phys. Rev. Letters 64 (1990) 2156. [5] P. Guyot-Sionnest, P. Dumas and Y.J. Chabal, J. Electron Spectry. Relat. Phenom. 54/55 (1990) 27. [6] P. Guyot-Sionnest, Phys. Rev. Letters 66 ( 1991 ) 1489. [7] P. Guyot-Sionnest, Phys. Rev. Letters 67 (1991) 2323. [8] P. Dumas, Y.J. Chabal and G.S. Higashi, Phys. Rev. Letters 65 (1990) 1124. [9] P. Jakob and Y.J. Chabal, J. Chem. Phys. 95 ( 1991 ) 2897. [ 10] P. Jakob, Y.J. Chabal, K. Raghavachari, P. Dumas and S.B. Christman, Surface Sci. 285 (1993) 251. [ 11 ] P. Jakob, Y.J. Chabal, K. Raghavachari and S.B. Christman, Phys. Rev. B 47 (1993) 6839. [ 12] P. Jakob, YJ. Chabal, K. Raghavachari, R.S. Becker and A.J. Becke, Surface Sci. 275 (1992) 407. [ 13] Y.J. Chabal, J. Mol. Struct. 292 (1993) 65. [14] M.A. Hines, Y.J. Chabal, T.D. Harris and A. Harris, Phys. Rev. Letters 71 (1993) 2280. [15] M. Morin, P. Jakob, N.J. Levions, Y.J. Chabal and A.L. Hams, J. Chem. Phys. 96 (1992) 6203. ll6] K. Kuhnke, M. Morin, P. Jakob, N.J. Levions, YJ. Chabal and A.L. Hams, J. Chem. Phys. 99 (1993) 6114. [17] Y.R. Shen, Nature 337 (1989) 519. [18] X.D. Zhu, H. Suhr and Y.R. Shen, Phys. Rev. B 35 (1987) 3047. [ 19] K. Kuhnke, A.L. Hams, Y.J. Chabal, P. Jakob and M, Morin, J. Chem. Phys. 100 (1994) 6896. [20] H. Gai and G.A. Voth, J. Chem. Phys. 99 (1993) 740. [21] EE. Figueirido and R.M. Levy, J. Chem. Phys. 97 (1992) 703. [22] J.M. Jean, R.A. Friesner and G.R. Fleming, J. Chem. Phys. 96 (1992) 5827. [23] H. Kobayashi, K. Edamoto, M. Onchi and M. Nishijima, J. Chem. Phys. 78 (1983) 7429.
368
Y.-C. Sun et al./ Chemical Physics 200 (1995) 357-368
[24] H. Froitzheim, U. KOhler and H. Lammering, Surface Sci. 149 (1985) 537. [25] K. Raghavachari, P. Jakob and Y.J. Chabal, Chem. Phys. Letters 206 (1993) 156. [26] X.-E Li, D. Vanderbilt and R.D. King-Smith, Phys. Rev. B 50 (1994) 4637. [27] R. Car and M. Parrinello, Phys. Rev. Letters 55 (1985) 2471. [28] E Ancilotto and A. Selloni, Phys. Rev. Letters 68 (1992) 2640. [29] E Ancilotto, W. Andreoni, A. Selloni, R. Car and M. Parrinello, Phys. Rev. Letters 65 (1990) 3148. [30] I. Stich, M.C. Payne, R.D. King-Smith, J.-S. Lin and L.J. Clarke, Phys. Rev. Letters 68 (1992) 1351. [31] K.D. Brommer, M. Needles, B.E Larson and J.D. Joannopoulos, Phys. Rev. Letters 68 (1992) 1355. [32] H. Gai and G.A. Voth, J. Chem. Phys. 101 (1994) 1734. [33] R.K. Wangsness and E Bloch, Phys. Rev. 89 (1953) 728. [34] A.G. Redlield, IBM J. Res. Devel. 1 (1957) 19. [35] C.E Slichter, Principles of magnetic resonance (Harper and Row, New York, 1963). [361 D. Oxtoby, Advan. Chem. Phys. 47 (1981) 487. [37] J.S. Bader and B.J. Berne, J. Chem. Phys. 100 (1994) 8359. [38] J. Cao and G.A. Voth, J. Chem. Phys. 99 (1993) 10070. [39] J. Cao and G.A. Voth, J. Chem. Phys. 100 (1994) 5106. [40] J. Cao and G.A. Voth, J. Chem. Phys. 101 (1994) 6157.
[41] J. Cao and G.A. Voth, J. Chem. Phys. 101 (1994) 6168. [42] K. Laasonen, private communication. [43] G.B. Bachelet, D.R. Hamann and M. Schluter, Phys. Rev. B 26 (1982) 4199. [44] L. Kleinman and D.M. Bylander, Phys. Rev. Letters 48 (1992) 1425. [45] J. Perdew and A. Zunger, Phys. Rev. B 23 ( 1981 ) 5048. [46] D.M. Cepedey and B.J. Alder, Phys. Rev. Letters 45 (1980) 566. [47] G. Gaili and A. Pasquarello, in: Computer simulation in chemical physics, NATO ASI Series C, Vol. 397, eds. M.P. Allen and D.J. Tildesley (Kluwer Academic Publishers, Dordrecht, 1993) pp. 285-287. [48] D.K. Remler and P.A. Madden, Mol. Phys. 70 (1990) 921. [49] X.-P. Li and D. Vanderbilt, Phys. Rev. Letters 69 (1992) 2543. [50] E. Kaxiras and J.D. Joannopoulos, Phys. Rev. B 37 (1988) 8842. [511 F.J. Harris, Proc. IEEE 66 (1978) 51. [52] E.B. Wilson Jr., J.C. Decius and P.C. Cross, Molecular vibrations (Dover, New York, 1980). [53] E.L. Sibert III, W.P. Reinhardt and J.T. Hynes, J. Chem. Phys. 81 (1984) 1115. [54] D.L. Thompson and L.M. Raft, in: Theory of chemical reaction dynamics, Vol. 3 (CRC Press, Boca Raton, 1984).