The dissociative chemisorption of water on Ni(111): Mode- and bond-selective chemistry on metal surfaces Azar Farjamnia and Bret Jackson

Citation: J. Chem. Phys. 142, 234705 (2015); doi: 10.1063/1.4922625 View online: http://dx.doi.org/10.1063/1.4922625 View Table of Contents: http://aip.scitation.org/toc/jcp/142/23 Published by the American Institute of Physics

THE JOURNAL OF CHEMICAL PHYSICS 142, 234705 (2015)

The dissociative chemisorption of water on Ni(111): Mode- and bond-selective chemistry on metal surfaces Azar Farjamnia and Bret Jacksona) Department of Chemistry, University of Massachusetts, Amherst, Massachusetts 01003, USA

(Received 2 April 2015; accepted 5 June 2015; published online 18 June 2015) A fully quantum approach based on an expansion in vibrationally adiabatic eigenstates is used to explore the dissociative chemisorption of H2O, HOD, and D2O on Ni(111). For this late barrier system, excitation of both the bending and stretching modes significantly enhances dissociative sticking. The vibrational efficacies vary somewhat from mode-to-mode but are all relatively close to one, in contrast to methane dissociation, where the behavior is less statistical. Similar to methane dissociation, the motion of lattice atoms near the dissociating molecule can significantly modify the height of the barrier to dissociation, leading to a strong variation in dissociative sticking with substrate temperature. Given a rescaling of the barrier height, our results are in reasonable agreement with measurements of the dissociative sticking of D2O on Ni(111), for both laser-excited molecules with one or two quanta of excitation in the antisymmetric stretch and in the absence of laser excitation. Even without laser excitation, the beam contains vibrationally excited molecules populated at the experimental source temperature, and these make significant contributions to the sticking probability. At high collision energies, above the adiabatic barrier heights, our results correlate with these barrier heights and mode softening effects. At lower energies, dissociative sticking occurs primarily via vibrationally nonadiabatic pathways. We find a preference for O–H over O–D bond cleavage for ground state HOD molecules at all but the highest collision energies, and excitation of the O–H stretch gives close to 100% O–H selectivity at lower energies. Excitation of the O–D stretch gives a lower O–D cleavage selectivity, as the interaction with the surface leads to energy transfer from the O–D stretch into the O–H bond, when mode softening makes these vibrations nearly degenerate. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4922625]

I. INTRODUCTION

The dissociation of water on metal surfaces has been well studied because of the important role it plays in a broad array of chemical processes.1,2 For example, in the well-known steam reforming of natural gas, methane and water react over a Ni catalyst to form 3H2 + CO. This CO molecule can then react with a second water molecule to form CO2 + H2 (the watergas shift reaction). Both steps require water dissociation on the surface of a catalyst, and this can be a rate-controlling step of the water-gas shift reaction.3–6 Early theoretical studies of water dissociation on metals focused on electronic structure and kinetic modeling, and the barriers for water dissociation on a large number of transition metal surfaces have been computed via DFT (Density Functional Theory).1–10 Only very recently have the dynamics of the direct dissociative chemisorption of H2O come under close scrutiny. In this activated process, the molecule collides with the surface of the metal, breaking an O–H bond and leaving chemisorbed H and OH products. Tiwari and coworkers studied this reaction on Cu(111), treating the molecule as a pseudo-diatomic and fitting a model potential energy surface (PES) to DFT data.11 Their quantum approach treated a)Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2015/142(23)/234705/11/$30.00

3 + 1 molecular degrees of freedom (DOFs). In a later study, they examined how the dissociative chemisorption of water on Ni(100) and Ni(111) varied with the substrate temperature,12 using a sudden model to include the effects of lattice motion.13,14 Guo and co-workers have constructed a global PES for this reaction on Cu(111), fitting to tens of thousands of DFT energies.15–17 Considering impact only on the top site of a rigid flat surface, they treated 6 + 1 molecular DOFs, using quantum scattering methods to examine both the mode- and bond-selective chemistry of H2O and HOD.15–17 Only in the past year has the first molecular beam study of the dissociative chemisorption of water appeared.18 This study of D2O on Ni(111) was accompanied by a quantum scattering calculation by the Guo group, following a similar approach as in their Cu(111) studies, but using sudden models to average over impact sites19–21 and to include the effects of lattice motion.13,14 The theory was able to reproduce the observed large enhancement in reactivity upon excitation of the antisymmetric stretch.18 More recently, Jiang and Guo have developed a 9-DOF PES for this reaction on a rigid Ni(111) surface and used quasiclassical trajectories to explore the dynamics.22 We have developed a quantum approach for computing the dissociative sticking probabilities of molecules on metal surfaces, based on the reaction path Hamiltonian (RPH).23–25 This approach treats all of the internal motions

142, 234705-1

© 2015 AIP Publishing LLC

234705-2

A. Farjamnia and B. Jackson

of the molecule using a vibrationally adiabatic basis set. Various sudden models describe the other molecular DOF, as well as the effects of lattice motion.19,20,26 All input into the model comes from DFT. In studies of methane dissociation on Ni(100),19 Pt(110)-(1 × 2),27 and Ni(111),20,21,28 this approach has been able to reproduce and elucidate experimental trends in reactivity with respect to collision energy, surface temperature, and vibrational state, including modespecificity and bond-selectivity. In this paper, we use these new techniques to examine the dissociative chemisorption of H2O, HOD, and D2O on Ni(111), with several goals. First, we will examine and compare the variation in dissociative sticking with energy, vibrational state, and substrate temperature for these three isotopologues. Among other things, we are interested in how differences in the zero point energies (ZPEs) of the different isotopologues can modify activation barriers for reaction. We will also examine mode-specific behavior for these molecules and consider how vibrational excitation can lead to bond-selective cleavage in HOD. In gas-phase reactions of H with HOD, it was demonstrated that laser excitation of either the O–H or O–D stretch would selectively, and almost exclusively, break that bond, forming H2 + OD or HD + OH products, respectively.29,30 Guo and coworkers saw similar behavior in their theoretical study on Cu(111).17 We note that our vibrationally adiabatic formulation provides a useful framework from which to examine mode- and bond-selective chemistry, though we stress that it can also describe vibrationally sudden behavior, as we include all of the non-adiabatic couplings. We will also examine how surface temperature can modify this selectivity. Finally, this is the first application of our methods to a molecule other than methane, and while methane and water dissociations on metals share some similarities, this is an important test of our approach. Moreover, we can compare not only just with the experiment but also with the more standard approach of Guo and co-workers, whose PES is constructed using a similar level of DFT as ours. In Sec. II, we briefly review our RPH model. In Sec. III, we present results for the dissociative sticking for all three molecules, and we summarize our results in Sec. IV.

