Kinetics of Catalytic OH Dissociation on Metal Surfaces - The Journal

Feb 2, 2012 - A Tunnel Model for Activated Hydrogen Dissociation on Metal Surfaces. Ernst D. German , Hadas Abir , and Moshe Sheintuch. The Journal of...
2 downloads 0 Views 1013KB Size
Article pubs.acs.org/JPCC

Kinetics of Catalytic OH Dissociation on Metal Surfaces Ernst D. German* and Moshe Sheintuch Department of Chemical Engineering, Technion - Israel Institute of Technology, Haifa 32000, Israel ABSTRACT: A tunnel mechanism of OH dissociation reaction on metal surfaces is suggested and analyzed. Expressions for the rate constant and the activation energy are obtained by calculating the transition probability on semiempirical potential energy surfaces in the framework of the nonadiabatic perturbation theory. The analysis of the transition probability leads to a linear Brønsted−Evans− Polanyi (BEP) relationship between the activation energies and the surface reaction energies (ΔI), correlated by Ea = 0.73ΔI + 18 (kcal/ mol). This relationship is compared with published literature data. The effect of the electron coupling on the pre-exponential factor is discussed, and the kinetic isotope effect of hydrogen for the reaction is predicted.

1. INTRODUCTION Steam is one of the reactants in many commercial reactions1−3 like water−gas-shift (WGS) (CO + H2O → CO2 + H2) or steam reforming (CH4 + H2O → CO + 3H2) reactions. Steam activation includes the following steps of adsorption and stepwise dissociation of an isolated water molecule H2O + * → H2O* (1a) H2O* + * → OH* + H*

(1b)

OH* + * → O* + H*

(1c)

H2O* + O* → 2OH*

(1d)

energy barrier when proceeding either from a ground or from other excited vibration levels may be significant. Thus, accurate account of H tunneling in the OH dissociation reaction is the first objective of this paper. The tunnel mechanism has been employed in ref 6 for studying the kinetics of oxygen-assisted H2O dissociation on a catalytic surface. The potential energy of this system was considered as a function of quantum coordinate (O−H bond distance) and of classical coordinates characterizing vibrations of metal atoms near the reactant and products. In the present work, we introduce additional degrees of freedom (coordinates) that describe the distances of O and H atoms from a surface and the slope of the OH group relative to that in the initial state in which the OH line is perpendicular to the surface or is slightly declined to it. The second objective of this paper is to use the derived theory to calculate a Brønsted−Evans−Polanyi (BEP) relationship between activation energies and surface reaction energies and to test it versus published results. Supercell DFT calculations3 have shown such a BEP relationship. The correlations between kinetic and thermodynamic characteristics were demonstrated for many catalytic reactions13 using DFT supercell calculations for both the reaction energies and the transition states. The BEP relations are useful as they allow us to estimate the activation energy or rate constant for unstudied surfaces (for example, bimetallic surfaces) using only adsorption energies of reactants and products; calculation of the latter is significantly more simple than that of the transition states. Thus, a theoretical prediction of the numerical BEP coefficients without performing lengthy DFT calculations on the transition states seems to be important. Both problems are solved using the approach described below assuming the reaction to proceed by H tunneling, while a classical approach is applied for the barrier along other

The last three reactions are of similar type, i.e., catalytic A−H dissociation (where A is H, C, N, etc.); these may be described by a H transfer mechanism in analogy with reactions of proton and atom transfer in the gas and condensed phase.4,5 We have recently addressed steps 1b and 1d, showing the existence of strong tunneling effects.6,7 Here we focus on step 1c. Note that due to step 1d the system has an autocatalytic character as OH* dissociation produces O* (1c) which, in turn, results in more OH* (1d). The kinetics of step 1c has been examined in a number of supercell DFT simulations.2,3,8−11 In these works, the activation energy of the reaction was calculated using a classical description of all degrees of freedom, including the OH bond, along the reaction path. A quantum mechanical tunneling correction coefficient (based on the Wigner theory12) was included9 in the classical rate constant expression. The barrier to tunneling in this work was considered to be fixed and time-independent. This tunneling model is not justified for surface reactions in which the activation barrier is modulated by a potential field due to thermal vibrations of the metal atoms catalyzing the surface of the reaction. Furthermore, consideration of the tunneling in the framework of the Wigner theory assumes that the tunneling factor is small. However, the O−H frequency is very high (≫kBT/h), and the H atom tunneling under a potential © 2012 American Chemical Society

Received: November 6, 2011 Revised: January 30, 2012 Published: February 2, 2012 5700

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

The simplest form of the Ui and Uf is that obtained under the assumption that the vibration modes along the coordinates x, r, φ, and {qk} are not coupled. In this case, the PESs of the initial and final states may be written as a sum of the following contributions

degrees of freedom. The nonadiabatic perturbation theory is used in the present work to calculate the quantum transition probability on a semiempirical potential energy surface (PES) which is then averaged over all initial vibrational levels and summed up over all final ones. This paper is organized as follows. In the next section, a geometrical and potential energy model of reaction 1c is described. In Section 3, expressions for the rate constant, the activation energy, and the tunnel factor are derived. Numerical results in terms of the developed theory are presented in Section 4.

Ui(x , r , {qk }) = ui(r ) + vi(x) + Φ(φ) +

∑ wik(qk) + Ii k=1

(2a)

Uf (x , r , {qk }) = u f (r ) + vf1(x1) + vf2(x2)

2. MODEL The geometrical picture of adsorbed reactant and products follows published theoretical data,1−3,8−11,14−27 and it is shown in Figure 1. The OH bond is either perpendicular or tilted to

+

∑ wfk(qk) + If k=1

(2b)

The first terms in eqs 2a and 2b describe the potential energy of the O−H bond and the repulsive interaction of the nonbonded O and H atoms, respectively; the second term in eq 1 and both the second and third terms in eq 2b describe the interaction of the adsorbed OH fragment and adsorbed O and H atoms with the surface, respectively, in the perpendicular direction; the third term in eq 2a is the potential energy of the deformation vibrations of the H atom of the OH fragment due to the change in angle φ (Figure 1); the fourth terms in eqs 2a and 2b represent the potential energies of metal atom vibrations relative to their initial and final equilibrium values, qk0i and qk0f, which correspond the adsorbed OH fragment and adsorbed H and O atoms, respectively. The last terms in eqs 2a and 2b are the minimum energies of the corresponding diabatic PESs at the equilibrium values of all coordinates. Coordinates x and φ of the initial state are connected with coordinates x1 and x2 by equations (following from the conservation law for angular moments)

