Femtosecond-Laser Desorption of H2 (D2) from Ru(0001): Quantum

7790

J. Phys. Chem. C 2009, 113, 7790–7801

Femtosecond-Laser Desorption of H2 (D2) from Ru(0001): Quantum and Classical Approaches Tijo Vazhappilly, Tillmann Klamroth,* and Peter Saalfrank Institut fu¨r Chemie, UniVersita¨t Potsdam, Karl-Liebknecht-Strasse 24-25, D-14476 Potsdam-Golm, Germany

Rigoberto Hernandez Center for Computational Molecular Sciences and Technology, School of Chemistry & Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332-0430 ReceiVed: December 5, 2008; ReVised Manuscript ReceiVed: February 23, 2009

The femtosecond-laser-induced, substrate-mediated associative desorption of molecular hydrogen and deuterium from a Ru(0001) surface in the so-called DIMET limit is studied theoretically. Two widely used models, a “quantum nonadiabatic” approach and a “classical adiabatic” one are employed and compared to each other. The quantum model is realized by the Monte Carlo wave packet (MCWP) method in the framework of opensystem density matrix theory. The classical approach is realized with the help of (frictional) Langevin dynamics with stochastic forces. For both models the same ground-state potential energy surface is used and the same two-temperature model adopted to describe the hot-electron-driven desorption dynamics. Apart from these common features both models are different. Still, both account well for the main experimental findings (Wagner et al. Phys. ReV. B 2005, 72, 205404). In particular, an isotope effect in desorption probabilities, the energy content of the desorbing molecules, and the scaling of these observables with laser fluence are reproduced and explained. The similarity of the results obtained with both models is traced back to the fact that, in the present case, the photodynamics takes place dominantly in the ground electronic state because the electronically excited state is rapidly quenched. The short lifetime of the excited state has also the effect that photoreaction cross sections are typically very small. An IR+vis hybrid scheme, by which the adsorbate is vibrationally excited by IR photons prior to the heating of metal electrons with the vis pulse, is shown to successfully promote the reaction even for strongly coupled adsorbate-surface systems. I. Introduction The rate and mechanism of chemical reactions are known to be strongly affected by the proximity of the reactants to a surface. The canonical case occurs when the surface catalyzes the reaction, but the suppression of undesirable products is equally useful. In fact, any alteration in the relative abundance and rates of production of the products and their states (viz. exit channels) presents the possibility of control. The photochemical route for selective control of chemical reactivity at surfaces is a particularly promising one. In this context, the challenge to theoretical and computational chemistry is the prediction of the composition and rates of reactive adsorbate systems based on fundamental principles. Ideally, this requires the use of a direct and exact quantum mechanical treatment of all the reactants′ molecular degrees of freedom and the surface. In practice, such a system is too large to be treated exactly. Instead one must take recourse in the use of some combination of the following approximations: the surface is represented by a small and finite particle model, the surface is treated as a classical or quantum mechanical “bath”, the reactant molecules are treated in a reduced dimensional representation, some degrees of freedom are decoupled from the reaction because of a mismatch in time scales, and some of the components are described by classical equations of motion rather than the fully quantum mechanical one. The question to be explored in this work is the degree to which these approximations can be taken while still obtaining quantitative accuracy, or perhaps less stringently, qualitative insight. Specif-

ically, in this article, we study the femtosecond-laser-induced recombinative desorption of hydrogen or deuterium molecules from a Ru(0001) surface. This reaction has been studied by Frischkorn et al. using femtosecond laser techniques.1-3 Perhaps not surprisingly, they find that the photochemical reaction has a completely different mechanism than that which was seen in thermal desorption experiments. Among their key observations the following are relevant to this work: (i) the recombinative desorption has a large isotope effect with a desorption probability for H2 that is typically 10 times higher than the one for D2 at low fluences; (ii) desorbing molecules have high translational energy compared to vibrational energy; (iii) the desorption probability increases nonlinearly upon increasing the laser fluence; at the same time, the isotope effect in the desorption probabilities goes down and energy partitioning to different degrees of freedom changes; (iv) the reaction goes through a hot-electron-mediated mechanism with a short response time of ∼1 ps as seen in their two-pulse correlation experiments. Frischkorn et al.1-3 performed these experiments using 120-130 fs laser pulses at 800 nm center wavelength (1.55 eV) and different fluences. These femtosecond laser pulses initially heat the electrons of the ruthenium surface to produce high electron temperatures, Tel(t), which then transfer energy to the adsorbate-substrate bond to break it. Thus the visible (vis) excitation is indirect, i.e., substrate-mediated, and therefore hard to control. The photodesorption of adsorbates from a metal surface is mainly nonadiabatic in nature because of the involvement of a

10.1021/jp810709k CCC: $40.75  2009 American Chemical Society Published on Web 04/09/2009

Laser Desorption of H2 (D2) from Ru(0001) continuum of electronic excitations and interstate, nonadiabatic couplings. In an effective two-state model, laser-induced transitions take place between a ground state, |g〉, and an excited state, |e〉. The latter relaxes within an ultrashort lifetime τel due to coupling with the electron-hole pair continuum of the metal surface. Desorption occurs if enough energy is gained to break the surface-adsorbate bond. When femtosecond laser pulses are used with high fluences, there may be multiple electronic excitations and relaxations. As a result, a vibrational ladder climbing in the ground electronic state occurs. If the rate of this process is higher than the rate of vibrational relaxation, desorption may take place. This is called desorption induced by multiple electronic transitions (DIMET).4,5 DIMET leads to higher desorption yields compared to the DIET case (desorption induced by electronic transitions), which is realized under singleexcitation conditions and typically enforced by nanosecond lasers. In DIMET, one usually finds a power law dependence of the yield Y on the fluence F: Y ∝ Fn, where n is typically in the range, 2 e n e 10.6 In contrast, for DIET the yield increases linearly with F. In this paper, we investigate the performance and relation of two dielectic approaches in surface femtochemistry hinging on the use of either adiabatic or nonadiabatic representations in the treatment of electronic excited-state dynamics. In the nonadiabatic representation, at least part of the dynamics explicitly takes place in one or several excited states, while in the adiabatic case, nonadiabatic effects are only implicitly accounted for. Both representations can be realized either classically or quantum mechanically. In the classical adiabatic picture, molecular dynamics with electronic frictions can be treated in the framework of the Langevin formalism7,8 using a single (ground-state) potential energy surface. This approach will be used in section III.B. The method has already been chosen by Luntz et al. for photodesorption of H2 from ruthenium.9 The interaction of the laser pulse with metal electrons was treated indirectly within the twotemperature model.10 The excited electrons couple to the nuclear adsorbate degrees of freedom through random forces that drive the reaction. As an alternative, DIMET can also be treated within a nonadiabatic representation using ground and excited potential energy surface(s).5 In ref 5, the electronic excitation/deexcitation of classical trajectories was realized stochastically by assuming Franck-Condon transitions between two states |g〉 and |e〉. The electronic excitation and deexcitation rates were calculated by the principle of detailed balance from the time-dependent electronic temperature, Tel, and a phenomenological excitedstate lifetime, τel. This type of stochastic dynamics can also be realized in quantum dynamical fashion using stochastic or jumping wave packets,11,12 e.g., open-system density matrix theory.13 This approach is followed here. Beyond the question of which of these methods is accurate and/or practical, we also wish to use them to study mechanisms with which to control/enhance the photodesorption in H/Ru(0001). In particular, we wish to realize conditions under which we can control/enhance the all-too-often small photodesorption yields. For this purpose, we suggest a hybrid scheme to control the reaction in which the adsorbate vibrations are initially excited by an IR pulse prior to a vis pulse (or UV/vis pulse in general). One can preexcite certain adsorbate modesse.g., the adsorbatesurface bondsbecause IR photons couple directly with it. In contrast, UV/vis photons “heat” the metal surface instead. Thus, the IR excitation is optically controllable14 and prepares the system for the subsequent photoreaction by the UV/vis laser

J. Phys. Chem. C, Vol. 113, No. 18, 2009 7791