J. Chem. Phys. 142, 234705 (2015)

state, passing over the transition state (TS). The distance along 9  (dx i )2, the xi are the massthis path is s, where (ds)2 = i=1

weighted Cartesian coordinates of the water nuclei, and s = 0 at the TS. At 37 points along s, we compute the total energy, V0(s). We also compute and diagonalize the force-projected Hessian to find the 8 normal vibrational coordinates Qk and corresponding frequencies ωk(s) that describe displacements orthogonal to the reaction path at s, in the harmonic approximation. We write our PES, in the reaction path coordinates s and {Qk}, V = V0 (s) +

8  1 k=1

2

ω2k (s) Q2k .

(1)

In Fig. 1, we see that at the TS, the oxygen atom is roughly over a top site, 1.9 Å above the plane of the surface, with the dissociating H angled towards the surface, 118◦ from the surface normal. The barrier height is V0(0) = 0.71 eV, relative to the energy of the molecule and metal slab at infinite separation, and there is a physisorption well 0.23 eV deep. The reactive O–H bond length is 1.55 Å, nearly 0.6 Å larger than the gas phase value, and vibrational excitation is expected to enhance dissociation for this late barrier system. Our energies and geometries are consistent with earlier work.7–10,18 In Fig. 2, we plot the ~ωq(s) along the reaction path. When the molecule is far above the surface (large negative s), 3 of the frequencies are non-zero. These asymptotically bound modes are the symmetric stretch, νss, the antisymmetric stretch, νas, and the bending mode, νb. At large negative s, the remaining 5 modes (labeled 4–8) correspond to center of mass translation parallel to the surface (X and Y in Fig. 1), and rotation of the molecule. Closer to the metal, these asymptotically unbound modes describe hindered types of motion, and the frequencies become non-zero. The eigenvectors from the normal mode calculations, Li,k(s), define the transformation between the xi and our reaction path coordinates, x i = ai (s) +

14 

L i,k (s) Q k ,

(2)

k=1

II. REACTION PATH MODEL FOR DISSOCIATIVE CHEMISORPTION

All total energy calculations are performed using the DFTbased Vienna ab initio simulation package (VASP), developed at the Institut für Materialphysik of the Universität Wien.31–35 A 3-layer 3 × 3 supercell with periodic boundary conditions represents the metal as an infinite slab, with a water coverage of 1/9 mono-layer. Adding an additional layer or using a larger 4 × 4 supercell only lowers our barrier by 20 and 10 meV, respectively. A 15 Å vacuum space separates the slab from its repeated images. The interactions between the ionic cores and the electrons are described by fully nonlocal optimized projector augmented-wave (PAW) potentials,35,36 and exchange-correlation effects are treated using the PerdewBurke-Ernzerhof (PBE) functional.37,38 To construct our PES, we first locate the reaction path, or minimum energy path (MEP) from our reactant to our product

where the ai(s) give the geometry of the molecule on the reaction path at point s. Transforming to these new coordinates, our Hamiltonian has the following form:20,21,25,26
1 1 H = p2s + V0 (s) + Hvib − bss p2s + 2ps bss ps + p2s bss 2 4 1 − (ps πs + πs ps ) , (3) 2

FIG. 1. The transition state for water dissociation on Ni(111).

234705-3

A. Farjamnia and B. Jackson

J. Chem. Phys. 142, 234705 (2015)

)  1 2 ps + Veff ;k χ k + Fˆk χ0 + Gˆ k,q χ q i~ χ˙ k (s; t) = 2 q   √ + Fˆq χ q,k + Fˆq χ k,q + 2 Fˆk χ k,k , (9) (

q k

) 1 2 ps + Veff ;k,q χ k,q + Fˆk χ q + Fˆq χ k 2   Gˆ k, j χ j,q + Gˆ q, j χ j,k + (

i~ χ˙ k,q (s; t) =

j q

j >k

FIG. 2. Energies of the normal modes of H2O on Ni(111), along the reaction path for dissociation. The MEP, V0, and the zero point energy-corrected MEP, Veff,0, are also plotted, scaled by 1/2.



corresponding to the set of quantum numbers {nk}. These vibrationally adiabatic Φn are products of one-dimensional harmonic oscillator eigenfunctions that depend parametrically on s. Given our Hamiltonian, the coupled equations of motion for the wave packets, χn (s; t), are defined by the following:21,26,28 ( )  1 2 i~ χ˙ 0 (s; t) = ps + Veff ;0 χ0 (s; t) + Fˆk χ k (s; t) , (8) 2 k

√

+ 2



Gˆ q,k χ k,q ,

(11)

k 0 for displacements away from the bulk).43 We also compute how the location of the barrier (or the repulsive wall of the molecule-surface potential) shifts along the Z-axis as this metal atom vibrates. We have found that this is equal to αQ for reasonable values of Q. On the “flat” (100) and (111) surfaces of Ni and Pt, no other types of lattice motion significantly modify the barrier height or location. For methane dissociation on Ni(111), α = 0.70 and β = 1.16 eV/Å. For water dissociation on Ni(111), we find that several types of lattice motion can modify the height of the barrier to dissociation. The largest coupling, β = 0.61 eV/Å, is due to the normal motion of the Ni atom over which the water dissociates. The next largest coupling is much smaller, 0.10 eV/Å, and we find that α = 0.78. Our studies have shown that lattice motion is slow on collision time scales, and that a sudden treatment of Q is very accurate.13,14 We estimate the dissociation probability for a particular value of Q using our energy-shifting approximation, Eq. (15), where ∆V = −βQ. We then average P0 over many values of Q, each with the proper Boltzmann weighting at the substrate temperature T, using DFT to determine the energy required to distort the lattice by Q.19,20 After averaging over the type of lattice motion with the largest β coupling (Q), averaging over additional types of lattice motion (with smaller β) did not significantly change the results. We found similar behavior in our studies of CH4 dissociation on Pt(110)-(1 × 2), where several types of lattice motion can also modify the barrier height.27 In both cases, the additional couplings are typically weaker and tend to add out of phase. The α coupling is treated using an approach similar to the surface mass model.44 The relative water-Ni collision velocity depends upon α and the Ni atom momentum conjugate to Q, and we average over all values of this momentum at the substrate temperature T.14 We note that identical methods for averaging over X, Y, and lattice motion were used in the earlier quantum study of D2O dissociation on Ni(111).18 III. RESULTS AND DISCUSSION