Figure 1. Pictorial view of an adsorbed OH fragment dissociation.

the surface by an angle π/2-φ depending on the adsorption site and on the metal: the tilted OH configurations are usually found for the ontop and bridge adsorption. The published data referenced above suggest that the preferred OH adsorption site is the fcc position on the (111) surfaces except for a platinum for which the preferred position is the bridged one; still, the energy difference between the bridge and ontop positions is rather small and seems to be comparable to computation accuracy (0.2 eV) of the DFT slab methods.28 The preferred adsorption site for both dissociated H and O atoms is the fcc one. The initial (Ui) and final (Uf) diabatic PESs that describe the transition from the initial stationary state (adsorbed OH fragment) to the final one (adsorbed O and H atoms) are introduced according to the reactants and products geometry portrayed above. The expressions for these PESs use wellaccepted functions with degrees of freedom that will be subject to optimization. In the suggested reaction model, the initial PES depends on the distance r between bonded O and H atoms, on the distance x of the center of mass of the OH group from a surface, on an incline angle π/2−φ and on a set of coordinates {qk} describing the vibrations of metal atoms near adsorption sites. The final PES depends on the distance r between dissociated atoms, on the distances of O and H atoms from the surface (x1 and x2, respectively), and on the coordinates {qk}.

x1 = x − k2r0i cos φ , k2 = mH /(mO + mH)

(3)

x2 = x + k1r0i cos φ , k1 = mO/(mO + mH)

(4)

cos φ = (x2 − x1)/r0i

(5)

where r0i is the equilibrium length of an O−H bond (∼1 Å). Now, we specify the terms in eqs 2a and 2b in explicit form. The potential energy ui(r) of the adsorbed OH fragment is approximated by a Morse function ui(r ) = Di(1 − exp[ −α i(r − r0i)])2

(6)

where Di and αi are the dissociation and anharmonicity parameters, respectively. The interaction between dissociated O and H atoms is assumed to have a repulsive character and described by an exponential function29 u f (r ) = Df exp[− 2α f (r − r0i)] + Δj

(7)

The mutual arrangement of the potential curves given by eqs 6 and 7 is shown in Figure 2; the Δj is identified with the transition energy along the coordinate r. The interaction energies of adsorbed particles OH, O, and H with the catalyst, vi(x), vf1(x1), and vf2(x2), are approximated by the Morse functions of the type vs = Bs(1 − exp[ −βs(x − x0s)])2

(8)

where subscript s is equal to i, f1, or f2 (see eqs 2a and 2b); numerical values of the parameters of these potentials are discussed in Section 4 below. Finally, the potential energy 5701

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

εf ,{μk}

Ψf

μ

εf

= φf (ξ; r , {qk })χf (r ; {qk }) ∏ ψ f k(qk ) k

where φγ(ξ; r, {qk}) are the electron wave functions (ξ are a set of the electron coordinates); χi,ft (r;{qk}) are the proton vibrational wave functions; and ψivk(qk) and ψfμk(qk) are the wave functions describing collective vibrations of lattice metal atoms in the initial and in the final state, respectively, and assumed to be harmonic ones; they are characterized by vibrational numbers νk and μk. Brackets in superscripts of eqs 11 and 12 denote sets of the vibration numbers of k- harmonic oscillators, and those in argument of the electron wave function denote sets of the normal lattice coordinates. 3.3. Transition Probability. In this subsection the average transition probability per unit time from the Ui to the Uf well (the rate constant) is calculated assuming the reaction to be nonadiabatic. The reason for this assumption is that sub-barrier transitions (tunneling) of the H atom are expected to be important for the considered reaction; see, for example, ref 9 in which the tunnel effect H atom for this reaction was shown to play a significant role. Due to this effect, the total electron−proton matrix element may be small; i.e., strong total nonadiabaticity is expected to be observed for this system even in the case of the large electronic coupling. The most appropriate method for describing nonadiabatic reactions is the Fermi Golden Rule.5 We calculate the average transition probability per unit time, P(x, φ), at a fixed distance x of the OH group from the surface and at fixed value of angle φ characterizing a slope of the OH to the surface

Figure 2. Potential energy curves along the r coordinate for the hydroxyl fragment, ui(r), and for dissociated hydrogen and oxygen atoms, uf(r), adsorbed on a metal surface (see eqs 6 and 7 in the text), uad(r) is the adiabatic curve; Δj is the transition energy along the r coordinate, and arrows show positive (up) and negative (down) energy values.

Φ(φ) along the deformation mode φ is approximated by a harmonic function Φ(φ) = (K /2)(φ − φ0i)2

(9)

where K is the force constant of this vibration and φ0i is the initial decline angle of the OH perpendicular to the surface; it may be equal to zero, for example, for the fcc adsorption.

3. RATE CONSTANT, ACTIVATION ENERGY, AND H TUNNELING 3.1. Physical Mechanism of OH Dissociation. In the initial state, the proton of an OH group is localized in the potential well ui(r) shown in Figure 2, and a set of discrete quantum vibration levels {εin} describes the states of the proton in this well. The proton occupies the vibrational levels in the well according to the Gibbs distribution. In the final state, the hydrogen atom may have any energy {εf} in the energy continuum spectrum of the decay potential uf(r). Therefore, we may describe the OH dissociation considering quantum transitions of a H atom from the nth vibrational level in the ui potential well to any level εf in the potential uf with subsequent averaging over the initial energy distribution and integration over the continuum of the final levels. For the quantum transition of the proton from the nth vibrational state in the OH to occur, it is necessary, according to the Franck−Condon principle,13,29 that the conservation of energy was fulfilled, i.e. ε f = εni − Δj ≡ εni − Δj(n , ε f )

2π P(x) = ℏZi

n ,{νk}

(13)

Ein,{vk}

where states, are

and

Efεf ,{vk},

E in ,{νk} = εni +



εf ,{μk}

Ef

= εf +

the total energies for the initial and final

ℏωki (νk + 1/2) + vi(x) + Φ(φ) + Ii

k=1

(14)

∑ ℏωkf (μk + 1/2) + vf1(x1) + vf2(x2) + If k=1

(15)

ωki,f