Figure 1. Ground-state potential Vg(r, Z). Contours start at 0.5 eV with an increment of 0.5 eV. Bullets indicate the minimum at Z0 ) 1.06 Å and r0 ) 2.75 Å and the transition state at Z‡ ) 2.24 Å and r‡ ) 0.77 Å. The dotted curve shows the approximate minimum energy path S for associative desorption.

pulse. Similar IR+UV/vis strategies to control the photodesorption or other reactions of molecules at surfaces have been suggested elsewhere.13-17 The strategy is inspired from “vibrationally mediated chemistry” in gas-phase dynamics.18-20 It seems to be much more difficult to attain such a goal for complex systems including condensed phases due to the dissipative environment. The paper is organized as follows. In section II, theoretical methods and models for quantum + nonadiabatic and classical + adiabatic descriptions are introduced. In the latter, molecular dynamics on a single potential with electronic frictions and random forces will be employed (section II.C). The former makes use of electronically excited states and open-system density matrix theory (section II.B). In both cases, a twodimensional ground-state potential is central as will be described in section II.A. Section III.A describes the quantum mechanical treatment of DIMET of H2/D2 from a Ru(0001) surface. Among other things, the nonlinear fluence dependence of the reaction yield on the incoming vis laser pulse is examined. In section III.B, the DIMET regime is studied using MD simulations in the adiabatic representation. The role of IR+UV/vis control scheme for different delay times is investigated. The two approaches, nonadiabatic quantum and adiabatic classical, will be compared in section III.C. II. Methods and Models A. Two-Mode Model: The Ground-State Potential. In the present work, we revisit the two-mode model used to study DIET in ref 21 to consider the more general process of femtosecond-laser-induced DIMET. The H-H distance, r, and the center-of-mass distance of the molecule perpendicular to the surface, Z, are included explicitly in this model. The 2D ground-state potential surface Vg(r, Z) is shown in Figure 1. Vg(r, Z) is a modified cut through the six-dimensional potential of Luppi et al.,22,23 i.e., Vg(r, Z; X0, Y0, θ0, φ0). Here, X and Y are the lateral, Cartesian center-of-mass coordinates of H2 relative to a Ru atom in the first layer (see ref 21), which are frozen to their equilibrium values X0 and Y0. Further, θ0 and φ0 are the frozen rotational degrees of freedom (polar angle and azimuthal). We choose X0 ) 0, Y0 ) (31/2)d/6, θ0 ) 90°, and φ0 ) 0°. Here, d ) 2.75 Å is the nearest neighbor distance on the Ru(0001) surface.21 The diffusion to neighboring fcc hollow sites and subsurface diffusion are blocked by additional barriers as we are interested in desorption from a (1 × 1) covered surface. The two-dimensional PES exhibits an adsorption minimum at Z0 ) 1.06 Å and r0 ) 2.75 Å with a binding energy

7792

J. Phys. Chem. C, Vol. 113, No. 18, 2009

Vazhappilly et al.

of 0.85 eV. This minimum corresponds to a situation where the two hydrogen atoms are located in adjacent fcc hollow sites. There is a small barrier of 0.18 eV for the incoming H2 molecule located in the entrance channel for dissociative adsorption and hence in the exit channel for associative desorption. This barrier (transition state) is located at Z‡ ) 2.24 Å and r‡ ) 0.77 Å. Apart from dynamics, Vg(r, Z) is also used to determine the bound vibrational states of the adsorbate, which are needed as initial states in the quantum model. For this purpose, the timeˆ gφn ,n (r, Z) ) independent vibrational Schro¨dinger equation, H r Z εnr,nZφnr,nZ(r, Z) is solved for 2H/Ru(0001) and 2D/Ru(0001) using the Fourier grid Hamiltonian (FGH) method.24 The lowest vibrational states can be classified according to their quantum numbers in the r and Z modes, i.e., nr and nZ. For convenience, we also denote vibrational states φnr,nZ as |nr, nZ〉 for short. The Hamiltonian for the two-mode model is given by 2 2 2 2 ˆ g ) - p ∂ - p ∂ + Vg(r, Z) H 2µr ∂r2 2µZ ∂Z2

is solved for the system density operator Fˆ s ) Fˆ gg|g〉〈g| + Fˆ ee|e〉〈e| + Fˆ eg|e〉〈g| + Fˆ ge|g〉〈e|. An uncoupled Hamiltonian

ˆs ) H ˆ g |g〉〈g| + H ˆ e |e〉〈e| H

(5)

ˆ i ) Tˆ + Vi (i ) g, e) is employed to describe the system, where H is the Hamiltonian of state |i〉. Since there is no coherent (direct) coupling between the two electronic states included in the model for the hot-electron mechanism (see below), the electronic coherences Fˆ eg and Fˆ ge of the reduced matrix remain zero and can be eliminated. In eq 4, the Cˆk are Lindblad operators and k denotes two different “dissipation” channels, i.e., excitation to and relaxation from the adiabatic representative excited state:

Cˆ1 ) √Γefg |g〉〈e|

(6a)

Cˆ2 ) √Γgfe |e〉〈g|

(6b)

(1)

Here, µr and µZ are vibrational masses along the respective coordinates, with µr ) mH/2 and µZ ) 2mH in the case of 2H/ Ru(0001) (mH ) hydrogen mass). B. Nonadiabatic Quantum Model. 1. Two-State Density Matrix Model of DIMET. In contrast to Langevin dynamics, where only the ground state Vg(r, Z) is needed, a two-state model with ground and excited states |g〉 and |e〉 is used in the nonadiabatic quantum model. The excited-state potential Ve(r, Z) is derived from results of cluster calculations (H2Ru3 cluster) at the B3LYP/LANL2DZ25,26 level of theory. Based on these cluster calculations,21 the excited-state potential is taken as a shifted ground-state potential

Cˆ1 accounts for the relaxation from the electronic excited state to the electronic ground state, with a rate Γefg ) 1/τel for the downward process. The excited-state lifetime is chosen as τel ) 2 fs in accordance with ref 21. Cˆ2 describes the hot-electronmediated excitation from ground to excited state, where Γgfe is the rate for the upward process. The upward rate Γgfe is timedependent and determined by the downward rate and a timedependent electronic temperature, Tel(t), by the principle of detailed balance28

Ve(r, Z) ) Vg(r - ∆r, Z + ∆Z) + Eex

Here, ∆V is the energy difference between Ve and Vg, which is coordinate-dependent in general. Since the lifetime is very short, however, the coordinate dependence of ∆V can be neglected, without altering the desorption considerably.13 Thus, we employ a constant ∆V ) Eex in what follows. 2. Monte Carlo WaWe Packet Method. The numerical solution of eq 4 is difficult even for a two-mode model because of large memory requirements. Here we employ the Monte Carlo wave packet (MCWP) approach12 to solve the LvN equation indirectly. The main advantage of this method, which has already been applied for photodesorption,29-31 is that the memory required for the propagation scales only linearly with the number of basis functions. In this method, one has to run different quantum wave packet propagations under the influence ˆ s - ((ip)/2)∑kCˆ†k Cˆk. An ˆ s′ ) H of a nonhermitian Hamiltonian H auxiliary wave function at time t + ∆t is obtained from ψ′(t + ˆ s′ is not hermitian, there is ˆ s′∆t)ψ(t). Since H ∆t) ) exp(-(i/p)H a loss of norm during the time interval ∆t, given as δp ) 1 〈ψ′(t + ∆t)|ψ′(t + ∆t)〉 ) ∑kδpk. (δpk is the loss of norm due to dissipative channel k.) Depending on the value of a random variable ε ∈ [0, 1] which is drawn in every time step, the wave function ψ(t + ∆t) is obtained either by simply renormalizing the auxiliary wave function (if ε > δp) or by enforcing a discontinuous “jump” between ground and excited states |g〉 and |e〉 (if ε < δp). In the MCWP method as realized heresi.e., ˆ ssthe wave packet is either without direct coupling term in H in the |g〉 or in the |e〉 state. A jump to the other state occurs by defining the new wave function as ψ(t + ∆t) ) Cˆkψ(t) and renormalizing it, where channel k is chosen with a probability pk ) δpk/δp in general. Since Cˆk contains the rates, the jump