When H2O is far above the metal, the two equivalent O–H stretches combine to make two non-local normal modes with similar energy: the symmetric and antisymmetric stretches. As the molecule approaches the TS, the reactive O–H bond, angled towards the metal, weakens and its vibrational frequency decreases significantly. The other O–H vibration is relatively unperturbed. This energy splitting decouples the vibrations, giving two normal modes that are each localized on one of the O–H bonds: the lower-frequency vibration on the reactive O–H bond correlates asymptotically with the lower frequency symmetric stretch, while the antisymmetric stretch correlates with the unreactive O–H stretch. We see this in Fig. 2 where the symmetric stretch significantly softens as we approach the TS. One consequence is that the activation energy for dissociation,

234705-5

A. Farjamnia and B. Jackson

J. Chem. Phys. 142, 234705 (2015)

TABLE I. Barrier heights along different adiabatic potentials for H2O, D2O, and HOD dissociation over Ni(111) for the ground state (gs) and the vibrationally excited states indicated. Vibrational state

H2O

D2O

Vibrational state

gs vas vss vb

0.536 0.514 0.182 0.428

0.590 0.572 0.335 0.510

gs vOH vOD vb

HOD, O–D cleavage

HOD, O–H cleavage

0.589 0.574 0.338 0.491

0.535 0.397 0.301 0.445

Ea, is lowered relative to V0(0). Ea is simply the ZPE-corrected barrier height equal to Veff,n(s = 0) − Veff,n(s = −∞), where n is the vibrationally adiabatic ground state. In Table I, we list Ea for the three isotopologues of water. For H2O, these ZPE corrections lower the barrier by about 0.17 eV. We find similar behavior for CH4 dissociation on Pt and Ni surfaces, where softening of the symmetric stretch lowers the activation energy by more than 0.1 eV relative to V0 (0).43 This softening also modifies the barrier heights, Veff,n(s = 0) − Veff,n(s = −∞), along the vibrationally excited adiabatic pathways.45 In Table I, we see that for molecules in the 1νss state, this is particularly dramatic. Thus, in the absence of any vibrationally non-adiabatic effects, we expect that excitation of the symmetric stretch should significantly enhance reactivity more than exciting any other mode. This is illustrated in Fig. 3 where we plot P0, the single-site rigid-surface reaction probability, for two levels of coupling. In the vibrationally adiabatic limit, corresponding to Bq,k = 0 and Bk,9 = 0, the reaction probability drops precipitously for energies below the vibrationally adiabatic barrier heights in Table I, as tunneling becomes the only mechanism for reaction. However, when we include all of the vibrationally non-adiabatic couplings, the reactivity at lower energies can increase dramatically. The curvature coupling, Bk,9, links states that differ by one quantum of vibration. Its primary effect is to lower reactivity near saturation, as energy flows to more excited states. It can also enhance the reactivity of vibrationally excited molecules at energies below the activation energy, by coupling to the ground state. This converts vibrational energy into motion along the reaction path. However, as in our CH4 studies,

FIG. 3. Single-site rigid-surface reaction probabilities for H2O on Ni(111) for molecules initially in the ground state (gs) or with one quantum of vibrational excitation in the mode shown. For the adiabatic case, all couplings Bq,9 and Bq,k are zero.

this effect is weak.28 The dominant pathways for reaction at low energies involve the Coriolis couplings Bq,k. Consider a molecule initially in the 1νss state. As it approaches the surface, the vibrational energy becomes localized on the reactive O–H bond. This localization is preserved as we transition between vibrationally adiabatic states at the avoided crossings, determined by the Bq,k. Transferring to states of lower vibrational energy also converts the excess vibrational energy into translational motion along the reaction path. Thus, starting in the 1νss state, a molecule can hop to the 1νb state, then to the 1ν4 state, and so on, preserving energy localization on the reactive bond and converting the vibrational energy into motion along the reaction path, which corresponds to O–H cleavage at the TS. This allows for over-the-barrier reaction down to Ei ≈ 0.1 eV. Molecules in the 1νb and 1νas states exhibit similar behavior. The enhancement due to the Coriolis couplings (relative to the adiabatic limit) is particularly large for the 1νas state; because there is no νas mode softening, the adiabatic reactivity is similar to that of the ground state. It is not obvious in Fig. 2, but Bas,ss is large in the entrance channel, mixing the two stretch states. Thus, molecules initially in the 1νas state can “hop” onto the 1νss state and follow a pathway like that described above. The reactivities of the two stretch states are thus very similar. We found similar behavior for CH4 on Ni(100) and Ni(111), for one of the three antisymmetric stretch modes. However, for CH4, the other two components of the antisymmetric stretch do not couple as strongly, leading to a lower overall efficacy for the antisymmetric stretch promoting reaction. In Fig. 4, we plot the dissociative sticking probabilities, S0, for the cases in Fig. 3, for a 300 K Ni substrate. To get S0 from P0, we average over the less reactive impact sites, which lowers the reaction probability. However, adding the effects of lattice motion increases it, particularly at lower Ei. We see that while the adiabatic approximation might work reasonably well for the softened vibrational modes, particularly at Ei ≥ Veff,n where the adiabatic reaction pathways are open, it fails completely for the antisymmetric stretch state, for the reasons discussed. All three vibrations significantly enhance dissociative sticking, and this is often described in terms of a