are the vibration frequencies of k-lattice oscillators; Zi is the total vibrational partition function for the initial state; V̂ is the perturbation; ρf is the density of states for free one-dimensional motion in the decay potential which may be written as30 ρf = 1/ℏΩ f

(16)

where Ωf is the vibration frequency in the normalization “box” and δ-function, which ensures the energy conservation law is equal to5a εf , μk

δ(E f =

ν

k

0

⎡ n ,{ν } ⎤ Ei k ⎥ ⎢ dε f ρ(ε f )exp − ⎢ kBT ⎥⎦ ⎣

εf ,{μk} ̂ n ,{νk} 2 εf ,{μk} n ,{ν } |V |Ψi ⟩| δ(E f − Ei k )

(10)

= φi(ξ; r , {qk })χ nf (r ; {qk }) ∏ ψ i k(qk )

∑ ∑∫ νk, μk n



× |⟨Ψf

We note that the initial and final energy states of the considered system include also vibrational energies of k classical metal lattice oscillators (phonon “bath”), so that the H-atom transfer is assisted by transitions from a set of these initial oscillatory levels to a set of the final oscillatory levels. The physical picture presented is the basis for the analytic expression for the rate constant (W) derived below. 3.2. Wave Functions. The double adiabatic approximation5,6 is used to determine the wave functions of reactant OH and products H and O, Ψin,{vk} and Ψfεf,{vk }; details may be found in ref 29. In this approximation, Ψi

(12)

(11)

− E in , νk )

1 2πikBT

i∞

∫−i ∞ dθ exp[−θ(Efεf , μk − Ein, νk)/kBT ] (17)

5702

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

It is assumed that ωki = ωkf = ωk. In this assumption, eq 13 may be transformed to the form5 P(x , φ , θ) = ∞

1 iℏkBT



since the integrand in eq 24 acquires a sharp maximum at x = x̂ and φ = φ̂ . Introducing notation H = vi(x) + Φ(φ) + F(x , φ , θ̂)

exp[− (εni − ε0i )/kBT ]

expanding function 25 at the saddle point in powers x and φ to quadratic terms, H(x,φ,θ̂) ≈ H(x̂,φ̂ ,θ̂) + (1/2)Hxx″(x − x̂)2 + (1/2)Hφφ″(φ − φ̂ )2 + Hxφ″(x − x̂)(φ − φ̂ ), diagonalizing the result to a quadratic form, and integrating the corresponding Gaussians we obtain the following equation for the probability

n=0

n , εf |2 /Z dε f ρf |Vep r

×

∫0

×

∫−i ∞ dθ exp[−{θ[ΔJ(n , εf ) + vf1(x1) + vf2(x2)

i∞

− vi(x) − Φ(φ)] + θ(1 − θ)Er}/kBT ]

(18)

W=

n,εf

where Vep is the electron−proton matrix element, which will be considered in detail in the next subsection and Zr is the proton vibrational partition function which is very close to unit and will be omitted in the future. Er is the reorganization energy of classical oscillators of the metal lattice which is equal to5 1 Er = 2

∑ mk ωk2(qk 0f

− qk 0i)2

(20)

where Δj(n,εf) is determined by eq 10 and ΔI is the total reaction energy. The integrand in eq 18 has a sharp maximum over θ and may be integrated approximately by the saddle point method.5,29 The integration yields 1 ℏkBT

P(x , φ) =

×



exp[−(εni − ε0i )/kBT ]



(21)

where

(27)

⎡ dv (x ) dv (x ) ⎤ dΦ =0 θ̂⎢ f1 1 + f2 2 ⎥ + (1 − θ̂) ⎣ dφ dφ ⎦ dφ

(28)

ωeff 2π

r

×

and the saddle point θ̂ is determined from the condition dF/dθ = 0 which gives

∑∫ z

n=0

z max min

4πγnep(z)ρf kBT

exp[− (Ĥ (z) + εni − ε0i )/kBT ] ZxZφ d1d2 /πkBT

dz (29)

where z = Δj/kBT is a new dimensionless coordinate; Ĥ (z) is defined by eq 32 below; and γepn denotes (referring to ref 5) the electron−proton transmission coefficient

θ̂ = [ΔJ(n , ε f ) + Er + vf1(x1) + vf2(x2) − vi(x) (23)

The quantity θ̂ has a geometric interpretation according to that it characterizes slopes of the initial and final PESs along the classical coordinates {qk} at the transitional configuration.5 The total transition probability W is obtained by the weighted integration of P(x, φ) over x and φ 1 ZxZφ

(26)

⎡ dv (x ) dv (x ) ⎤ dv (x) θ̂⎢ f1 1 + f2 2 ⎥ + (1 − θ̂) i =0 ⎣ dx dx ⎦ dx

W=

F(x , φ , θ̂) = θ̂[ΔJ(n , ε f ) + vf1(x1) + vf2(x2) − vi(x) − Φ(φ)] + θ̂(1 − θ̂)E (22)

W=



Equation 26 may be used for numerical calculations. For practical application, it is more convenient to transform it into another form. First, we note that the energy transition along the r coordinate is equal to Δj(n,εf) (Figure 2). According to eq 10, different εf values correspond to different Δj(n,εf). Therefore, for a given n the integration over εf is equivalent to the integration over Δj. Then, following ref 5, we introduce an effective frequency of fluctuation vibrations of classical subsystem (ωeff) multiplying and dividing eq 26 by (ωeff/2π) and rewrite W as follows

n=0 ̂ ∞ n , εf |2 exp[− F(x , φ , θ)/kBT ] dε f ρf |Vep 0 ″ θ̂ /2πkBT Fθθ

− Φ(φ)]/2Er

⎡ εn − ε 0 ⎤ ∞ n , εf i ⎥ exp⎢ − i dε f ρf |Vep ⎢⎣ kBT ⎥⎦ 0 n=0



In this equation, Jacobian |D|, which appears while the H(x, φ) is transformed to the diagonal form over x and φ, is equal to unity, and it is omitted. d1 and d2 are the eigenvalues of the diagonal form. The saddle point coordinates x̂ and φ̂ are determined from a set of eqs 27 and 28 in which x1 and x2 are the functions of x and φ (see eqs 3−5)

(mk is the mass of the kth oscillator). The quantity ΔJ(n,εf) in this equation is the transition energy along only the classical coordinates ΔJ(n , ε f ) = If − If − Δj(n , ε f ) ≡ ΔI − Δj(n , ε f )

1 ZxZφℏkBT