(2)

where ∆Z ) ∆r ): ∆ ) 0.2 Å according to ref 21. Eex ) 1.55 eV is chosen in accordance with the experimental wavelength of 800 nm.3 In the present work, we used a slightly different functional form in the exit channel to avoid complications in the analysis of wave packets with “late jumps” during the stochastic wave packet propagation (see section II.B.2). For this purpose, ∆r is varied as follows

∆r (Å) )

{

0.2; 0.2 - a(Z - 2.1); 0.2;

Z < 2.1 Å 2.1 Å e Z e 5.2 Å Z > 5.2 Å

(3)

with a ) 0.0645. By this choice we get the same ground-state and excited-state potentials at the line Z ) Zdes, where the wave packets are analyzed (see below). To treat DIMET within the nonadiabatic quantum model, the single “representative” excited state |e〉 with potential Ve(r, Z) and finite lifetime τel is assumed to promote the desorption. The hot-electron-mediated excitation and simultaneous relaxation is described by an open-system density matrix approach. In the latter, a Liouville-von Neumann (LvN) equation in Lindblad form27

∂Fˆ s i ˆ , Fˆ ] + ) - [H ∂t p s s

2

∑ (CˆkFˆ sCˆ†k - 21 [Cˆ†kCˆk, Fˆ s]+)

k)1

(4)

Γgfe(t) ) Γefge-∆V/(kBTel(t))

(7)

Laser Desorption of H2 (D2) from Ru(0001)


TABLE 1: Parameters of Wave Packet Propagation and Analysis Used in the MCWP Method for the DIMET Model time propagation time step for propagation total propagation time

0.0121 fs 2.5 ps grid

number of grid points along r grid spacing along r number of grid points along Z grid spacing along Z

128 0.149 a0 512 0.039 a0

analysis number of states χVr to analyze wave packet time energy analysis line along Z (Zdes) absorbing potential starts at (Zabs) value of imaginary potential at last Z point

20 10.57 a0 10.6 a0 -i 0.0027 Eh

frequency is driven by the upward and downward rates, Γgfe(Tel(t)) and Γefg, respectively. The procedure is repeated for N times. For each trajectory, n, the expectation value of an operator, Aˆ, is calculated as An(t) ) 〈ψn(t)|Aˆ|ψn(t)〉 and overall expectation values of this operator are obtained as the average over all trajectories

1 〈Aˆ〉(t) ) N

N

∑ An(t)

(8)

n

This method solves eq 4 exactly with less computational effort when the number of realizations needed is not too large. The method can also be seen as a quantum version of the model of ref 5 with jumping classical trajectories. For each realization we run a 2D wave packet propagation on a grid using the split operator method with the fast Fourier transform technique as done in ref 21 for the DIET case. The analysis of observables is carried out in a similar fashion. To prevent unbound parts of the wave packet from reaching the grid boundaries, we use a negative imaginary absorbing potential that linearly increases after a point Zabs toward the end of the grid.32 We analyze the outgoing wave packet in a state-resolved fashion Via the time-energy method33 at a line Z ) Zdes, by projecting on vibrational eigenstates χVr(r) of H2 (D2). The parameters for the wave packet propagations and the analysis are given in Table 1. The quantum trajectories are propagated to a final time of 2.5 ps, which is longer than in the DIET case. This is done to include the substantial drop of Tel(t) toward longer times. The computation of observables is done by averaging over all the stochastic trajectories using eq 8. As initial states for the stochastic propagation we usually use the vibrational ground state ψ(0) ) φ0,0. In a few cases, we took vibrationally excited states φnr,nZ instead so as to model the effects of vibrational preexcitation. Based on computed dipole functions dz(r, Z) (with z denoting the z component of the dipole vector) and optimized pulses, it has been shown in ref 21 that the Z mode, in particular, can easily be excited with z-polarized infrared lasers (see below). 3. Two-Temperature Model. The time-dependent electronic temperature, which determines the upward rate Γgfe for electronic excitation according to eq 7, is calculated with the so-called two-temperature model (TTM) developed by S. I. Anisimov et al.10 The TTM is applied using a Gaussian-shaped laser pulse consistent with experiment,3 i.e., with a fwhm ) 130 fs and a wavelength of 800 nm. The TTM gives the

Figure 2. Electronic temperatures calculated for Ru from the twotemperature model with different laser fluences. The dotted curve represents the phonon temperature, Tph, for a fluence of 60 J/m2. Gaussian pulse parameters: fwhm ) 130 fs, λ ) 800 nm. The employed parameters entering the TTM are the following: electron-phonon coupling constant, g ) 185 · 1016 W m-3 K-1; electron specific heat constant, γ ) 400 J m-3 K-2; thermal conductivity, κ0 ) 117 W m-1 K-1; Debye temperature, θ ) 600 K; density, F0 ) 12 370 kg m-3; optical penetration depth at λ ) 800 nm, ζ ) 15.6 nm.

TABLE 2: Parameters Used in the Functional Form of ηZZ(r, Z) and ηrr(r, Z) parameter

value

a1 ) a 2 b1 b2 c Z1 Z2 r1 ) r2 σZ σr

0.3 meV ps Å-2 0.6 meV ps Å-2 0.1 meV ps Å-2 4 Å-1 2.36 Å 2.25 Å 0.74 Å 0.2 Å 0.3 Å

electronic temperature Tel(t) and the phonon temperature Tph(t) of the metal surface in response to the absorbed laser fluence, F, under the assumption that the excited metal electrons thermalize rapidly. The TTM is realized here for laser fluences which were also used experimentally, ranging from 60 J/m2 up to 120 J/m2. The material constants required for the TTM for a ruthenium surface are taken from refs 34 and 35. The calculated electronic temperatures for all fluences considered in the quantum mechanical model are shown in Figure 2. The figure caption lists the TTM parameters. All calculations were performed assuming that the surface temperature is 100 K prior to the vis pulse excitation, in accordance with the experimental conditions.3 In Figure 2, one can see that Tel(t) reaches a maximum around 350 fs before the onset of equilibration with the lattice. The maximum electronic temperature also increases with increasing ∝ F, in fluence. It roughly follows a power law, Tmax el agreement with earlier findings.36,37 C. Adiabatic Classical Model. 1. Classical Equations of Motion: LangeWin Dynamics. The generalized Langevin approach has been successful in treating the dynamics for many adsorbate-surface systems.38,39 The approach can also be used for nonadiabatic dynamics at metal surfaces, including substratemediated photochemistry.9,40,41 In this case, the method is also referred to as molecular dynamics with electronic frictions. In this model, one assumes that the metal electrons respond nearly instantaneously to the slow adsorbate motion. As a consequence, the dynamics occurs on a single potential energy surface42 and nonadiabaticity can be indirectly included in the form of electronic frictions. At finite electron temperature, Tel, random forces drive the nuclear motion.

7794


Vazhappilly et al.

The two-mode model described in section II.A with one vibrational and one translational coordinate is also employed here. The two-dimensional Langevin equations of motion for the dynamics on the ground potential energy surface are written as

µrr¨ ) -

∂Vg(r, Z) - ηrr(r, Z)r˙ + Rr(t) + Frfield ∂r

µZZ¨ ) -

∂Vg(r, Z) - ηZZ(r, Z)Z˙ + RZ(t) + FZfield (9b) ∂Z

(9a)

The first term in the right-hand side of both equations accounts for the force arising from the ground-state potential Vg(r, Z). The second term is a friction term, containing diagonal friction elements, ηrr and ηZZ, respectively. Off-diagonal elements ηrZ and ηZr of the friction tensor have been neglected. The terms ηqq (q ) r, Z) account for electronic friction and thus for vibrational relaxation in the electronic ground state (see below). The third term is a random force along mode q which must be present at finite electronic temperatures to maintain detailed balance. Finally, the last term is a term by which the direct vibrational excitation by an external (infrared) field can be modeled. 2. The Friction Terms. The friction coefficients ηZZ(r, Z) and ηrr(r, Z) for H/Ru(0001) have been determined by Luntz and Persson using perturbation theory and periodic density functional theory in ref 9. Accordingly