FIG. 4. Dissociative sticking probabilities of H2O on Ni(111), at 300 K, for molecules initially in the ground state (gs) or with one quantum of vibrational excitation in the mode shown. For the adiabatic case, all couplings Bq,9 and Bq,k are zero.

234705-6

A. Farjamnia and B. Jackson

J. Chem. Phys. 142, 234705 (2015)

vibrational efficacy, η=

∆Ei Ei (0, S0) − Ei (v, S0) = , ∆Ev ∆Ev

(16)

where ∆Ei is the increase in Ei necessary to give the same S0 as increasing the vibrational energy by ∆Ev. In Eq. (16), Ei(v, S0) is the incident energy giving a sticking probability of S0 for an initial vibrational state v, where v = 0 is the ground state. For S0 = 10−4, η = 0.95, 0.97, and 0.93 for the antisymmetric stretch, symmetric stretch, and bend modes, respectively. These results vary weakly with the choice of S0. Thus, unlike our results for CH4 dissociation on Ni(111) and Ni(100), we observe fairly statistical behavior; all efficacies are close to 1.0, and the stretch efficacies are very similar to the bend efficacies. We have found the symmetric stretches of CH4 on Ni(100)19 and Ni(111),20 and CHD3 on Ni(111)28 to be more efficacious for promoting dissociation than the antisymmetric stretches, though this results primarily from the degeneracy of the antisymmetric stretch, and the fact that only one of the degenerate modes couples strongly to the reaction path. This obviously does not play a role for the singly degenerate H2O stretches, which are mixed in the entrance channel by the Coriolis coupling. In their fully quantum fixedsite rigid-surface study of H2O dissociation on Cu(111), Guo and co-workers found similar behavior.15 For a dissociation probability of 10−4, they report η = 1.00, 1.08, and 1.17 for the antisymmetric stretch, symmetric stretch, and bend modes, respectively. While these efficacies are a bit larger, they are also close to statistical, perhaps surprising for this (very) late barrier system. In Fig. 5, we plot dissociative sticking probabilities for D2O on Ni(111) at two substrate temperatures (all full coupling). The ωq(s) for this √ reaction are very similar to those in Fig. 2, with a rough 1/ 2 scaling of the internal vibrations. As a result, S0 at 300 K is smaller for D2O than for H2O for all 4 initial states and all energies. For the ground state, the mode softening and ZPE effects are smaller, giving the increased Ea in Table I. For the excited states, the vibrationally adiabatic barriers are similarly larger, but more importantly, there is simply less energy in the excitations. At S0 = 10−4, η = 1.03,

FIG. 5. Dissociative sticking probabilities of D2O on Ni(111), at 100 and 300 K, for molecules initially in the ground state (gs) or with one quantum of vibrational excitation in the mode shown.

1.10, and 1.03 for the antisymmetric stretch, symmetric stretch, and bend modes, respectively, at 300 K. We again see roughly statistical behavior. Our efficacies are a bit larger for D2O than H2O, perhaps because the sticking probabilities are smaller and further from saturation. In addition, at higher Ei, above the corresponding vibrationally adiabatic barriers, the symmetric stretch state is more reactive than the antisymmetric stretch state, due to the mode softening (see Table I). One can also see this in Fig. 4 for H2O, but the effect is weaker. Decreasing the substrate temperature from 300 K to 100 K greatly decreases S0, particularly at low Ei where reaction is only possible at or near thermally puckered lattice sites. Note also that for fixed values of S0, the vibrational efficacies increase a bit as we go to lower temperatures, as puckered sites with lower barriers become less available. At S0 = 10−4 and 100 K, η = 1.14, 1.20, and 1.16 for the antisymmetric stretch, symmetric stretch, and bend modes, respectively. Our results in Fig. 5 for 300 K are in rough agreement with the calculations of Guo and co-workers for the same reaction and temperature.18 To construct their global PES, they use DFT with a similar functional (PW91), cell size, and slab thickness as ours. In addition, they use our methods for averaging over impact sites and for including the effects of lattice motion, so these parts of the two calculations are the same. Their reported values for S0 for the three excited states in Fig. 5 are similar to our values (see the supplementary material for Ref. 18). For example, for the antisymmetric stretch, the mode excited in the experiment, both methods give S0 ≈ 0.0003 at Ei = 0.3 eV, while at Ei = 0.75 eV, S0 ≈ 0.01 in their calculation and S0 ≈ 0.03 in ours. S0 for the symmetric and antisymmetric stretches are roughly equal to each other in both calculations. On the other hand, our ground state results are much larger than theirs, typically by over an order of magnitude. Both calculations compute values for S0 that are much larger than measured in the experiments of the Beck group.18 Guo and co-workers thus rescale their PES using the following function: ( ) Vscaled = V 1 + e−λZ , (17) where λ was set equal to 1.0 Å−1. This increases their barrier height by about 0.1 eV, significantly improving agreement with the experimental data. We take a similar approach here, but use Eq. (12) to rescale V0(s) only, where we compute Z as a function of s. We find that using λ = 0.6 Å−1 works better for our PES, corresponding to a barrier height increase of about 0.2 eV. We note that both the PBE and PW91 functionals give similar results37 and are known to overbind. This leads to barriers that are too small by 0.1 eV or so, based on studies of H2 dissociation on Cu(111)46 and CHD3 dissociation on Pt(111).42 In Fig. 6, we plot our results using the rescaled PES, along with the experimental data.18 While the experimental results for the 1νas and 2νas states are independent of the nozzle temperature, the “laser off” sticking probability contains contributions from vibrationally excited molecules and is larger than the true ground state result. The vibrational state distribution in the laser off beam is roughly Boltzmann and described by the nozzle temperature, which is large, up to 773 K at the highest Ei. We have computed S0 for all vibrational