⎡ H(x ̂, φ̂ , θ̂) ⎤ 2πkBT πkBT ⎥ exp⎢ − × |2 ″ |θ̂ d1d2 kBT |Fθθ ⎦ ⎣

(19)

k

(25)

γnep(z) =

n , εf [Vep (z)]2

″ |θ̂ /π ℏωeff 2kBT |Fθθ

(30)

Equation 30 is applied if the γepn ≤1 (nonadiabatic transition); otherwise, it is taken to be equal to unity (adiabatic transition). To describe the smooth transition from the nonadiabatic to the adiabatic case, the quantity 4πγepn in eq 29 is replaced with the transmission coefficient Kn which is determined as5

∫ exp{−[vi(x) + Φ(φ)]/kBT }P(x , φ)dxdφ (24)

K n(z)̂ =

where Zx and Zφ are the configuration integrals. The probability W may be calculated approximately by the saddle point method 5703

1 − exp[−2πγnep] 1 − (1/2)exp[− 2πγnep]

(31)

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

One can see from eq 31 that Kn ≈ 4πγep at γep ≤ 1, and Kn = 1 if the opposite condition is met. The quantity Ĥ in eq 29 is the activation energy for a local n → εf transition which may be represented in a compact form using eqs 22 and 23

describing the tunneling H atom between the initial (εin) and the final (εf) levels within the resonance (see Figure 2); it depends, particularly, on the electronic coupling Ve between these states. The theory of the Vepn,εf for the electron nonadiabatic and adiabatic transitions was developed in ref 31 and according to the results which are employed below

Ĥ (z) ≡ H(x ̂, φ̂ , θ̂, z) = vi(x(̂ z)) + Φ(φ̂ (z)) + [θ̂(z)]2 Er E

(32)

n , εf Vep = Cn2 κ nℏ(Ωni Ω f )1/2 exp[−σn/2]

Here the first term is the repulsive contribution to the local activation energy due to the approaching OH center of mass to the surface to form the transition configuration (TS); the second one is the energy of its rotation to decrease the slope of the OH fragment to the surface at the TS; and the last term is the energy required to change configuration of metal atoms near the adsorption sites at the TS compared to the initial state. The value of Ĥ depends explicitly on reaction heat via the θ̂ (eq 23) and implicitly via the vi and Φ. These dependencies are considered in next section. Integration limits in eq 29 are determined by the mutual disposition of the curves ui and uf (Figure 2). The lower limit zmin=Δjmin/kBT is given by the crossing point of the curves ui and uf that coincides with the crossing point of the level εin and the curve ui(x). The upper integration limit zmax=Δjmax/kBT is given by such value of Δj at that the level εin is tangential to the curve uf. Preliminary numerical analysis shows that the contributions of excited terms (n > 0) in eq 29 are negligible as compared to that for n = 0 because the excited levels of the H atom in the initial potential well are sufficiently above the ground one. Therefore, one may approximately retain in eq 29 only the term corresponding to n = 0 and represent this equation in the identity form omitting subscript n = 0 W=



̂





̂

B

B

where, Ωin is the O−H frequency for nth vibrational level in the initial state; exp[−σn/2] is the tunnel factor; Ωf is a formal quantity for nonadiabatic reactions; and the product ρf(Vepn,εf)2 (see eq 16) really does not depend on Ωf. Coefficient Cn is equal to

Cn = π−1/4 2n/n! exp[(2n + 1)/4[ln(2n + 1) − 2ln2 − 1]

ρf kBT ωeff K (z ) 2π ZxZφ d1d2 /πkBT

Ĝ (z) = Ĥ (z) − kBT ln a(z)

p=

MVe2 ℏ2kcΔFc

(40)

where M is the reduced mass of the vibration in the ui potential, and kc and ΔFc are equal to kc =

2M(ui(rc) − εin) /ℏ

ΔFc = [dui(r )/dr − du f (r )/dr ]c

(41) (42)

Subscript “c” denotes the crossing point of the curves ui(r) and uf(r) (Figure 2), and Γ(z) is the gamma-function. In quasiclassical approximation and at strong electron coupling, the tunneling σn in eq 37 is calculated as

(33)

rr 2M [ dr[uad(r , {qk *}) − εni ]1/2 ] ℏ rl



σn/2 =

(34)

(43)

where rl and rr are the left and right turning points of the classical trajectories of energy εin on both sides of the barrier; the adiabatic potential curve uad is constructed from the diabatic potentials ui(r) and uf(r) by the usual way

(35)

and Ĥ (z) is given by eq 32. The dependence of the integrand in eq 33 on z is usually of the Gaussian form due to the ascending nature of K(z) vs the descending exp[−Ĥ (z)/kBT] shape when z varies from a lower to upper integration limit (the right curve in Figure 2 moves up relative to the left one). If the electron coupling is very strong and K(z) ∼ 1, the integrand is an exponential decreasing function. In the last case, integration in eq 33 is performed numerically; however, in the first case, integral 33 may be calculated approximately by the saddle point method W = a(z)̂

(39)

In eq 39, p is the adiabaticity parameter defined as

where a(z ) =

(38)

2πp exp[p ln p − p]/Γ(p + 1)

κn =



∫ a(z)exp⎢⎣− Hk (Tz) ⎥⎦dz ≡ ∫ exp⎢⎣− Gk (Tz) ⎥⎦dz

(37)

uad =

1 [(ui + u f ) − 2

(ui − u f )2 + 4Ve2 ]

(44)

If the electron coupling is small and the parameter p ≪ 1 (electron nonadiabatic transition), the tunneling is calculated under the diabatic potential curves ui(r) and uf(r) from the point rl to the crossing point rc and from the point rc to the point rr

⎡ Ĥ (z)̂ ⎤ 2πkBT exp⎢ − ⎥  A exp[−Ea /kBT ] ″ |z ̂ |Gzz ⎣ kBT ⎦

rc 2M [ dr[ui(r , {qk *}) − εni ]1/2 ℏ rl



σn/2 =

(36)

+

where |Gzz″|ẑ is the second derivative of G(z) at the saddle point ẑ and Ea Ĥ (ẑ). We can see from eqs 36 and 32 that the effective activation energy is determined as the value of the Ĥ (ẑ) at the optimal positions of the curves ui and uf. 3.4. Electron−Proton Matrix Element. In this subsection, we consider the electron−proton matrix element Vepn,εf

r

∫r r dr[uf (r , {qk*}) − εni ]1/2 ] c

(45)

Therefore, the factor κn in eq 37 describes both electron nonadiabatic and adiabatic reactions, so that eq 36 for the transition probability may be really used in a wide range of the electron couplings Ve. 5704

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

4. RESULTS AND DISCUSSION The methodology and origin of various parameter values used in the estimation process is described below. Then, the H atom tunneling and its effect on pre-exponential factor of the rate constant are discussed in Section 4.2. Finally, this information is used in Section 4.3 to derive the parameters of the BEP relation. 4.1. Input Parameters. In this subsection, the parameters of equations used for numerical calculations of the kinetic characteristics of the reaction are discussed. Some of them are given in Table 1. These parameters describe the adsorption

Figure 2) and therefore the pre-exponential factor via the shape of the potential barrier for H tunneling which is approximately considered to be constant for the OH fragment adsorbed on different metals. The Di value is assumed to be lower than that for gas phase OH dissociation (equal to 111 kcal/mol35). On the basis of results of ref 36 we believe that the decrease of this parameter for an adsorbed OH fragment as compared with the gas phase is proportional to a square root of the corresponding decrease of the stretching O−H frequency (3400 vs 3700 cm−1)35,37 that follows from the assumption that the form of the potential curve (i.e., anharmonicity parameter (αi)) does not change under OH adsorption. This leads to the Di = 105 kcal/mol. The equilibrium bond length OH is taken to be 0.99 Ǻ , and the deformation frequency of a H atom in an adsorbed fragment OH along the φ coordinate (see Figure 1) is taken to be 1000 cm−1.37 The values of both Df and αf for the dissociated potential uf (eq 7) are taken equal to the Di and αi, respectively, since the repulsive potential uf is assumed to be the same as the repulsive part of the potential ui (eq 8). The catalyst reorganization energy is taken to be equal to 6 kcal/mol as it is justified in our previous work;7 the effective frequency for the classical subsystem is ∼kBT/h, and the electron coupling is considered as a variable parameter. 4.2. H Atom Tunneling. A quantitative estimate of the preexponential factor is as important as that about the activation energy. However, in most works a quantum chemical DFT analysis of elementary steps is only reduced to calculation of the Ea assuming the pre-exponential factor to follow the classical ∼kBT/h value. This assumption is not obvious and may lead to overestimation of the pre-exponential factor for the H transfer reactions. We show below that for the considered reaction the pre-exponential factor is lower than kBT/h even at a rather strong electron coupling due to H tunneling, and only at the very strong Ve ∼ 1 eV the pre-exponential factor reaches the classical value. The calculations were performed for the OH dissociation on nickel and copper(111) surfaces at different values of the electron coupling; these systems were taken as model reactions characterized by small and high reaction energies, respectively (Table 1).