ηqq ) 2πp

δV |ψβ〉 | 2δ(εR - εF)δ(εβ - εF) (10) ∑ |〈ψR| δq R,β

where ψR(β) and εR(β) are Kohn-Sham orbitals and energies for combined band and k-point index R(β). Further, εF is the Fermi energy and (δV)/(δq) the functional derivative of the exchangecorrelation functional with respect to nuclear coordinate q. We have used data points published in ref 9 to create approximate analytical representations for ηrr(r, Z) and ηZZ(r, Z) as follows:

ηZZ(r, Z) ) a1

1 + 1 + ec(Z-Z1)

[

b1 exp -

ηrr(r, Z) ) a2

(r - r1)2 2σZ2

-

(Z - Z1)2 2σZ2

1 + 1 + ec(Z-Z2)

[

b2 exp -

(r - r2)2 2σr2

-

(Z - Z2)2 2σr2

] ]

slightly different ground-state potential.) Figure 3 shows that the analytical representation of the friction terms reproduces the main features of the DFT data reasonably well: (i) around the minimum of the potential (S ) 0, corresponding to r ) r0 and Z ) Z0), both friction coefficients are similar, ηrr ∼ ηZZ ∼ 0.3 meV ps Å-2; (ii) both functions show a strong coordinate dependence; (iii) around the transition state (S ∼ 2.2 Å, corresponding to r ∼ r‡, Z ∼ Z‡), ηrr has a maximum with ηrr ∼ 1/3ηZZ; (iv) for large S where H2 is far away from the surface both friction terms drop to zero as they should. 3. Random Forces. The relation between the friction coefficients and the random forces is given by the well-known fluctuation-dissipation theorem.43 For the sake of simplicity, the random forces are considered as a Gaussian process with infinitesimally small correlation time. The random force which follows from the above assumptions is Gaussian white noise with the following properties:44

〈Rq(t)〉 ) 0

(12a)

〈Rq(t)Rq(t′)〉 ) 2kBTelηqqδ(t - t′)

(12b)

Here, Tel is the electronic temperature of the bath. In the present case, Tel is calculated from the TTM again for various fluences and therefore is time-dependent. In practice, we use the Box-Muller algorithm45 in accordance with the fluctuation-dissipation theorem to construct the forces44,46

[

]

The parameters used in the above expressions are given in Table 2. The above analytical forms are compared with the electronic friction coefficients calculated by Luntz et al. In Figure 3 we show the analytical functions in comparison to the DFT data along the minimum energy path S of the ground-state potential as indicated in Figure 1. (Note, however, that the minimum energy path of ref 9 is slightly different from ours due to a

(13)

where ∆t and ηqq are the time step and electronic friction coefficient along the respective coordinates, and the random numbers, e and f, are uniformly distributed on the interval [0, 1]. 4. Field Terms. For IR-driven vibrational excitation, the field terms in eqs 9a and 9b are needed. In the semiclassical dipole approximation and under the assumption of a z-polarized IR field, they are given as

Frfield )

∂dz(r, Z) IR Ez (t) ∂r

(14a)

FZfield )

∂dz(r, Z) IR Ez (t) ∂Z

(14b)

(11a)

(11b)

1

2kBTel(t)ηqq) 2 Rq(t) ) (-2ln e)1/2 cos 2πf ∆t

Here, dz(r, Z) is the z component of the dipole moment as mentioned in section III.B. A functional form for dz(r, Z)

dz(r, Z) ) aZ2[tanh(b(r - c)) + 1]

(15)

based on first-principles cluster calculations, was suggested in ref 21. The parameters are a ) -0.0473 ea0-1, b ) 0.7 a0-1, and c ) 1.195 a0. The IR excitation of adsorbate modes is enforced by so-called π-pulses of sin2 shape, similar to those employed in the quantum mechanical model of ref 21


EzIR(t) )

( )

πt 2πp sin2 cos(ωt) tfindif tfin


(16)

In eq 16, tfin is the pulse length (chosen as 500 fs), ω the IR laser frequency, and dif ) 〈φi|dz|φf〉. In a quantum mechanical two-state model, π-pulses with ω ) ωif induce a perfect population inversion from initial vibrational state |i〉 to final state |f〉. In the classical model, this resonance condition is not strictly needed but appears to be beneficial as well. In the following we use IR parameters which are optimized for the Z mode, which is then dominantly excited. For this purpose we choose pω ) pωZ ) 136 meV, and dif ) d(0,0);(0,1) ) 60.8 · 10-3 ea0 for 2H/Ru(0001).21 With a 500 fs pulse, then, the IR laser fluence is 164 mJ/cm2, which was employed below. The IR excitation of the r mode is harder due to the smaller transition dipole moment, dif ) d(0,0);(1,0) ) 6.9 · 10-3 ea0 and was not attempted here. The transition dipole for Z excitation is larger because the dipole functions changes significantly along Z, but not along r.21 5. Integration of Equations of Motion. The 2D Langevin equations are propagated using the Ermak and Buckholz algorithm47 for many stochastic trajectories. All the trajectories are initialized on the ground-state surface at the adsorption minimum (r0, Z0) of the ground-state potential. In the absence of IR excitation, the initial velocities are chosen from a Boltzmann distribution such that the kinetic energy in each mode is kBTsurf, where kB is the Boltzmann constant. The surface temperature Tsurf ) 100 K is taken in accordance with experiment.3 Propagation for 1 ps on the ground-state surface was seen to be sufficient to achieve an equilibrated distribution in both the kinetic and potential energies of the subsystems variables, r and Z. After the initial thermalization period, the vis pulse comes in, heating the metal electrons and the lattice according to the TTM. As in the quantum mechanical model, the vis pulse was always chosen to have the form of a Gaussian temporal shape with 800 nm wavelength and 130 fs fwhm. The fluence was varied between 80 and 160 J/m2 in the classical model. In the case of the IR+vis strategy, both IR and vis pulses were applied. For simplicity, we also chose an initial surface temperature of Tsurf ) 0 K in this case. Although all (N) trajectories in the ensemble are initialized at the adsorption minimum, as before, the initial velocities are now set to zero so as to be consistent with the temperature, Tsurf ) 0 K. The results are converged with an integration time step of 0.01 fs which is used for all calculations. The desorption probability and the energy partitioning are analyzed at a point in the exit channel, Zdes, where the electronic friction is

negligibly small and the potential is flat. A trajectory crossing the line at Zdes is considered as desorbed. In the MD simulations, a larger value Zdes ) 9.5 Å was chosen than that in the quantum wave packet model because the additional computation time was negligible while allowing for better discrimination of the product channel. The expectation values for translational energy, Etr, and vibrational energy, Evib, are calculated by averaging over all (Nd) desorbing trajectories. The desorption probability is Nd/ N. III. Results and Discussion A. Nonadiabatic Quantum Model. 1. DIMET at a Fixed Fluence. We first consider the quantum mechanical, nonadiabatic MCWP model to describe the desorption of H2 or D2 from an initial ground vibrational state, φ0 ) φ0,0 at 100 K. The assumption that only the ground vibrational state is initially populated is appropriate for low initial surface temperatures. The effect of the variation of the laser fluence will be considered in the next section; in the present section, we use a typical value of 100 J/m2 for the fixed laser fluence. The convergence of the desorption probability is shown in Figure 4. The desorption probability is dominated by trajectories with multiple quantum jumpsswhich are relatively rare events (see below)sto a larger extent in the case of the deuterated isotope than it is for H2. This, in turn, requires the sampling of a larger number of quantum trajectories for D2 than that required for H2 to achieve convergence. Indeed, 8000 and 10000 quantum trajectories were required for D2 and H2 to obtain good accuracy at the fluence of 100 J/m2. The observed ratio of 5 in the isotope effect of the desorption yield is in reasonable agreement with experiment.3 Specifically, we find the desorption probabilities of H2 and D2 to be Y(0, 0) ) 3.78 × 10-3 and 7.59 × 10-4, respectively, at a fluence of 100 J/m2. Here, the notation Y(nr, nZ) is used to denote the desorption probability out of an initial state with quantum numbers |nr, nZ〉. The desorption yields computed here are larger than the desorption probabilities as found for the DIET case in ref 21. Therein, the same two-state model and the same lifetime τel ) 2 fs have been used, but only a single, sudden excitation was permitted. That treatment resulted in desorption probabilities of Y(0, 0) ) 8.07 × 10-4 for H2 and Y(0, 0) ) 8.81 × 10-5 for D2. Note that in the DIET model, this probability is in units per successful excitation event, whereas in the DIMET model the finite excitation probability is accounted for (see below). Thus, DIMET probabilities are desorption yields per pulse, which is why we will use both the notions “yield” and “desorption probability” in the case of DIMET below. In contrast, the