234705-7

A. Farjamnia and B. Jackson

J. Chem. Phys. 142, 234705 (2015)

FIG. 6. Dissociative sticking probabilities of D2O on Ni(111), at 300 K, for molecules initially in the ground state (gs), or with one or two quanta of excitation in the antisymmetric stretch (νas) The lines are from theory, using the rescaled potential, and the circles are from experiment.18

states that contribute to the laser off reactivity, averaging them to get the laser off results in Fig. 6. We also plot the true ground state reactivity for comparison, and we see that the enhancement in S0 due to the vibrationally excited molecules is large. Guo and co-workers found much the same behavior for this reaction,18 and we have seen similar behavior in studies of CH4 dissociation on Ni(111).20 Overall, the agreement with experiment is reasonable with regard to the variation with Ei and vibrational state, given the scaling, which essentially shifts our results along the energy axis by about 0.2 eV. Our values for S0 tend to be too large at saturation, which was typical for our methane studies.26 As they could not measure the ground state reactivity directly, Beck and co-workers reported a vibrational efficacy of η = 1.1 for the 1νas to 2νas “excitation” and commented that the efficacy for the 1νas excitation is likely to be larger.18 At S0 = 10−4, we compute η = 1.00 and 0.62 for the 1νas state and for the 1νas to 2νas excitation, respectively, for our rescaled PES. Interestingly, while both our scaled and unscaled PESs roughly reproduce the shifts between the laser off, 1νas and 2νas experimental sticking curves, our efficacies are too small when compared with experiment at this value of S0. Much of this stems from our too-large saturation values for S0, which limits the efficacy of the second quantum of vibrational energy, and leads to slopes for our sticking curves that are different from experiment. On the other hand, Guo and co-workers18 report an efficacy for the antisymmetric stretch excitation of 1.65–1.78. This calculation, likely to be more accurate than ours, has a somewhat better overall agreement with the data (the saturation values are better), but the efficacy would appear to be a bit too large. Further study will be needed to sort this out, but overall, both calculations are in reasonable agreement with the experimental data, particularly considering the complexity of the system, and given some scaling of the barrier height. There are two possible transition states for the dissociative chemisorption of HOD (see Fig. 1): one where the O–H bond breaks and one where the O–D bond breaks. In Fig. 7, we plot the vibrationally adiabatic frequencies along the MEP for these two configurations. Because of the mass asymmetry, the two stretching vibrations, νOH and νOD, are almost completely localized on the O–H and O–D bonds of the gas

FIG. 7. Energies of the normal modes of HOD on Ni(111), along the reaction path for dissociation on Ni(111), for O–H cleavage and O–D cleavage.

phase molecule, respectively, making the mode- and bondselective behavior of HOD relatively simple. The dissociating bond weakens (softens) as the molecule approaches the metal. Thus, for the O–H cleavage configuration, νOH softens up to the νOH-νOD avoided crossing. At this avoided crossing, the modes exchange character. As a result, to the right of this avoided crossing, the blue curve, asymptotically labeled “νOD”, now corresponds to the weakening reactive O–H bond, and this continues to soften up to the TS. For the O–D cleavage configuration, the curves asymptotically labeled “νOH” and “νOD” do not cross, and only the “νOD” mode softens. The νOH mode is thus only weakly connected to the reaction pathway. As a result, the ground state activation energy for O–H cleavage is similar to that for H2O, and the activation energy for O–D cleavage is similar to that for D2O (see Table I). We thus expect that the dissociative sticking probability for HOD will be roughly the average of that for H2O and D2O. We see this behavior in Fig. 8, where we plot S0 for H2O, HOD, and D2O on Ni(111) at 300 K, using the scaled potential. For molecules in the ground state, the sticking curves for H2O and D2O are shifted relative to each other along the energy axis by about 0.06 eV, the difference in Ea. Computing S0 for HOD using only either the O–H or O–D cleavage configuration of Fig. 7, we find sticking curves that are almost identical to that for H2O or D2O, respectively. To get the total sticking probability for HOD, we take the average for these two equally probable configurations, and the resulting S0 is thus smaller than for H2O. Thus, the difference in reactivity of H2O, HOD, and D2O in the ground state can be explained entirely in terms of Ea, and the differences in Ea are due entirely to differences in the ZPE corrections to the barrier. Guo and co-workers found similar behavior comparing the reaction probabilities of H2O and HOD on a rigid flat Cu(111) surface,17 as have we,

234705-8

A. Farjamnia and B. Jackson

FIG. 8. Dissociative sticking probabilities of H2O, HOD, and D2O on Ni(111) at 300 K, for molecules initially in the ground state, using the scaled potential. For HOD, also shown are the sticking curves for the O–H and O–D cleavage configurations of Fig. 7.

comparing S0 for the isotopologues of methane dissociating on Ni(111).28 In Fig. 9, we compare the sticking probabilities of H2O, HOD, and D2O in more detail, for the scaled PES. We have seen that the ground state S0 is consistent with the vibrationally adiabatic barriers in Table I. However, as illustrated in Figs. 3 and 4, this is not necessarily so for the excited states, as non-adiabatic reaction pathways are important at lower Ei, below the vibrationally adiabatic barrier heights. The Coriolis coupling in the entrance channel mixes the symmetric and antisymmetric stretch states of H2O, leading to similar values of S0 for the two states. We observe the same behavior for D2O, though the values for S0 are smaller than for H2O, as there is less energy in the stretching excitations. There is a somewhat larger enhancement in S0 for the symmetric relative to the antisymmetric stretch excitation at higher energies, where the adiabatic reaction pathways are open and the mode softening (lower barrier) can play a role. The vibrational efficacies for H2O and D2O for the shifted PES are nearly identical to those given earlier for the unmodified PES, as the morphology of the transition state is not changed. The efficacies of the bending states are large for all three molecules,