Table 1. Notation of Parameters and Their Numerical Valuesa A

OH

O

H

x0(A) ω⊥x (metal−A) nickel(111) −Eads ΔI copper(111) −Eads ΔI

1.5 440

1.25 450

1.15 1000

74.5

127.3

64.3

OHads → Oads + Hads

1.34 68.7

107.2

54.9 26.5

A is the adsorbed particle. Parameters of Morse functions: x0i  x0(OH), x01  x0(O), x02  x0(H); Bi  −Eads(OH), Bf1  − Eads(O), Bf1  −Eads(O); adsorption energies and ΔI in kcal/mol from ref 38; frequencies in cm−1 and equilibrium distances in Ǻ as explained in the text. a

properties of OH, H, and O species on specific metal surfaces (adsorption energies, adsorption heights, and metal species vibration frequencies in direction perpendicular to the surface); other parameters are the reorganization of metal atoms near adsorption sites and characteristics of the potential barrier for H tunneling. All these parameters determine the shape of the PESs, Ui and Uf. The parameters may be often taken from published experimental or computational works since the adsorption of the species has been investigated in many papers; otherwise, if we want to estimate the activation energy and the pre-exponential factor for the reaction proceeding on an unstudied metal, they should be calculated or measured. Note that the adsorption heights and the frequencies for each type of the adsorption sites vary in a rather narrow range of values for many metals and may be approximately considered to be constants in our calculations. The typical adsorption height of the fcc/hcp or bridge adsorbed OH on (111) metals is 1.5 Å,15,32 and those for the fcc/hcp adsorbed O and H atoms are approximately 1.2518,19,33 and 1.1 Å,15,21,34 respectively. The vibration frequencies of adsorbed OH, O, and H relatively metal surface in perpendicular direction may be approximately taken equal to 440, 450, and 1000 cm−1.1,11,14,21 While we mix here data from various sources (experimental, slab, or cluster), the good agreements between these sources and their limited effect on the W justifies our approach. The dissociation parameters Bi, B1f, and B2f of the potentials vi(x), v1f(x1), and v2f(x2) are identified with the corresponding OH, O, and H negative adsorption energies (see eq 8), and the anharmonicity parameters of these potentials are calculated using the frequencies in Table 1. The parameter Di is identified with the dissociation energy of an adsorbed OH fragment and does not affect directly the activation barrier of the reaction; it affects the transmission coefficient (see Subsection 3.4 and

Table 2. H Tunneling Characteristics for the Cusp Form and Adiabatic Potential Barriers (Figure 3) on Nickel for Different Values of Electron Coupling Vea Ve barrier type Δj;̂ eq 40 p̂; eq 40 κ̂; eq 40 σ̂/2; eq 43 or 45 γ̂ep; eq 30 K̂ ; eq 31 A; eq 36 Ea; eq 36 W; eq 33

5 cusp form 0.392

7.467 3.55 4.46 1.30 20.9 7.55

adiabatic

0.628 0.0284 0.377 7.216

× 10−7 × 10−6 × 108 × 10−8

5.87 7.38 2.10 20.7 1.63

× 10−7 × 10−6 × 108 × 10−7

10

15

adiabatic

adiabatic

1.819 0.115 0.625 6.344 9.21 1.16 3.40 19.8 1.20

× 10−6 × 10−4 × 109 × 10−5

4.045 0.262 0.769 5.298 1.13 1.42 7.93 19.1 1.63

× 10−4 × 10−3 × 1010 × 10−3

Ea, Ve, and Δj ̂ in kcal/mol; W and A in s−1; other characteristics are dimensionless; T = 300 K. a

Details of these calculations are listed in Tables 2 and 3. The values given in these tables are the kinetic characteristics for quantum transitions at specific values of the Δj ̂ at which the integrand in eq 33 has a maximum. The corresponding mutual 5705

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

Gaussian form (although it has a maximum), and it is not used; instead the activation energy and the pre-exponential factor are obtained from the rate constants calculated by numerical integration of eq 33 at two temperatures assuming the Arrhenius form of the rate constant temperature dependence. Therefore, the rate constant (W), the pre-exponential factor (A), and the activation energy (Ea) at this electron coupling listed in the last column of Tables 2 and 3 are the characteristics averaged over all quantum transitions; other data in the columns correspond to the specific transition at the Δj = Δj.̂ As it follows from these results, at the Ve = 15 kcal/mol, the reaction proceeds also by the nonadiabatic mechanism (A ∼ 1011 s−1 for both metals). For the electron coupling equal to 20 kcal/mol, the preexponential factor for the reaction on Ni reaches the classical value, and that for Cu is slightly lower; the integrand in eq 33 is an exponentially decreasing function. The result of calculation of the pre-exponential factors at this electron coupling is demon-

Table 3. H Tunneling Characteristics for the Cusp Form and Adiabatic Potential Barriers (Figure 3) on Copper for Different Values of Electron Coupling Vea Ve barrier type Δj;̂ eq 40 p̂; eq 40 κ̂; eq 40 σ̂/2; eq 43 or 45 γ̂ep; eq 30 K̂ ; eq 31 A; eq 36 Ea; eq 36 W; eq 33

5 cusp form

10 adiabatic

adiabatic

15 adiabatic

0.638

0.926

2.245

4.624

7.607

0.0284 0.377 7.390

0.117 0.626 6.602

0.263 0.769 5.699

2.69 × 10−7

4.16 × 10−7

5.50 × 10−6

5.08 × 10−5

3.38 × 10−6 7.94 × 107 39.6

5.23 × 10−6 1.20 × 108 39.4

6.91 × 10−5 1.69 × 109 38.3

6.38 × 10−4 2.70 × 1010 37.6

1.12 × 10−21

2.48 × 10−21

2.02 × 10−19

2.85 × 10−17

Ea, Ve, and Δj ̂ in kcal/mol; W and A in s−1; other characteristics are dimensionless; T = 300 K. a

locations of the potential curves (ui and uf) at these values of Δj ̂ are such that the uf curves are slightly higher than the r axis at the TS for both metals (since Δj ̂ > 0) (Figure 3). For the electron

Figure 4. Effect of the electron coupling on the effective preexponential factor.

strated in Figure 4 where a correlation between the pre-exponential factors and the electron coupling is presented. A similar correlation between the activation energies and the electron coupling is given

Figure 3. Nonadiabatic and adiabatic barriers for H tunneling from the ground level ε0i in the initial potential well to continuum levels of the final potential at different values of the electron coupling Ve (see eq 37); α = i, f, or ad; 1: Ve = 5 kcal/mol; 2: Ve = 10 kcal/mol; 3: Ve = 15 kcal/mol; 4: Ve = 20 kcal/mol; Di is the dissociation parameter of the initial potential ui.

coupling Ve = 5 and 10 kcal/mol, the kinetic characteristics may be considered as effective values of the dissociation reactions because in these cases the integrand in eq 33 is of the Gaussian form. At these values of the Ve, the electron adiabaticity parameter p̂ (eq 40) is very small; i.e., the electron is transferred on both nickel and copper surfaces by the nonadiabatic mechanism. The quantity σ̂/2 is considerable at this position of the potential curves which leads to the very small electron− proton transmission coefficient γ̂ep (or K̂ ) (given by eqs 30 and 31). As a result, the pre-exponential factor is much lower than its classical value (equal to ∼6 × 1012 s−1 at T = 300 K, Tables 2 and 3). We can see from these tables that only for the Ve = 5 kcal/mol the results obtained for the nonadiabatic cusp form, and the adiabatic forms of the potentials are close to each other. At the Ve = 15 kcal/mol, the barrier height for tunneling is rather small (Figure 3, curves 3 and 4). In this case, the approximation (eq 36) is poor since the integrand is not of the

Figure 5. Effect of the electron coupling on the effective activation energy.

in Figure 5. It is important that although at the Ve = 20 kcal/mol most of the quantum transitions are adiabatic there are a set of the nonadiabatic transitions near the upper integration limit which provide an important contribution to the total rate constant. 5706

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

Namely, due to these nonadiabatic transitions, we may expect the kinetic isotope effect (KIE) in the dissociation reaction. We calculated KIEs for the reaction on two considered metals at two temperatures, 300 and 500 K. Our estimates of KIEs are

which depends generally on ΔI may be represented according to eqs 23, 32, and 36 as

Table 4. Values of KIE Predicted for Hydrogen Dissociation on Cu and Ni(111) Surfaces at Two Temperatures and at the Electron Coupling Ve = 20 kcal/mol

(46)

dEa d[v (x)̂ + Φ(φ̂ )] = θ̂(ΔI ) + [1 − θ̂(ΔI )] i d(ΔI ) d(ΔI ) d[v (x ̂ ) + vf2(x2̂ )] + θ̂(ΔI ) f1 1 d(ΔI )

where all values are taken at the point ẑ which characterizes the maximum contribution of the integrand in eq 36. One can see from this equation that the linear relationship between the activation and reaction energies might be evident only if θ̂ is weakly dependent on the reaction energy and a sum of the second and third terms in eq 46 is close to zero. We found that slop θ̂ determined by eq 23 is really slightly varied in a wide range of the reaction energies from ∼0.753 for Ni to ∼0.792 for Cu at the Ve =10 kcal/mol due to mutual compensation of the components, and the result was found to be the same at Ve =5 kcal/mol. This finding may be considered as a justification of the linear BEP relationship (note that these θ̂ values are the effective ones which really correspond according to eq 23 to the reaction energies ΔJ = ΔI − Δj.)̂ The linear character of the Ea vs ΔI correlation allows a quantitative description of the BEP in the framework of the above theory using the activation energies only for any two metals; e.g., we chose nickel and copper. On the basis of these data and on the ΔI values of Table 1, we plotted in Figure 6 the correlation dependencies Ea vs ΔI for two values of the electron coupling (10 and 20 kcal/mol), by line (a), Ea = 0.723ΔI + 19.3, and by line (b), Ea = 0.735ΔI + 16.9 kcal/mol, respectively, where the ΔI value corresponds to reaction between nearest neighbors. Taking average values of the line parameters we obtain the approximate BEP relation in the form

T, K metal

ΔI, kcal/mol

300

500

nickel copper

1.34 26.5

15.5 30.9

8.3 12.3

given in Table 4. One should note that since OH dissociation is one of the steps of the total WGS process the effective value of the KIE which might be measured in the WGS reaction depends on the contribution of this step to the total rate constant. In any case, the experimental detection of the predicted KIE might be considered as a confirmation of our model. 4.3. Activation Energies and the Brønsted−Evans− Polanyi Relationship. While one objective is to estimate the rate, another, less ambitious objective, is to analyze trends within a group of metals. We are able to do it in an analytic form for the case of small or moderate electron coupling which allows introducing effective activation energy by eq 36. It follows from eqs 23, 32, and 36 that the activation energy would be of a quadratic (or in a special case a linear) function of the reaction energy ΔI (or ΔJ) if other components of the Ea are considered to be independent of the ΔI. However, these components depend really on the reaction energy which enters in eqs 23, 27, and 28 for determination of the saddle point coordinates (θ̂, x̂, φ̂ ) on the PES. Thus, the slop dEa/d(ΔI)

Ea = 0.73·ΔI + 18.1

(47)

where both Ea and ΔI are in kcal/mol. One can see that both lines in Figure 6 are compared favorably with published numerical results except the results of ref 38 for platinum which

Figure 6. Theoretical BEP relationships (in kcal/mol) between the activation energy (Ea) and the energy of the OH surface dissociation reaction (ΔI): Ve = 10 (a) and Ve = 20 kcal/mol (b); points in the figure correspond to slab DFT values calculated for the reaction on (111) surfaces of the following metals. 1: Pt, ref 10; 2: Ni, ref 3; 3: Ni, ref 8; 4: Pt, ref 3; 5: Pd, ref 3; 6: Pd, ref 9; 7: Pt, ref 1; 8: Rh, ref 8; 9: Cu, ref 3; 10: Au, ref 3; 11: Cu, ref 11; 12: Au, ref 2; 13: Ag, ref 38; 14: Cu, ref 38; 15: Pt, ref 38; 16: Pd, ref 38; 17: Ni, ref 38; 18: Ir, ref 38. 5707

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

transmission coefficient is investigated in a wide range of the electron couplings using the theory.31 It is shown that the H atom tunnels mainly from the ground vibrational level in OH to a small group of levels of energy continuum of the dissociated state. We found that the effective pre-exponential factor reaches the classical value at a strong electron coupling of ∼1 eV. However, even in this case part of the quantum transitions providing the contributions to the integral rate constant are of nonadiabatic character due to which a considerable hydrogen kinetic isotope effect (KIE) may be observed. Quantum chemical calculations of the electron coupling or measurements of KIE are needed to obtain accurate conclusions about the pre-exponential factor. Calculated linear BEP relation agrees well with published computed results using DFT. It is important that the developed approach may be used for prediction of the activation energies of the reaction on unstudied surfaces (for example, on bimetallic surfaces) using only adsorption energies of OH, O, and H species without performing lengthy DFT calculations on the transition states.

fall considerably from that given by eq 47. It seems that it is due to the reaction energy (ΔI) calculated in this work being very high (22.6 kcal/mol) which is greater than those calculated by other authors1,3,10 which are in the range from −1 to 6 kcal/mol. Although the trend described by the published data corresponds to that predicted by eq 47, most of the activation energies are higher than those given by eq 47. It results from ignoring of the OH extension contribution; the H atom tunnels along this coordinate. Sometimes, for practical application, it may be more convenient to use another form of the BEP relation based on the fact that the surface reaction energy may be also expressed through the OH, H, and O adsorption energies and the OH



AUTHOR INFORMATION

Corresponding Author

*Tel.: +972 4829 3561. E-mail: [email protected].

Figure 7. Thermodynamical cycle used for calculation of the energy ΔI of the surface reaction HO(ads) → H(ads) + O(ads): Eads(k) are the adsorption energies of the species k, and D(OH, gas) is OH dissociation energy in the gas phase.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Hadas Abir of our group for the data in Figure 6 (ref 38). This work is supported by the Clean Coal Technology Project and by the I-CORE Program of the Planning and Budgeting Committee and The Israel Science Foundation

dissociation energy in the gas phase (D(OH)). Applying the thermodynamic cycle (Figure 7) we find that ΔI(∞) = D(OH, gas) + Eads(O) + Eads(H) − Eads(OH) ≡ D(OH, gas) + ΔEads



(48)

where ΔI(∞) is the reaction energy at infinite separation of H and O. This energy is different from the energy at the reaction distance by a repulsion energy Q of the products (O and H), i.e., ΔI(∞) = ΔI + Q. Equation 48 allows representing the BEP relationship as the correlation of the activation energy vs the difference of the adsorption energies products and reactant Ea = s ·ΔEads + c

REFERENCES

(1) Grabow, L. C.; Gokhale, A. A.; Evans, S. T.; Dumesic, J. A.; Mavrikakis, M. J. Phys. Chem. C 2008, 112, 4608−4617. (2) Kandoi, S.; Gokhale, A. A.; Grabow, L. C.; Dumesic, J. A.; Mavrikakis, M. Catal. Lett. 2004, 93, 93−100. (3) Phatak, A. A.; Delgass, W. N.; Ribeiro, F. H.; Schneider, W. F. J. Phys. Chem. C 2009, 113, 7269−7276. (4) (a) Bell, R. P. The Proton in Chemistry; Chapman and Hall: London, 1973. (b) Bell, R. P. The Tunnel Effect in Chemistry; Chapman & Hall: London, 1980. (5) (a) Kuznetsov, A. M. Charge Transfer in Physics, Chemistry and Biology; Gordon and Breach: Reading, UK, 1995. (b) Stochastic and Dynamic View of Chemical Reaction Kinetics in Solution; Press Polytechniques et Universitaires: RomandesCH-1015 Lausanne, 1999. (c) Charge Transfer in Chemical Reaction Kinetics; Press Polytechniques et Universitaires: Romandes CH-1015 Lausanne, 1997. (6) German, E. D.; Sheintuch, M. J. Phys. Chem. C 2011, 115, 10063−10072. (7) German, E. D.; Sheintuch, M. J. Phys. Chem. C 2010, 114, 3089− 3097. (8) Pozzo, M; Carlini, G.; Alfe, D. J. Chem. Phys. 2007, 126, 164706. (9) Cao, Y.; Chen, Z.-X. Surf. Sci. 2006, 600, 4572−4583. (10) Michaelides, A.; Hu, P. J. Chem. Phys. 2000, 112, 6006−6014. (11) Gokhale, A. A.; Dumesic, J. A.; Mavrikakis, M. J. Am. Chem. Soc. 2008, 130, 1402−1414. (12) Wigner, E. Z. Phys. Chem. B. 1932, 19, 203−212. (13) (a) Liu, Z.-P.; Hu, P. J. Chem. Phys. 2001, 115, 4977−4980. (b) Michaelides, A.; Liu, Z.-P.; Zhang, C. I.; Alavi, A.; King, D. A.; Hu, P. J. Am. Chem. Soc. 2003, 125, 3704−3705. (c) Cheng, J.; Hu, P.; Ellis, P.; French, S.; Kelly, G.; Lok, C. M. J. Phys. Chem. C 2008, 112, 1308− 1311. (d) Logadottir, A.; Rod, T. H.; Nørskov, J. K.; Hammer, B.; Dahl, S.; Jacobsen, C. J. H. J. Catal. 2001, 197, 229−231. (e) Dahl, S.; Logadottir, A.; Jacobsen, C. J. H.; Nørskov, J. K. Appl. Catal. A: Gen.

(49)

where s = 0.73 is the slope and c = s([D(OH),gas] − Q) + 18.1; the parameter c should be estimated using an experimental value of the D(OH) if the adsorption energies are taken from experiments; otherwise, the D(OH) should be calculated in the same model as the adsorption energies.



CONCLUSIONS The tunnel mechanism of OH dissociation on metal surfaces is analyzed. This step is of importance in a variety of processes like steam reforming and water−gas-shift. Previous theoretical (DFT) studies of this step used a classical approach in determining the pre-exponential factor (PEF). In terms of our model, the initial and final semiempirical potential energy surfaces are introduced, and the activation energy and the PEF are calculated here using the nonadiabatic perturbation theory. The activation energy includes contributions due to the repulsion interaction between the OH fragment and metal surface, due to the rotation of the OH fragment around its center of mass resulting in the increase of an angle between the OH line and perpendicular to a surface and due to the reorganization of metal atoms adjacent to the reaction species. The electron−proton 5708

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709

The Journal of Physical Chemistry C

Article

2001, 222, 19−29. (f) Bligaard, T.; Nørskov, J. K.; Dahl, S.; Matthiesen, J.; Christensen, C. H.; Sehested, J. J. Catal. 2004, 224, 206−217. (g) Christensen, C. H.; Nørskov, J. K. J. Chem. Phys. 2008, 128, 182503. (14) Ford, D. C.; Xu, Y; Mavrikakis, M. Surf. Sci. 2005, 587, 159−174. (15) Mavrikakis, M.; Rempel, J.; Greeley, J.; Hansen, L. B.; Nørskov, J. K. J. Chem. Phys. 2002, 117, 6737−6744. (16) van Grootel, P. W.; Hensen, E.J. M.; van Santen, R. A. Surf. Sci. 2009, 603, 3275−3281. (17) Wang, G.-C.; Tao, S.-X.; Bu, X.-H. J. Catal. 2006, 244, 10−16. (18) Ganduglia-Pirovano, M. V.; Scheffler, M. Phys. Rev. B 1999, 59, 15533−15543. (19) (a) Xu, Y; Mavrikakis, M. Surf. Sci. 2003, 538, 219−232;(b) Surf. Sci. 2001, 494, 131−144. (20) Wang, G.-C.; Jiang, L.; Cai, Z.-S.; Pan, Y.-M.; Zhao, X.-Z.; Huang, W.; Kie, K.; Sun, Y.; Zhong, B. J. Phys. Chem. B 2003, 107, 557−562. (21) Yang, H.; Whitten, J. L. Surf. Sci. 1997, 370, 136−154. (22) Tang, Q.-L.; Chen, Z.-X.; He, X. Surf. Sci. 2009, 603, 2138− 2144. (23) Jiang, L.; Wang, G.-C.; Cai, Z.-S.; Pan, Y.-M.; Zhao, X.-Z. J. Mol. Struct. (Theochem) 2004, 710, 97−104. (24) Greeley, J.; Mavrikakis, M. J. Phys. Chem. B 2005, 109, 3460− 3471. (25) Li, T.; Bhatia, B.; Sholl, D. S. J. Chem. Phys. 2004, 121, 10241− 10249. (26) Yamagishi, S.; Jenkins, S. J.; King, D. A. Surf. Sci. 2003, 543, 12−18. (27) Xu, Y; Mavrikakis, M. J. Chem. Phys. 2002, 116, 10846−10853. (28) Greely, J; Norskov, J. K.; Mavrikakis, M. Annu. Rev. Phys. Chem. 2002, 53, 319−348. (29) German, E. D.; Kuznetsov, A. M.; Sheintuch, M. J. Phys. Chem. A 2005, 109, 3542−3549. (30) Merzbacher, E. Quantum Mechanics, 2nd ed.; Wiley: N.Y., 1970. (31) Georgievskii, Y; Stuchebrukhov, A. A. J. Chem. Phys. 2000, 113, 10438−10450. (32) Wilke, S.; Natoli, V.; Cohen, M. H. J. Chem. Phys. 2000, 112, 9986−9995. (33) Wang, Z.; Tian, F. H. J. Phys. Chem. B 2003, 107, 6153−6161. (34) Watson, G. E.; Wells, R. P. K; Willock, D. J.; Hutchings, G. J. J. Phys. Chem. B 2001, 105, 4889−4894. (35) Handbook of Chemistry and Physics, 81st ed.; CRC press: Boca Raton, Fl, 2000−2001. (36) Berfolini, J. C.; Tardy, B. Surf. Sci. 1981, 102, 131−150. (37) Bedürftig, K.; Völkening, S.; Wang, Y.; Wintterlin, J.; Jacobi, K.; Ertl, G. J. Chem. Phys. 1999, 111, 11147−11154. (38) Adsorption energies were calculated by DFT slab method with PW91 functional. H. Abir, personal communication (work toward PhD thesis).

5709

dx.doi.org/10.1021/jp2106499 | J. Phys. Chem. C 2012, 116, 5700−5709