Figure 3. Electronic friction coefficients along the minimum energy path (S) obtained from the analytical form compared with the friction coefficients of Luntz et al.9 Left: ηZZ, Right: ηrr. S ) 0 corresponds to hydrogen adsorbed at Ru(0001) in the equilibrium position (r0, Z0), S ) 2.1 Å is the transition state, and S ) 3 Å corresponds to the situation where the molecule is about to leave the Ru surface.

7796


Figure 4. Desorption probability for H2 and D2 for a fluence of 100 J/m2 as a function of the number of quantum trajectories, N, using the MCWP method in the DIMET model.

computed DIET probabilities are typically much larger than DIET yields, because the excitation probability is typically much smaller than 1. Thus we find that DIMET is, by far, more efficient than DIET in agreement with earlier findings.5,4 The absolute desorption probabilities and other properties depend on the topology of the excited-state potential and on the excited-state lifetime, which are both not well-known. The assumed excited-state lifetime of 2 fs is in a reasonable range for adsorbates on metals and in line with experimental estimates. Larger lifetimes generally lead to larger probabilities but suppressed isotope effects in the desorption probabilities. For example, in the DIET case, the computed isotope effect in the desorption probability drops from a factor near 8 for τel ) 2 fs to about 3 for τel ) 10 fs.21 The excited-state potential is based on information from cluster calculations, with the same topology as the ground-state potential, but a shifted minimum. At least in the DIET case21 it turned out that most of the results do not strongly depend on the shift parameter ∆ if chosen from a reasonable range. Even a change of the sign of ∆ leads to no qualitative differences. Still, absolute desorption probabilities depend on excited-state potential and lifetime parameters. The energy partitioning of the desorbing molecules has also been analyzed within the DIMET model at a fluence of 100 J/m2. In H2, translations are “hotter”, Etr ) 454 meV, than the vibrations. In D2, a similar trend is seen, with a translational energy of Etr ) 403 meVand vibrational energy of Evib ) 146 meV. As in the experiment, we also obtain more energy transfer for H2 than for D2 in both coordinates. For example, the reported experimental values for D2 (at a fluence of 85 J/m2m) are Etr ) 430 meV and Evib ) 100 meV.2 The reason for the almost fourfold higher translational energy lies in the fact that the desorption dynamics is dominated by the topology of the ground-state potential. The excited state is only transiently populated, giving the wave packet merely a “kick” toward desorption. The ground-state potential has a “late” barrier along the desorption path. Thus the desorbing molecules should be translationally rather than vibrationally excited according to Polanyi’s rules.9,21 For DIET, similarly large differences between translation and vibration were observed in computational models.21 For H2, Etr ) 341 meV and Evib ) 43 meV, and Etr ) 311 meV and Evib ) 45 meV for D2.21 These figures show also that the internal energy distribution of the desorbing molecules depends quantitatively on whether femtosecond (DIMET) or nanosecond (DIET) lasers are used. 2. Dependence on Fluence. The desorption probabilities depend nonlinearly on the laser fluence in the DIMET regime. This peculiarity of DIMET has been investigated using different

Vazhappilly et al. fluences in the range F ∈ [60, 120] J/m2. At each fluence, we must run between 5000 and 20 000 trajectories to obtain the roughly 10% accuracy in the desorption probability comparable to that obtained in Figure 4. For D2 and low fluences, the yields are very small leading to substantially slower convergence and/ or lower accuracy. In Figure 5a, we show the desorption probabilities per pulse, with the excitation probability explicitly included, as computed with the MCWP approach as a function of F for H2 and D2. A power law fit Y ∝ Fn to the desorption probabilities results in n ) 4.9 for H2. For D2, a similar nonlinear behavior is found with an exponent n ) 6. The exponents are higher than the experimentally observed ones, i.e., n ) 2.8 for H2 and n ) 3.2 for D2,2 but in general agreement with experiment. The isotope effect, Ides, in the yields decreases with increasing laser fluence from ∼ 8 to 3 in the interval F ∈ [60, 120] J/m2. Similar trends were observed in experiments and previous MD simulations in which Ides shows a strong fluence dependence ranging from ∼10 to 5 in the same fluence range.2,9 The lowering in the ratio of Y(H2)/Y(D2) with increasing F results from the fact that the desorption probability increases more rapidly for D2 than for H 2. The nonlinear dependence of the desorption probability on fluence is best explained by the upward rate Γgfe. Compared to Tel(t), Γgfe increases dramatically with laser fluence as shown in Figure 6a. This results from detailed balance as stated in eq 7 according to which Γgfe is an exponential function of Tel(t). Therefore, an increase in the fluence greatly enhances the electronic excitation probability. As a result, more molecules reach higher vibrational levels on the ground electronic state after electronic relaxation, and this promotes desorption. In this context, it is important to enumerate the “multiple transitions” within DIMET and their effects on desorption. Within the MCWP method, the number of excitation/deexcitation cycles (double jumps) can easily be counted for each quantum trajectory. In Figure 6b, we show histograms giving the distribution of m-fold excitations, N(m), for different fluences F ∈ [60, 120] J/m2. Here, one can see that the most frequent excitation level depends on fluence showing a clear trend that for higher fluences the most probable excitation level m increases. We further calculate the aVerage number of excitations per pulse, 〈m〉 as

〈m〉(F) )

∑ m mN(m) N

(17)

where N ) ∑mN(m) is the total number of quantum trajectories. This average 〈m〉 also increases with increasing laser fluence. For example, at a fluence of 120 J/m2, the most probable excitation level is m ) 1, but higher excitations also occur and result in 〈m〉 ) 1.5. For F ) 60 J/m2, the adsorbate remains unexcited in most trajectories and 〈m〉 is only 0.13, with almost no multiple excitations taking place. The number of excitations determines the desorption yields per pulse. This is demonstrated in Figure 6c, where the desorption probability for H2 is given for cases when successively the multiple excitations are neglected by artificially removing the corresponding quantum trajectories from the averaging process. When taking all excitations for a given vis pulse into account (“all”), the yield increases as Y ∝ F4.6 as already seen in Figure 5. Once higher excitations are neglected, the exponent decreases down to the point in which the exponent is between 1 and 2 if only single excitations are kept. Thus, the



Figure 5. Dependence of the desorption yields for H2 and D2 on the applied laser fluence calculated with the quantum MCWP method, panel (a), and the classical Langevin case, panel (b). In the former, simulations are performed for F ∈ [60, 120] J/m2 and an excited-state lifetime of 2 fs. In the latter, fluences F are taken in the range [80, 160] J/m2. All curves are associated with a different value of the power law constant, n.

Figure 6. (a) Dependence of the upward rate Γgfe on the applied laser fluence calculated from eq 7. Parameters: Γefg ) (2 fs)-1, ∆V ) 1.55 eV. (b) The histograms give the number of realizations, N(m), for an m-fold excitation in the DIMET model using the MCWP method. (c) Desorption yields as a function of fluence if, artificially, higher excitations are successively neglected.

multiple excitations are responsible for the nonlinear increase of yield vs fluence. The DIET limit with an exponent 1 is not