FIG. 9. Dissociative sticking probabilities of H2O, HOD, and D2O on Ni(111) at 300 K for molecules initially in the ground state (gs; black) or with one quantum of vibrational energy in the antisymmetric stretch (red), symmetric stretch (blue), or bend (green). For HOD, the red and blue curves correspond to the O–H and O–D stretch excitations, respectively. The scaled potential is used.

J. Chem. Phys. 142, 234705 (2015)

and at S0 = 10−4, η = 0.91 for the HOD bend, similar to that for H2O. As expected, S0 for the localized O–H and O–D stretch states of HOD are about the same as the stretch excited states of H2O and D2O, respectively, with efficacies η = 1.01 and 0.97 at S0 = 10−4. In Fig. 10, we consider the sticking probabilities for the O–H and O–D cleavage configurations of HOD separately. The ground state S0 behaves as discussed, given the variation in ZPE and Ea, and we thus see a significant preference for O–H over O–D cleavage for ground state HOD, except at very high energies. The enhancement in O–H cleavage with excitation of the O–H stretch is substantial, though this is expected for a vibration localized on the reactive bond, given relatively minor IVR (intramolecular vibrational energy redistribution). S0 increases by several orders of magnitude relative to the vibrational ground state at lower Ei, as molecules follow the sudden (non-adiabatic) reaction pathways discussed earlier. Excitation of the O–D stretch also enhances O–D cleavage relative to that for ground state molecules. The magnitude of this effect is large and similar to that for the O–H stretch enhancement of O–H cleavage, keeping in mind the scaling of the excitation energy with mass. Interestingly, at low energies, we see an enhancement in O–D cleavage when the O–H stretch is excited. This surfaceinduced IVR results from a weak Coriolis coupling between the O–H and O–D stretch states in the entrance channel (for the O–D cleavage configuration). Moreover, when the O–D stretch is excited, there is a significant enhancement in O–H cleavage at all energies. Note (in Table I) that for the 1νOD state, mode softening lowers the barrier for both O–H and O–D cleavages by close to a quarter of an eV, while for the 1νOH state, it lowers the barrier only for O–H cleavage. At Ei below this barrier, simply staying in the 1νOD state after passing through the first avoided crossing leads to O–H cleavage. A better way to understand this is to simply examine the two configurations in Fig. 7. For the O–D cleavage configuration, the O–H bond is unperturbed and the frequency difference between the two stretch vibrations remains large. There is thus very little energy transfer between the two stretches. For the O–H cleavage

FIG. 10. Dissociative sticking probabilities of HOD on Ni(111) at 300 K, for the O–H and O–D cleavage configurations, for molecules initially in the ground state (gs; black) or with one quantum of vibrational excitation in the O–H stretch (red), O–D stretch (blue), or the bend (green). The scaled potential is used.

234705-9

A. Farjamnia and B. Jackson

configuration, the O–H bond weakens, and at s ≈ −1 amu1/2 Å, the two stretching modes become roughly degenerate, and energy transfer is much more probable. These behaviors may limit the amount of O–D bond-selectivity possible via laser excitation of HOD. This is seen more clearly in Fig. 11 where we plot the O–H cleavage percentage for the data in Fig. 10. As discussed, excitation of the O–H stretch leads to nearly 100% O–H cleavage at lower energies, while excitation of the O–D stretch enhances the amount of O–D cleavage, but not to the same degree. Our results are similar to those reported by Guo and co-workers in their quantum study of the dissociative chemisorption of HOD on Cu(111),17 with some exceptions. Guo and co-workers report that the % of O–H cleavage for ground state molecules varies from 85% to 43% as the incident energy increases from 1.0 to 1.3 eV. This is qualitatively consistent with what we observe. We stress that the energy scales in the two studies are different, as the barriers on the two metals differ. In addition, the Guo calculation is for impact on a single site of a rigid flat surface, and we have seen how impact site averaging and lattice motion effects can significantly modify the sticking behavior. However, both calculations predict a strong preference for O–H cleavage at low energies, approaching a statistical 50% at high energies. Similarly, they report that the % O–H cleavage decreases from 71% to 45% over this same energy range for molecules with one quantum of energy in the bending mode. As in our calculation, at lower Ei, the bending excitation increases O–H cleavage above the statistical 50% we observe at higher energies. It is not clear why Guo and co-workers find an O–H cleavage percentage below 50% at high energies. They also report that the O–H stretch excitation leads to 96% and 93% O–H cleavages, and the O–D stretch excitation gives 4% and 7% O–H cleavages for Ei = 1.0 and 1.3 eV, respectively. The trends with energy and excitation are consistent with our results, though their calculation gives more O–D selectivity with excitation of that bond than we observe. Finally, we note that the bond selective behavior of HOD is similar in many ways to that of partially deuterated methane on Ni(111) and Pt(111). In particular, bond selectivity is largest at lower energies where vibrational excitation can more significantly impact S0. At large Ei, there is a loss of selectivity, as the reactivity curves saturate and the bond closest

FIG. 11. Percentage of O–H cleavage for HOD on Ni(111) at 300 K, for the scaled potential, for molecules initially in the ground state (gs) or with one quantum of vibrational excitation in the mode indicated.

J. Chem. Phys. 142, 234705 (2015)

to the TS is likely to break. For methane, this behavior was observed in quasiclassical studies on Ni(111)47 and Pt(111),48 and more recently using our quantum reaction path approach on Ni(111).28 IV. SUMMARY AND CONCLUSIONS

We have examined both mode-specificity and bond-selectivity in the dissociative chemisorption of H2O, HOD, and D2O on Ni(111) using a quantum model based on the reaction path Hamiltonian. The effects of lattice motion are added using sudden models developed in earlier studies. The motion of lattice atoms near the dissociating molecule can significantly modify the barrier to dissociation for this reaction, similar to what we have found for the dissociative chemisorption of methane on metal surfaces. This leads to a strong variation in S0 with substrate temperature for both reactions. Our computed S0 are too large when compared with available experimental data, and we have rescaled our PES, increasing the rigid-surface barrier height by about 0.2 eV. With this modification, our results are in reasonable agreement with measurements for the dissociative chemisorption of D2O on Ni(111) at 300 K for molecules in the laser off, 1νas and 2νas states. For this late barrier system, excitation of both the bending and stretching modes can significantly enhance dissociative sticking. The vibrational efficacies, η, vary somewhat from mode-to-mode, but they are all relatively close to one. This is different from what we have observed for methane, where there is a larger variation in η between the different stretch normal modes, and the bends tend to be less efficacious than the stretches. Given this large enhancement with vibrational excitation for all modes, it is necessary to include vibrationally excited molecules to properly describe the “laser off” experimental data. Vibrational time scales and reaction/collision time scales are similar, and neither adiabatic nor sudden behavior is likely to dominate the overall dynamics. At energies above the vibrationally adiabatic barrier heights, the reaction probabilities are large. Vibrational modes that soften have lower adiabatic barriers and typically have larger S0 at these energies. At energies below the adiabatic barrier heights, over-the-barrier pathways for reaction exist via non-adiabatic transitions at the avoided crossings to states of lower vibrational frequency, converting vibrational energy into motion along the reaction path of the new state. As for our methane studies, the Coriolis couplings are far more important than the curvature couplings for facilitating reaction at lower energies. While the overall dynamics are not necessarily sudden, the small fraction of molecules that do react at energies below Veff,n(0) do so via vibrationally sudden pathways that preserve the localization of energy on the reacting bond. Our results are for the most part consistent with the quantum studies of water dissociation on Cu(111)15–17 and Ni(111)18 by the Guo group. Guo’s sudden vector projection (SVP) model49 is qualitatively consistent with much of the computed behavior, perhaps because at typical energies the reaction pathways are sudden. One exception is that the SVP model underestimates the efficacies of the bend excitations, at least when compared with our and Guo’s quantum results.50

234705-10

A. Farjamnia and B. Jackson

For HOD, the two stretch excitations are localized on the O–H and O–D bonds, and laser-induced bond-selectivity should occur trivially so long as IVR is minor. However, we do observe some surface-induced IVR in this study: excitation of the O–H stretch can increase O–D cleavage, and vise versa. For the O–H stretch excitation, the effect is weaker, occurring mostly at lower energies where the sticking probabilities are small. For the O–D stretch, on the other hand, enhancement is large at all energies. For the TS configuration where the O–H bond breaks, interaction with the surface weakens (softens) the O–H bond, and close to the TS, the two vibrations become roughly degenerate, leading to energy transfer between the O–H and O–D stretches. We predict a preference for O–H over O–D bond cleavage for ground state HOD molecules, except at very high collision energies, due simply to ZPE effects and differences in Ea. Given this and the IVR behavior, we find that excitation of the O–H stretch in HOD gives close to 100% O–H cleavage selectivity at lower energies, while excitation of the O–D stretch should give a lower O–D cleavage selectivity. Overall, our reaction path approach has worked well for modeling the dissociative chemisorption of methane on several metal surfaces.19–21,26,27 It has successfully reproduced the variation in sticking with collision energy and vibrational state, as well as mode selective and bond specific behaviors. In conjunction with our sudden models for treating lattice motion,13,14 we have fully explained the variation in reactivity with substrate temperature. Moreover, using our admittedly crude “energy shifting” approximations to average over impact site and molecular orientation, our results for methane sticking have been in semi-quantitative agreement with experiment. In this paper, when applied to water dissociation on Ni(111), our methods appear to work reasonably well given the complexity of the system, but the direct comparison with available experimental data is far less quantitative. We note that so far the only molecular beam data available are that plotted in Fig. 6. Thus, for example, while our lattice sudden model predicts a strong variation in S0 with T, using both our static lattice results and those of the Guo group,18 there are as yet no measurements to confirm this. We do reproduce the strong enhancement in S0 with both Ei and excitation of the antisymmetric stretch seen in the data of Fig. 6, but the agreement is only qualitative and requires a rescaling of the barrier by +0.2 eV. This is a bit larger than the +0.1 eV one might expect for a DFT-PBE based PES.42,46 One problem might be our treatment of X and Y and our average over surface impact sites. Our energyshifting approach (Eq. (15)) should be most accurate when the reaction is dominated by a single transition state, and where the barrier height increases harmonically as we move away from the TS along X and Y. We also assume that the form of the variation of P0 with Ei is similar at all surface impact sites that contribute to S0 at a given Ei and T. Simply put, the approach works well when most of the reactive trajectories remain relatively close to our single MEP. For methane, we have found that the barriers for dissociation become much larger as we move away from the top site, though we have seen evidence that dissociation near other sites might contribute at very high collision energies.21 In their recent paper, Jiang and Guo report a minimum barrier for dissociation very similar to ours in terms of energy and location.22 However, they show

J. Chem. Phys. 142, 234705 (2015)

that the barriers at many other sites on the surface unit cell are relatively small, much smaller than our Eq. (15) gives, and that these other impact sites contribute significantly to S0 at collision energies relevant to the experiment. Moreover, the morphology of the PES, and the corresponding form for P0(Ei), can be very different at these other sites. It is important to stress that our model is designed to describe, quantum mechanically, the reaction dynamics of polyatomic molecules and also to include the effects of lattice motion while maintaining maximum simplicity. Our goal is not to reproduce the experimental data but to explain the data using the simplest model possible in which the physics is laid bare. Our PES is based on dozens of total energy calculations, not tens of thousands, and the scattering calculation can be done quickly on a single processor. All of the data needed to average over impact sites and molecular orientation come from a single Hessian calculation at the TS, and the lattice coupling is described by two simple, but physically meaningful parameters. Overall, the approach works surprisingly well. However, it could be improved. We can replace our single MEP with several reaction paths describing dissociation on “representative” surface impact sites. It might be sufficient in some cases to simply replace ∆V (X,Y) with one fit to DFT computations of the barriers at several sites. It is important to note, however, that the biggest modification to P0 necessary to compute sticking probabilities that can be compared with experiment is the lattice averaging. Our experience is that this averaging often washes out the details of P0(Ei), and that the final form for S0(Ei, T) is in large part determined by the saturation value of P0(Ei) and the energy at which it saturates. ACKNOWLEDGMENTS