quite reached in Figure 6c for the case when only single excitations are accounted for. This is also not surprising because of the high intensity of the femtosecond pulses and the presumption that metal electrons are thermalized, both conditions which are not met in standard DIET processes. The translational and vibrational energies of desorbing H2 and D2 molecules as a function of F are shown in Figure 7a. The computed translational and vibrational energies increase approximately linearly with fluence. The increasing energy content in the product channels with increasing fluence arises because the latter leads to more energy absorption and hence more energy which must be emitted with the desorbates. The translation remains hotter compared to vibration, and H2 possesses more energy than D2 in both modes. These two findings are in good agreement with experimental observations.3 For example, the translational energy of D2 grows from around 300 to 500 meV in the fluence range of the experiments,3 while our theoretical values are 338 and 443 meV (at 120 J/m2), respectively. In further agreement with experiment, we also notice that the ratio between translation and vibration, Etr/Evib, decreases from 5.6 to 2.3 for H2 and 4.3 to 2.1 for D2, respectively. This is due to the fact that the additional energy gained by the desorbing molecule at higher fluences is approximately equally distributed between translation and vibration; i.e., the difference between energies remains nearly constant while their ratio changes. 3. Effect of Vibrational Preexcitation. We also studied the effect of vibrational preexcitation on the DIMET probability for H2 using the MCWP method. Calculations were performed with the first vibrationally excited states in r, |1, 0〉, and in Z, |0, 1〉, as initial states, for a fluence of 100 J/m2. The desorption yield increases from Y(0, 0) ) 3.78 × 10-3, by about a factor of 2 when starting from the |1, 0〉 vibrational level, Y(1, 0) ) 7.24 × 10-3. Similarly, one obtains Y(0, 1) ) 7.8 × 10-3 for the Z-excited level. Thus, excitations of both the r and Z modes are about equally efficient. This can be explained by the fact that our excited-state potential is equally shifted along r and Z, by an amount ∆, with similar gradient along both directions. Thus, after the FranckCondon excitation, the wave packet moves simultaneously in both directions on the excited-state potential. These results are similar to those obtained in the previous DIET calculations, where, however, the relative enhancement factor was bigger.21 To conclude, vibrational preexcitation appears to be beneficial also for DIMET.

7798


Vazhappilly et al.

Figure 7. (a) Fluence dependence of vibrational and translational energies obtained from the MCWP simulations for different fluences for H2 and D2, in the fluence range F ∈ [60, 120] J/m2. (b) The same for classical Langevin case, and a fluence range of F ∈ [80, 160] J/m2.

Figure 8. Upper panel: Trajectory number of desorbing classical trajectories (for 20 000 realizations) plotted as a function of time. Lower panel: Electronic temperature as a function of time. The vis pulse starts after 1 ps equilibration period.

B. Adiabatic Classical Model. 1. DIMET at a Fixed Fluence. The associative desorption triggered by femtosecond laser excitations can also be described by the 2D Langevin model as outlined above. Following the same paradigmatic scenario as in section III.A.1, we assume that the adsorbate is unexcited initially apart from thermal motion at a surface temperature of 100 K. After the equilibration time of 1 ps, a 800 nm laser pulse with a duration of 130 fs (fwhm) and a fluence of 120 J/m2 is applied. According to the TTM the electronic temperature then increases, which gives random kicks to the adsorbate through the random forces Rr(t) and RZ(t). As a result, the adsorbate gains energy and may desorb on the ground-state potential. The statistics for the escape time of 20 000 realizations of desorbing trajectories are illustrated in Figure 8. The upper panel of the figure shows which trajectories desorb at which time. For example, a cross at trajectory number 5000 at time t ) 2 ps means that the 5000th trajectory has crossed the desorption line at Z ) Zdes after that time. From Figure 8, one can also conclude that a total propagation time of 3 ps including 1 ps of equilibration time is reasonable, to converge desorption probabilities. These equilibration and propagation times are used throughout all of the adiabatic classical model calculations below. Meanwhile, at this fluence, convergence of averaged properties required over 60 000 trajectories for H2 and 80 000 trajectories for D2 to obtain reasonable accuracy. At a laser fluence of 120 J/m2, we find a desorption probability of Y ) 2.55 × 10-2 for H2 and Y ) 5.46 × 10-3 for D2. The larger yield for H2 compared to D2 is again consistent with experimental observations.2 The absolute desorption probability

is larger than in the quantum dynamical results mentioned in section III.A, but the isotope effect is similar. Most of the desorptions take place after the peaking of the electronic temperature Tel as shown in Figure 8. This can be easily explained from the fact that the random forces are large if Tel is large. There is a time delay between the maximal electronic temperature and the highest desorption rate of about 0.5-1 ps. This suggests that it takes a certain time to climb up the vibrational ladder into the desorption continuum. We also computed averaged translational and vibrational energies of the desorbing molecules. The translational energies are 363 and 324 meV for H2 and D2, and the vibrational energies are 218 and 163 meV, respectively. This unequal energy partitioning is again in agreement with experimental results, the MCWP results, and previous calculations.9 Analysis of the desorbing trajectories gives no evidence for desorption of atomic hydrogen. 2. Dependence on Laser Fluence. The dependence of both the reaction probability and the energy partitioning within the desorbing species on the fluence are also characteristic hallmarks of DIMET in the classical Langevin model. The reaction yield behavior can be seen in Figure 5b through the fluence range F ∈ [80, 160] J/m2. Here, convergence of the desorption probabilities required 60 000-40 000 trajectories. A power law fit to the desorption probability can be made readily and resulted in exponents of n ) 3.4 and 4.3 for H2 and D2, respectively. These values are in close agreement to those obtained from experiments and previous MD calculations, where n ∼ 3 for H2.2,9 Figure 5b also suggests that the lighter isotope H2 desorbs more easily than D2 at all fluences. The isotope effect, Ides ) Y(H2)/Y(D2), decreases with increasing fluence from ∼8 to 3 in the range F ∈ [80, 160] J/m2. These observations are in agreement with experimental findings2 and the quantum DIMET results of section III.A. The energy partitioning and its dependence on fluence are shown in Figure 7b. H2 has more translational and vibrational energy compared to D2 as in the MCWP model. Translations are hotter than vibrations, as before. The energy gained by desorbates in both modes again increases approximately linearly with the laser fluence. For example, the translational energy of H2 changes from 336 to 393 meV in the fluence range F ∈ [80, 160] J/m2, while for D2 it increases from 300 to 360 meV. Similarly, the vibrational energy of H2 increases from 172 to 264 meV while for D2 it changes from 104 to 199 meV for the above fluence range. In the classical calculations, the ratio between translation and vibration, Etr/Evib, decreases as before.

Laser Desorption of H2 (D2) from Ru(0001) Specifically, the ratio ranges from 2 to 1.5 for H2 and from 2.9 to 1.8 for D2. 3. Effect of Vibrational Preexcitation. We also studied desorption of vibrational preexcited states using the IR+vis strategy. If not stated otherwise, a 500 fs π-pulse with a fluence of 164 mJ/cm2 was used for direct vibrational excitation of the Z mode as outlined above. For the vis pulse, which heats the metal electrons, a 130 fs pulse with 800 nm wavelength and a fluence of 120 J/m2 was adopted. In the MD with electronic frictions model, the vibrational relaxation of adsorbate modes is included although this was neglected in the quantum mechanical model (section III.A). In fact, electronic frictions and vibrational relaxation are now decisive. Without friction, there is also no hot-electronmediated excitation because of the fluctuation-dissipation theorem. In contrast, in the nonadiabatic DIMET models vibrational relaxation in the ground-state acts, if not unrealistically strong, as a small perturbation only.48 Before addressing the IR excitation, therefore, let us consider the time scales for vibrational relaxation for a moment. Quantum mechanically, the vibrational lifetimes are defined as the exponential decay time after which the first vibrationally excited state, |1〉, has relaxed to the ground state |0〉. After τvib, the population of |1〉 has dropped to 1/e of it is initial value. This lifetime can be estimated from the friction coefficients ηqq r (q ) r, Z) at the equilibrium position as τvib (r0, Z0) ∼ µr/(ηrr(r0, Z Z0)) and τvib(r0, Z0) ∼ µZ/(ηZZ(r0, Z0)), respectively. Using our parametrization of the friction coefficients, ηrr ∼ ηZZ ∼ 0.3 meV r Z ∼ 175 fs for the r mode, and τvib ∼ 700 fs ps Å-2, we find τvib for the Z mode. (These values remain very similar if the DFT data of ref 9 are used.) In ref 21, we have calculated the vibrational lifetimes from a cluster model instead (H2Ru3 cluster, B3LYP/LANL2DZ), directly from a Golden-rule expression r Z ∼ 110 fs and τvib ∼ similar to eq 10. In this case one finds τvib 600 fs. Finally, we can estimate the vibrational relaxation time from MD with electronic frictions. For this purpose, we considered single trajectories with an initial energy equal to pωq, with initial motion directed along coordinate q. We found damped oscillations in mode q, and an energy loss which is, in excellent approximation, exponential (not shown here). From there we determine Tvib ∼ 215 fs and Tvib ∼ 685 fs, respectively. Thus, all theoretical estimates suggest a subpicosecond vibrational relaxation time for both modes, with the r vibration relaxing by about a factor of 3-5 faster. The vibrational relaxation of D2 is slower by about a factor of 2 compared to H2, because D is heavier than H by this factor. Keeping this in mind, we now IR-excite the 2H/Ru(0001) system along the Z mode by 500 fs IR pulses. We find that a π-pulse pumps more vibrational energy in the system than expected for an ideal, quantum mechanical two-level system: The energy gain is more than 300 meV, while pωZ ) 136 meV. This is not a special classical effect. Rather, it is due to overtone excitation caused by the strong harmonicity of the ground-state potential, which is also observed quantum mechanically.21 The second, visible pulse can cause desorption. In computing averages, 50 000 trajectories were run for all calculations, resulting in an estimated accuracy of ∼10%. Different delay times between the IR and vis pulses are applied. The delay times, ∆τ, between the onsets of the IR and the vis pulses are varied between -5 ps and +4 ps. Negative delay times indicate that the IR pulse arrives prior to the vis pulse, while for positive ∆τ the vis pulse arrives first, as indicated in the upper panel of Figure 9. The lower panel of Figure 9 gives the desorption yield for H2 as a function of ∆τ. The slight fluctuations in the


Figure 9. Upper panel: Schematic representation of IR+vis strategy. “IR” shows the electric field of IR pulses for different delay times relative to the vis pulse, which leads to an electronic temperature profile as indicated. The lower panel gives the desorption yield of H2 for variable delay times, ∆τ ∈ [-5, +4] ps. The IR excitation of the Z mode is done by a π-pulse with a fluence of 164 mJ/cm2. The vis pulse excitation is done by a laser pulse with a fluence of 120 J/m2. The dashed curve shows the desorption probability in case of no IR preexcitation.

desorption probabilities for large positive or negative delays are within the error bars of statistical convergence. Figure 9 suggests that the IR pulse in combination with a vis pulse influences the desorption probability. When the IR and vis excitations are close to each other in time (∆τ ∼ 0), the yield reaches a maximum, which is larger than the probability in the “uncorrelated case” by a factor of 2, i.e., for large positive or negative ∆τ. In the “uncorrelated case”, the desorption probability is the same as without IR pulse excitationssee the dashed line in the figure. Looking closer, we also find that the largest desorption yield is obtained when the IR pulse is timed such that the maxima in the vibrational excitation and the electronic temperature are coincident. As a consequence, the yield vs ∆τ curve is asymmetric, and it is maximum shifted to slightly negative delays. For large and negative ∆τ, vibrational preexcitation is inefficient because vibrational relaxation took place before the vis pulse hits the surface. For large and positive ∆τ, vibrational excitation is also inefficient because the electron temperature Tel is already too low to cause large fluctuating forces. (Note, however, that the electrons are still “hot”. They are just no longer effective for desorption, as can also be seen from Figure 8.) Analogous calculations were carried out for D2, and for H2 but with other IR or vis fluences. In particular, lower IR fluences were considered as they may be more experimentally feasible. In all cases, which are not shown here, the results were qualitatively similar to those shown in Figure 9, with a similar enhancement factor that can, however, be lower than 2. To summarize, the MD with frictions approach suggests that vibrational excitation can play a beneficial role in the hotelectron-mediated desorption of adsorbates from a metal surface. This is a nontrivial finding. Substantial computer resources are required for the calculations using either the quantum MCWP or the classical MD model. Nevertheless, both find that even though the vibrational quantum is very small compared to the formal electronic excitation energy of 1.55 eV, it can play a role in controlling the desorption. C. Comparison of Nonadiabatic Quantum and Adiabatic Classical Models. From the above observations, it is evident that the MD results are in good agreement with the MCWP calculations for most properties. Isotope effects in the desorption probabilities, the energy content and the energy distribution in the desorbing molecules, the trends of increasing fluence, and finally the effects of vibrational preexcitation are all very similar

7800


and sometimes in almost quantitative agreement. The one exception is the total desorption probability, which is larger by a factor of between 5.5 and 3 for H2 in the MD model. Since the MD calculations are in slightly better agreement with experiment, this suggests that the quantum DIMET yields are too low. This may be due to the chosen excited-state potential and lifetime parameters, which are both not well-known, and which can influence the desorption probabilities. On the other hand, the MD with frictions model does not account for zeropoint effects, quantum tunneling, and nonclassical reflections. Experimentally, the absolute values of the desorption probabilities are not known. To test if tunneling and zero-point effects play a role, we have carried out nonadiabatic, classical calculations in analogy to ref 5, with classical trajectories jumping back and forth between states |g〉 and |e〉. The same two-dimensional potentials, the TTM, and the same transition rates as in the MCWP method were used. Preliminary results indicate that for H2/Ru the nonadiabatic classical and nonadiabatic quantum results are not dramatically different, even though hydrogen is involved. The question of which combination of approximations is the optimal one thus hinges less on the choice of “classical” vs “quantum”, but rather of “nonadiabatic” vs “adiabatic.” Of course, a nonadiabatic calculation is to be preferred if potentials and nonadiabatic couplings are reasonably well-known. The adiabatic model does not give much detailed insight into the excitation process and on the excited-state dynamics. In the MCWP method or related approaches, an excited wave packet will move initially not necessarily along the desorption pathway but, for example, toward the surface as in the celebrated Antoniewicz model.49 Such directional effects are much harder to describe by frictional dynamics. Further, if the excited-state lifetimes are long, which is the rule for large-gap substrates, the nonadiabatic picture is mandatory because the excited adsorbates remain in the excited state(s) for a longer time, if they relax at all. In contrast, for metallic surfaces ultrashort lifetimes are characteristic, and the dynamics proceeds dominantly in the electronic ground state. This is why the present MCWP and Langevin models yield results so similar: They use the same ground-state potential. In passing, we note that the classical Langevin approach seems computationally more attractive because it should be easy to run classical trajectories even on a multidimensional potential energy surface, which may be calculated accurately and efficiently. On the other hand, if f modes are to be considered, one has to know the f2 elements of the friction tensor, each element being a function of f coordinates in general. Here, further approximations will be needed in practice. IV. Summary and Conclusions In this paper, the role of vibrationally and electronically excited adsorbates for the molecular desorption of H2 (D2) from a Ru(0001) surface by femtosecond lasers has been investigated. We have employed two approaches, a quantum, nonadiabatic and a classical, adiabatic one. Within a two-mode model we were able to reproduce qualitatively and semiquantitatively the experimental observations for the associative desorption of H2 (D2). The system was first treated with a quantum mechanical, nonadiabatic model in the DIMET regime. Here, the adsorbate is electronically excited indirectly by finite electronic temperatures which were calculated from the two-temperature model. Multiple electronic excitations can occur which are correlated in DIMET. This helps to climb the vibrational ladder. The