B. Jackson gratefully acknowledges support from the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, U.S. Department of Energy, under Grant No. DE-FG02-87ER13744. 1P.

A. Thiel and T. E. Madey, Surf. Sci. Rep. 7, 211 (1987). A. Henderson, Surf. Sci. Rep. 46, 1 (2002). 3C. V. Ovesen, P. Stoltze, J. K. Norskov, and C. T. Campbell, J. Catal. 134, 445 (1992). 4N. Schumacher, A. Boisen, S. Dahl, A. A. Gokhale, S. Kandoi, L. C. Grabow, J. A. Dumesic, M. Mavrikakis, and I. Chorkendorff, J. Catal. 229, 265 (2005). 5A. A. Gokhale, J. A. Dumesic, and M. Mavrikakis, J. Am. Chem. Soc. 130, 1402 (2008). 6C. Ratnasamy and J. P. Wagner, Catal. Rev. 51, 325 (2009). 7G. C. Wang, S. X. Tao, and X. H. Bu, J. Catal. 244, 10 (2006). 8M. Pozzo, G. Carlini, R. Rosei, and D. Alfe, J. Chem. Phys. 126, 164706 (2007). 9A. A. Phatak, W. N. Delgass, F. H. Ribeiro, and W. F. Schneider, J. Phys. Chem. C 113, 7269 (2009). 10S. C. Huang, C. H. Lin, and J. H. Wang, J. Phys. Chem. C 114, 9826 (2010). 11A. Mondal, H. Seenivasan, and A. K. Tiwari, J. Chem. Phys. 137, 094708 (2012). 12H. Seenivasan and A. K. Tiwari, J. Chem. Phys. 139, 174707 (2013). 13A. K. Tiwari, S. Nave, and B. Jackson, Phys. Rev. Lett. 103, 253201 (2009). 14A. K. Tiwari, S. Nave, and B. Jackson, J. Chem. Phys. 132, 134702 (2010). 15B. Jiang, X. F. Ren, D. Q. Xie, and H. Guo, Proc. Natl. Acad. Sci. U. S. A. 109, 10224 (2012). 16B. Jiang, J. Li, D. Q. Xie, and H. Guo, J. Chem. Phys. 138, 044704 (2013). 17B. Jiang, D. Xie, and H. Guo, Chem. Sci. 4, 503 (2013). 18P. M. Hundt, B. Jiang, M. E. van Reijzen, H. Guo, and R. D. Beck, Science 344, 504 (2014). 2M.

234705-11 19B.

A. Farjamnia and B. Jackson

Jackson and S. Nave, J. Chem. Phys. 135, 114701 (2011). Jackson and S. Nave, J. Chem. Phys. 138, 174705 (2013). 21B. Jackson, F. Nattino, and G. J. Kroes, J. Chem. Phys. 141, 054102 (2014). 22B. Jiang and H. Guo, Phys. Rev. Lett. 114, 166101 (2015). 23R. A. Marcus, J. Chem. Phys. 45, 4493 (1966). 24R. A. Marcus, J. Chem. Phys. 45, 4500 (1966). 25W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980). 26S. Nave, A. K. Tiwari, and B. Jackson, J. Phys. Chem. A 118, 9615 (2014). 27D. Han, S. Nave, and B. Jackson, J. Phys. Chem. A 117, 8651 (2013). 28H. Guo and B. Jackson, “Mode- and bond-selective chemistry on metal surfaces: The dissociative chemisorption of CHD3 on Ni(111),” J. Phys. Chem. C (in press). 29A. Sinha, M. C. Hsiao, and F. F. Crim, J. Chem. Phys. 92, 6333 (1990). 30M. J. Bronikowski, W. R. Simpson, B. Girard, and R. N. Zare, J. Chem. Phys. 95, 8647 (1991). 31G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). 32G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994). 33G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996). 34G. Kresse and J. Furthmuller, J. Comput. Mater. Sci. 6, 15 (1996). 20B.

J. Chem. Phys. 142, 234705 (2015) 35G.

Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). E. Blöchl, Phys. Rev. B 50, 17953 (1994). 37J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 38J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997). 39M. Mastromatteo and B. Jackson, J. Chem. Phys. 139, 194701 (2013). 40J. Q. Dai and J. Z. H. Zhang, J. Phys. Chem. 100, 6898 (1996). 41D. H. Zhang, Q. Wu, and J. Z. H. Zhang, J. Chem. Phys. 102, 124 (1995). 42F. Nattino, H. Ueta, H. Chadwick, M. E. van Reijzen, R. D. Beck, B. Jackson, M. C. van Hemert, and G. J. Kroes, J. Phys. Chem. Lett. 5, 1294 (2014). 43S. Nave, A. K. Tiwari, and B. Jackson, J. Chem. Phys. 132, 054705 (2010). 44A. C. Luntz and J. Harris, Surf. Sci. 258, 397 (1991). 45B. C. Garrett and D. G. Truhlar, J. Chem. Phys. 81, 309 (1984). 46C. Diaz, E. Pijper, R. A. Olsen, H. F. Busnengo, D. J. Auerbach, and G. J. Kroes, Science 326, 832 (2009). 47B. Jiang and H. Guo, J. Phys. Chem. C 117, 16127 (2013). 48X. J. Shen, A. Lozano, W. Dong, H. F. Busnengo, and X. H. Yan, Phys. Rev. Lett. 112, 046101 (2014). 49B. Jiang and H. Guo, J. Chem. Phys. 138, 234104 (2013). 50B. Jiang and H. Guo, J. Phys. Chem. C 118, 26851 (2014). 36P.