Vazhappilly et al. nonlinear fluence dependence of the desorption yield is well reproduced. It is a consequence of the dramatic increase in electronic excitation rates which favors multiple excitations. A linear increase in the calculated translational and vibrational energies with fluence was found, with the former being larger due to the “late barrier” on the ground-state potential. The model predicts that desorption cross sections can be enhanced by vibrational preexcitation. We also applied an adiabatic, molecular dynamics model with electronic friction to treat DIMET. In particular, the previous model9 was extended to include IR excitation of vibrations. The friction model shares parameters with the nonadiabatic, quantum DIMET model such as the 2D ground-state potential and electronic temperatures. The coupling of the system to electron-hole pair excitations is included by fluctuating forces. The electronic friction coefficients which account for the vibrational damping of adsorbate modes were derived from semiphenomenological models with input from DFT calculations.9 This model was also successful in reproducing the experimental findings and previous theoretical calculations, especially the large isotope effect, nonlinear fluence dependence, and unequal energy partitioning. It is also computationally amenable to the calculation of the excitation of the Z mode for both H2 and D2 using shaped IR pulses. This, in turn, allowed us to investigate an IR+vis strategy. The time delay between both pulses was varied. The results also show that for this classical model, IR excitation prior to vis excitation promotes the reaction. To conclude, agreement with experimental findings was successfully achieved for both our nonadiabatic and adiabatic representations. Even though both models are based on different representations, we obtained qualitatively and semiquantitatively similar results. This proves that dynamics takes place predominantly in the electronic ground state. This work can be extended in many directions. One goal is to include more degrees of freedom, notably the rotational modes of the adsorbate. This can most easily be done within the classical Langevin approach, if ways are found to circumvent the “frictional bottleneck”. In the quantum dynamical simulations, the electronic excited-state and its lifetime should be improved by a higher level of theory. The performance of classical, nonadiabatic models (e.g., classical stochastic trajectories in a two-state system), and of quantum, adiabatic models should be studied in detail. Finally, the problem that electrons are not thermalized on ultrashort time scales after laser impact is not yet solved satisfactorily. Still, we posit that, in its present stage, the theory already has predictive power. Acknowledgment. We acknowledge support of this work by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich 450, project C7, and by the Fonds der Chemischen Industrie. R.H. acknowledges the Alexander-von Humboldt foundation and the National Science Foundation, CHE 0749580, for support. Fruitful discussions with C. Frischkorn (Berlin) and G. Fu¨chsel (Potsdam) are gratefully acknowledged. References and Notes (1) Denzler, D. N.; Frischkorn, C.; Hess, C.; Wolf, M.; Ertl, G. Phys. ReV. Lett. 2003, 91, 226102. (2) Denzler, D. N.; Frischkorn, C.; Wolf, M.; Ertl, G. J. Phys. Chem. B 2004, 108, 14503. (3) Wagner, S.; Frischkorn, C.; Wolf, M.; Rutkowski, M.; Zacharias, H.; Luntz, A. C. Phys. ReV. B 2005, 72, 205404. (4) Prybyla, J. A.; Heinz, T. F.; Misewich, J. A.; Loy, M. M. T.; Glownia, J. H. Phys. ReV. Lett. 1990, 64, 1537.

Laser Desorption of H2 (D2) from Ru(0001) (5) Misewich, J. A.; Heinz, T. F.; Newns, D. M. Phys. ReV. Lett. 1992, 68, 3737. (6) Saalfrank, P. Chem. ReV. 2006, 106, 4116. (7) Brandbyge, M.; Hedegård, P.; Heinz, T. F.; Misewich, J. A.; Newns, D. M. Phys. ReV. B 1995, 52, 6042. (8) Newns, D. M.; Heinz, T. F.; Misewich, J. A. Prog. Theor. Phys. Suppl. 1991, 106, 411. (9) Luntz, A. C.; Persson, M.; Wagner, S.; Frischkorn, C.; Wolf, M. J. Chem. Phys. 2006, 124, 244702. (10) Anisimov, S. I; Kapeliovich, B; Perel’man, T. SoV. Phys. JETP 1975, 39, 375. (11) Gadzuk, J. W.; Richter, L. J.; Buntin, S. A.; King, D. S.; Cavanagh, R. R. Surf. Sci. 1990, 235, 317. (12) Mølmer, K.; Castin, Y.; Dalibard, J. J. Opt. Soc. Am. B 1993, 10, 524. (13) Saalfrank, P.; Kosloff, R. J. Chem. Phys. 1996, 105, 2441. (14) Nakagami, K.; Ohtsuki, Y.; Fujimura, Y. Chem. Phys. Lett. 2002, 360, 91. (15) Manz, J.; Saalfrank, P.; Schmidt, B. J. Chem. Soc., Faraday Trans. 1997, 93, 957. (16) Boendgen, G.; Saalfrank, P. J. Phys. Chem. B 1998, 102, 8029. (17) Paramonov, G. K.; Saalfrank, P. J. Chem. Phys. 1999, 110, 6500. (18) Letokhov, V. S. Science 1973, 180, 451. (19) Shapiro, M.; Brumer, P. J. Chem. Phys. 1993, 98, 201. (20) Daniel, C.; Heitz, M. C.; Lehr, L.; Manz, J.; Schroeder, T. J. Phys. Chem. 1993, 97, 12485. (21) Vazhappilly, T.; Beyvers, S.; Klamroth, T.; Luppi, M.; Saalfrank, P. Chem. Phys. 2007, 338, 299. (22) Luppi, M.; Oslen, R. A.; Baerends, E. J. Phys. Chem. Chem. Phys. 2006, 8, 688. (23) Vincent, J. K.; Olsen, R. A.; Kroes, G.-J.; Luppi, M.; Baerends, E.-J. J. Chem. Phys. 2005, 122, 044701. (24) Marston, C. C.; Balint-Kurti, G. G. J. Chem. Phys. 1989, 91, 3571. (25) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 270. (26) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (27) Lindblad, G. Commun. Math. Phys. 1976, 48, 119.

J. Phys. Chem. C, Vol. 113, No. 18, 2009 7801 (28) Blum, K. Density Matrix Theory and Applications; Plenum Press: New York, 1996. (29) Saalfrank, P. Chem. Phys. 1996, 211, 265. (30) Saalfrank, P.; Boendgen, G.; Finger, K.; Pesce, L. Chem. Phys. 2000, 251, 51. (31) Gao, S.; Stro¨mquist, J.; Lundqvist, B. I. Phys. ReV. Lett. 2001, 86, 1805. (32) Neuhauser, D.; Baer, M.; Kouri, D. J. J. Chem. Phys. 1990, 93, 2499. (33) Balint-Kurti, G. G.; Dixon, R. N.; Marston, C. C. J. Chem. Soc. Faraday Trans. 1990, 86, 1741. (34) Bonn, M.; Denzler, D. N.; Funk, S.; Wolf, M.; Wellershoff, S.-S.; Hohlfeld, J. Phys. ReV. B 2000, 61, 1101. (35) Weaver, J. H. Physics Data; Fachinformationszentrum: Karlsruhe, 1981. (36) Corkum, P. B.; Brunel, F.; Sherman, N. K.; Srinivasan-Rao, T. Phys. ReV. Lett. 1988, 61, 2886. (37) Nest, M.; Saalfrank, P. Phys. ReV. B 2004, 69, 235405. (38) Tully, J. C. J. Chem. Phys. 1980, 73, 1975. (39) Adelman, S. A.; Doll, J. D. J. Chem. Phys. 1976, 64, 2375. (40) Tully, J. C. Dynamics of Molecular Collisions; Plenum, New York, 1976. (41) Head-Gordon, M.; Tully, J. C. J. Chem. Phys. 1995, 103, 10137. (42) Luntz, A. C.; Persson, M. J. Chem. Phys. 2005, 123, 074704. (43) Kubo, R. Rep. Prog. Phys. 1966, 29, 255. (44) Dohle, M.; Saalfrank, P.; Uzer, T. J. Chem. Phys. 1998, 108, 4226. (45) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, 2003. (46) Dohle, M.; Saalfrank, P.; Uzer, T. Surf. Sci. 1998, 409, 37. (47) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (48) Nest, M.; Saalfrank, P. J. Chem. Phys. 2002, 116, 7189. (49) Antoniewicz, P. R. Phys. ReV. B 1980, 21, 3811.

JP810709K

Femtosecond-Laser Desorption of H2 (D2) from Ru(0001): Quantum

Recommend Documents